Trapped in the morphospace: The relationship between morphological integration and functional performance.

The evolution of complex morphological structures can be characterized by the interplay between different anatomical regions evolving under functional integration in response to shared selective pressures. Using the highly derived humeral morphology of talpid moles as a model, here we test whether functional performance is linked to increased levels of evolutionary integration between humerus subunits and, if so, what the strength is of the relationship. Combining two-dimensional geometric morphometrics, phylogenetic comparative methods, and functional landscape modeling, we demonstrate that the high biomechanical performance of subterranean moles' humeri is coupled with elevated levels of integration, whereas taxa with low-performance values show intermediate or low integration. Theoretical morphs occurring in high-performance areas of the functional landscape are not occupied by any species, and show a marked drop in covariation levels, suggesting the existence of a strong relationship between integration and performance in the evolution of talpid moles' humeri. We argue that the relative temporal invariance of the subterranean environment may have contributed to stabilize humeral morphology, trapping subterranean moles in a narrow region of the landscape and impeding any attempt to reposition on a new ascending slope.

The evolutionary history of vertebrates is characterized by extensive morphological changes. Such anatomical transformations are responses to altered selective regimes, spurred by changes in environmental conditions, generally coupled with the reorganization of genetic and developmental pathways (Wainwright 2007; Losos 2011; Brocklehurst and Benson 2021). Most of these transitions lead to, or were made possible by, the emergence of key morphological innovations ultimately allowing organisms to explore novel areas of morphospace (Losos 2010; Raia and Fortelius 2013; Dial et al. 2015).

A central concept in our understanding of the evolution of morphofunctional variation is how different anatomical regions of an organism integrate into a coherent, functionally apt structure (Klingenberg 2013; Orkney et al. 2021). Such integration generates covariance in anatomical structures, which is useful in maintaining functional association between traits and is not expected to be partitioned between regions (modules) evolving semi‐independently from one another.

However, fine‐tuning of a structure to a specific function can occur either way: it may require that the different parts be finely integrated toward optimal performance of a tasks, or the task be fulfilled by a single, aptly evolved module (Monteiro and Noguera 2010; Kane and Higham 2015). Understanding how different anatomical regions interconnect to achieve a particular performance outcome will provide a deeper understanding of the mechanisms that allow populations to shift toward different adaptive peaks along their evolutionary history.

Previous studies have identified an association between ecological shifts and changes in patterns and magnitudes of phenotypic integration (Martín‐Serra et al. 2015; Yuan et al. 2019; Randau et al. 2019; Rhoda et al. 2020, among others). Yet, to the best of our knowledge none of these coupled integration analyses with empirical measures of biomechanical performance. Herein we provide an example of such analysis by studying talpid moles (Talpidae, Mammalia). The mammalian family Talpidae includes ambulatorial (Uropsilini), semiaquatic (Desmanini and Condylurini), semifossorial (Urotrichini and Neurotrichini), and fully subterranean species (Scalopini and Talpini). Each of these clades display different locomotory strategies and have colonized independently a wide array of habitats: terrestrial, aquatic, and subterranean (Sánchez‐Villagra et al. 2006). Such ecological diversity has been achieved through the combination of different behavioural, physiological, and morphological adaptations, especially evident in the forelimbs, which display a clear gradient toward fossoriality along the clade's evolution (Gambaryan et al. 2002; Piras et al. 2012).

The humerus in particular is widely recognized as a strong indicator of mole ecomorphology and is tightly linked to locomotor performance (Gambaryan et al. 2002; Meier et al. 2013; Piras et al. 2015; Sansalone et al. 2020). The humerus of subterranean species displays a pattern of coordinated early ossification of the distal (muscular) and proximal (articular) regions when compared to other nonsubterranean taxa (Bickelmann et al. 2014). Furthermore, the humeri of moles are generally well‐preserved in the fossil record, allowing the inclusion of extinct taxa that are key to the accurate reconstruction of evolutionary transitions within the clade and the coverage of the realized morphospace (Ziegler 1999; Schwermann et al. 2019; Castiglione et al. 2020).

Coupling analyses of evolutionary integration with empirical measures of biomechanical performance is likely to yield a more accurate understanding of the intrinsic and extrinsic factors shaping the response of morphological structures to shifts in selective regimes (Piras et al. 2014; Maiorino et al. 2017; Evans et al. 2021).

In this study, we provide a novel framework for the investigation of the relationships between performance and morphological integration, achieved by combining two‐dimensional geometric morphometrics (GM), phylogenetic comparative methods, and theoretical morphospace modeling. The principles of quantitative trait evolution allow us to predict expected changes in the mean phenotype of a population under different selective scenarios (Fisher 1930; Wright 1932). Simpson (1944) proposed the concept of the adaptive landscape: a multidimensional space with the height of the landscape corresponding to an adaptive peak where fitness is maximized, so that a population mean phenotype is expected to move toward local peaks. Hence, the magnitude and direction of selection are defined by the shape of local adaptive peaks, whereas the rate and magnitude of phenotypic change are direct consequences of the strength of selection (Polly et al. 2016). Recently, different studies have applied the adaptive landscape concept to investigate the evolution of morphological structures within a functionally informed framework (Dumont et al. 2009, 2014; Tseng 2013; Dickson and Pierce 2019). Here, we combined this strategy with GM and Finite Elements Analysis (FEA) to measure von Mises (VM) stress, as a proxy for humeral performance, to create a form‐function landscape based on hypothetical functional morphospace. Applying this method, we tested whether functional specialization is paired with increased levels of morphological integration along the evolutionary history of the family Talpidae. By comparing integration magnitudes between realized forms and a range of theoretical phenotypes, we further provide new insights into how functional performance may shape the evolution of anatomical structures.

Material and Methods

MATERIAL

Our sample included 1015 left humeri belonging to 84 species (47 extinct and 37 extant) encompassing the known morphological variability of the family Talpidae (see Tables 1 and S3). We classified the species using the taxonomic groups defined by Sansalone et al. (2019) and partitioned talpid clades into distinct ecological categories following Sánchez‐Villagra et al. (2006) and Bickelmann et al. (2014): Uropsilini (ambulatorial), Desmanini (semiaquatic), Urotrichini and Neurotrichini (semifossorial), and Talpini and Scalopini (subterranean).

Table 1

Details on the sample and ecological categorization

Tribe	N species	N humeri	Lifestyle
Uropsilini	7	35	Ambulatorial
Desmanini	10	58	Semiaquatic
Urotrichini	10	76	Semifossorial
Neurotrichini	5	49	Semifossorial
Scalopini	22	343	Subterranean
Talpini	30	454	Subterranean

PHYLOGENETIC TREE

The phylogenetic relationships of Talpidae have been the subject of ongoing debate for decades (Whidden 2000; Shinohara et al. 2003; Bannikova et al. 2015; He et al. 2017; Schwermann et al. 2019). However, recently a consensus appears to have been reached following He et al. (2021). Because we included fossil taxa, we built our phylogenetic hypothesis based on maximum parsimony trees produced by cladistic analyses of previously published matrices of osteological characters (Sánchez‐Villagra et al. 2006; Schwermann and Thompson 2015; Hooker 2016; Sansalone et al. 2018, 2019; Schwermann et al. 2019). The matrix includes 171 total characters (five unordered; Tables S1 and S2), all characters were equally weighted, and the topology was constrained following the hypothesis proposed by He et al. (2021). Some fossil species were too poorly represented to facilitate codification of a sufficient number of characters, hence we coded characters only for the eight taxa represented by complete or near‐complete material (Domninoides mimicus, Mygalea jaegeri, Proscapanus sansaniensis, Mioscalops isodens, Mioscalops ripafodiator, Leptoscaptor robustior, Yanshuella primaeva, and Yunoscaptor scalprum) following the coding in Schwermann et al. (2019). The comparative analyses presented here were performed on a strict consensus tree generated from the four most parsimonious trees produced from a heuristic search and stepwise addition, with a random addition sequence of 1000 replicates. The cladistic analyses were performed using PAUP 4.0 a147 (Swofford 2002). The final topology was time calibrated following an extensive review of the literature including the first and last occurrence of extinct species and molecular clock estimates, when available, for extant species (see also Sansalone et al. 2018, 2019; Schwermann et al. 2019; Table S4).

SHAPE ANALYSIS

The humeri have been photographed at a fixed distance of 50 cm using a Nikon D100 camera with Micro‐Nikkor 105‐cm macro lens. We digitized 22 homologous landmarks and 14 sliding semilandmarks on homologous curves using the software tpsDig2 (Rohlf 2015), following previously published protocols (Sansalone et al. 2019) (Fig. 1a). A Generalized Procrustes Analysis (Bookstein 1986) has been performed, using the function procSym from the R package Morpho (Schlager 2014), to rotate, translate, and scale the landmark coordinates to unit centroid size (CS = the square root of the sum of squared distances of the landmarks from the centroid; Bookstein 1991). We performed a standard principal component analysis to visualize the multivariate ordination of the aligned coordinates. Shape changes associated with the PC (Principal Component) axes were visualized using a vector field between the reference (mean shape) and the target (extreme of PC axes), generated by interpolating the shape changes within the considered configuration (Claude 2008). The significance of morphological differences between clades has been assessed using a Phylogenetic Generalized Least Square approach implemented in the function PGLS_fossil from the R package RRphylo (Castiglione et al. 2018), designed for nonultrametric trees. These tests have been repeated after accounting for the effect of size measured as CS using residual shape variation.

Figure 1

(a) Landmarks and semilandmarks configuration. The dashed white line separates the two functional modules (proximal muscular region, distal articular region). (b) Finite element model constrained on the region corresponding to the clavicle articular facet and loaded on the area corresponding to the trochlea where the ulna articulates with the humerus. (c) The theoretical shapes mapped on the morphospace. (d) The reconstructed theoretical functional landscape.

HUMERAL INTEGRATION

The humerus can be divided into two functionally distinct units. The proximal region provides surfaces (pectoral crest and teres tubercle) for the attachment of the main muscles involved in locomotion (Pectoralis major, Subscapularis, Teres major, and Latissimus dorsi) and the distal part is where the loadings are applied to the trochlear area and the ulna articulates (Fig. 1b; see also Piras et al. 2012). We assessed the extent of the correlation between these two subunits using a Partial Least Squares (PLS) analysis within a phylogenetic context implemented in the function phylo.integration from the package geomorph (Adams and Otárola‐Castillo 2013). These tests were repeated after accounting for the effect of size measured as CS using residual shape variation.

We evaluated differences in the magnitude of integration between the clades and the ecological categories using a “global integration” strategy to measure the scales of shape variation, as recently proposed by Bookstein (2015). This method offers the means to obtain different measures and perspectives on the levels of integration “intrinsic” to a structure by testing the null hypothesis of “self‐similarity” or the absence of interpretable changes at any spatial scale. This analysis aims to determine whether global (large) scale or local (small) effects are the primary influences on shape variation. The framework is based on a linear regression of the logarithm of partial warps’ variance against their logarithm of the inverse of bending energy (i.e., the log of eigenvalues of the bending energy matrix computed on the mean shape; Mitteroecker et al. 2020). In this context, a regression slope larger than 1 indicates an integrated pattern at global scale (large‐scale shape changes), whereas a slope less than 1 indicates integration at local scale (small‐scale shape changes). If the regression slope is exactly 1, data can be considered “self‐similar” (see also Mitteroecker et al. [2020] for a different contextualization of the integration at local scale). Finally, we compared the resulting slopes between subterranean and nonsubterranean species using the R package RRPP (Collyer and Adams 2018).

RATES OF MORPHOLOGICAL EVOLUTION

We applied the recently developed RRphylo framework (Castiglione et al. 2018) to compute the rates of phenotypic change along each branch of the phylogeny. The RRphylo strategy is based on the phylogenetic implementation of ridge regression to estimate rates values at nodes and along branches. After each value has been assigned, the function search.shift can be used to locate evolutionary rates that are significantly lower or higher than the tree average by means of randomization. To account for the effect of size, we repeated the RRphylo analysis using the residuals of multivariate regression between Procrustes aligned coordinates and the CS.

To account for uncertainty in tree topology, branch lengths, and node heights, we used the function overfitRR. This randomly swaps species positions on the tree (15% of the total tips) and removes a number of tips corresponding to the 30% of the tree size. In addition, the function “moves” in time 15% of the tree nodes by changing the age of the node in between that of its direct ancestor and the age of its oldest descendant. At each overfitRR iteration, the function search.shift is applied on the pruned tree and the new pattern is compared with the original data to assess the robustness of the results. We performed 100 overfitRR iterations.

Lastly, to test for statistically significant differences among the evolutionary rates of the two modules, a Friedman rank sum test, a nonparametric version of a repeated measures ANOVA, was carried out on the logarithms of the evolutionary rates identified with the RRphylo strategy. As a post hoc test, a Wilcoxon signed‐rank test was performed. Both these tests were implemented in the package rstatix (Kassambara 2020) and significance was assessed following 1000 permutations.

FUNCTIONAL PERFORMANCE

Locomotor performance of the theoretical and realized shapes was predicted by means of FEA (Fig. 1b). We applied FEA in a comparative fashion removing size information scaling all models to the same surface area (SA), in keeping with the recommendations for biomechanical shape analyses of Dumont et al. (2009).

The models were treated as homogeneous solids and subjected to the same loads. The same material properties have been attributed to all specimens, that is, isotropic material properties corresponding to Haversian bone (Young modulus: 10 GPa, Poisson's ratio: 0.41; Rayfield et al. 2001), although in this comparative context the material properties themselves will not affect interpretations provided that they are uniform across the sample. We used the inverse of VM stress as a proxy for performance. VM stress was averaged across the entire humeral shaft after removing the top 2% values to account for artefactually high values at the boundary and on the site of force application (Dumont et al. 2014).

The three‐dimensional volume meshes were generated by extruding the two‐dimensional landmark configurations along the orthogonal plane. Although this represents a simplification of the humeral shape, this approach allows for a faster solution of the Finite Element Models, the inclusion of fossil species, and simpler comparison with the theoretical shapes. It has been demonstrated that three‐dimensional extruded models are a viable alternative to three‐dimensional models generated from CT scans (Morales‐García et al. 2019). However, it must be stressed that two‐dimensional models, or homogeneously three‐dimensional extruded models, are limited in capturing the three‐dimensionality and the muscle arrangement of the structure as well as the forces acting on it. These models necessarily represent a first approximation due to their simplicity and further studies should be aimed at incorporating more complex models to allow finer simulation of structures’ biomechanical behavior.

We also note that previous FEA analyses based on two‐dimensional and three‐dimensional reconstruction of talpid humeral meshes have revealed largely similar interspecific pattern of stress distribution (with subterranean species characterized by less stress than nonsubterranean taxa), demonstrating that extruded surfaces can be good predictors of the humeral biomechanical behavior, at least in a comparative context (Piras et al. 2012, 2015; Sansalone et al. 2020).

As maximum thickness values measured along the humeral shaft are only available for the observed sample, thickness of the theoretical meshes was determined by using the coefficients derived from the regression between thickness (dependent) and maximum width of the proximal region of the humerus (independent), following the protocol defined by Piras et al. (2012). The models were fully restrained in all directions on the area corresponding to the clavicular facet (articulation with the clavicle) and on the pectoral ridge where the pectoral muscles perform both stabilizing and propulsive functions and an elastic boundary condition (modeled using a spring constant) has been applied (Gambaryan et al. 2002; Piras et al. 2012, 2015). To evaluate the VM stress distributions on the humerus, we applied the same loads (22 N; Gambaryan et al. 2002) on the area corresponding to the trochlear region in the orthogonal direction, simulating the direction of the force applied during locomotion (Scott and Richardson 2005; Piras et al. 2012).

THEORETICAL MORPHOSPACE

The theoretical morphospace was produced by sampling 80 warps from a 10 × 8 grid from a space defined by the first two PC scores as they maximize morphological variability (89% of total variation). The landmark configurations corresponding to each of the 80 warps have been reconstructed using the function showPC from the package Morpho (Fig. 1c).

MORPHOFUNCTIONAL LANDSCAPE

A functional surface for locomotion performance was constructed by fitting a thin plate spline (TPS) regression using restricted maximum likelihood to estimate the smoothing parameters (Silverman and Green 1994). We used the Tps function from the fields package (Nychka et al. 2015) using the first two PC scores as x and y axes and the VM stress values measured on the theoretical warps as z‐axis (Fig. 2b). As in an adaptive landscape, where adaptive peaks and valleys are represented on the third axis by fitness, the morphofunctional landscape can map the functional performance measured from biomechanical simulations over a bivariate plot representing the phenotype (shape) under investigation (McGhee 1999; Tseng 2013). Consequently, we mapped the realized shapes on the previously generated functional landscape to visualize landscape occupation in relationship to areas of maximum and minimum performance. Finally, to determine the relationships between functional performance and morphological integration, we randomly sampled 45 theoretical shapes occupying areas of performance values that exceeded those displayed by the realized shapes and compared the relative magnitudes of integration. The number of selected theoretical shapes was chosen to avoid having a sample too different from the largest clade (30).

Figure 2

Phylomorphospace of the actual species based on PC1 and PC2. Deformations are visualized through field of vectors and reflect mean clade shape changes along axes of variation. Color scale is indicative of shape deformation (expressed in distance in mm). Warmer colors indicate larger deformations; cooler colors indicate smaller deformations.

It must of course be borne in mind that locomotion kinematics constitutes only part of an organism's functional repertoire and not its entire fitness. However, for lineages that have experienced a directional trend toward increased specialization to perform a particular task, such as talpids, exploring a specific function helps to conceptualize the evolution of relevant traits.

Results

The recovered topology reflects the hypothesis published by He et al. (2021), specifically Uropsilini being basal to the other talpid clades, whereas the North‐American subterranean moles (Scalopini) are sister to the remaining clades. Desmans are recovered as sister of a clade composed of the semifossorial shrew‐moles (Neurotrichini and Urotrichini) and the Eurasian subterranean taxa (Talpini).

SHAPE VARIATION

The phylomorphospace based on the first two PC scores (89% of total variation) evidenced the separation of subterranean taxa (Talpini and Scalopini) from the other clades, with a clear shift from a slender and elongated humeral shape to a robust and more rounded one, typical of the highly specialized mole species (Fig. 2). Along the second axis, the ambulatorial taxa (Uropsilini) are distinguished from the semiaquatic (Desmanini) and semifossorial groups (Urotrichini and Neurotrichini) along a positive to negative gradient.

MORPHOLOGICAL INTEGRATION

Results from PLS and global integration analyses are summarized in Table 2. Overall, we found the two humeral subunits (muscular and articular) to be tightly integrated (r‐PLS = 0.910; slope = 1.03; P = 0.001) and similar levels of integration were also found when accounting for size (r‐PLS = 0.904; slope = 1.03; P = 0.001). However, subterranean tribes (Talpini and Scalopini) displayed significantly higher levels of integration compared to the other, less specialized taxa. The semifossorial Neurotrichini were characterized by a lower slope (although not statistically significant) than the other nonsubterranean clades (Uropsilini, Desmanini, and Urotrichini).

Table 2

Results from phylogenetically corrected PLS and global integration analyses

PLS pairwise comparisons
	Uropsilini	Desmanini	Urotrichini	Neurotrichini	Scalopini	Talpini
Uropsilini	1	0.728	0.863	0.675	0.03	0.018
Desmanini		1	0.504	0.93	0.002	0.004
Urotrichini			1	0.434	0.018	0.008
Neurotrichini			1	0.001	0.001
Scalopini					1	0.977
Talpini						1

Note: Significant P‐values are in bold.

EVOLUTIONARY RATES

Overall, the subterranean taxa displayed significantly higher rates of morphological evolution, whereas the semifossorial groups exhibited the lowest (Fig. 3a). The two humeral modules were characterized by significantly different evolutionary rates, with the muscular region showing faster evolution than the articular subunits (Fig. 3b).

Figure 3

(a) Distribution of evolutionary rates along the talpid phylogeny. Branches are colored according to the magnitude of the evolutionary rates. (b) Boxplot representing the difference in rates of evolution between the proximal muscular region (faster, blue box) and the distal articular region (slower, orange box).

We found that our results are not due to the tree topology adopted after applying overfitRR to account for uncertainty in tree topology, branch lengths, and node height; we found similar results, with rate shifts still located at the same nodes. These results also remained consistent when size (CS) was included as a covariate.

The morphofunctional landscape identified predictable trends in space occupation. Fully subterranean taxa (Talpini and Scalopini) were tightly clustered in a relatively small area of the landscape characterized by high biomechanical performance within the isocline defined by fitness values between 0.9 and 1 (Fig. 3a,b). Semifossorial species (Urotrichini and Neurotrichini) showed a broader adaptive landscape occupation pattern, ranging from areas marked by low‐performance (isocline 0.4–0.5) to areas characterized by high‐performance values (isocline 0.8–1). Lastly, ambulatorial and semiaquatic species (Uropsilini and Desmanini) occupied areas characterized by low‐performance values (isocline 0.1–0.4) (Fig. 4a,b).

Figure 4

(a) Two‐dimensional visualization of the theoretical functional landscape. The circles represent the actual species mapped over the functional landscape; the triangles indicate the theoretical morphs sampled within isoclines characterized by higher performance values. Morphological deformations of the theoretical models are visualized through field of vectors and reflect mean variation of each of the three considered isoclines versus the consensus. (b) Three‐dimensional visualization of the theoretical functional landscape.

The areas characterized by the highest performance values (three isoclines ranging from 1 to 1.3) are not occupied by any taxon. These theoretical shapes are characterized by extreme enlargement of the pectoral ridges, teres tubercles, and clavicle articular facets, whereas the distal region of the humerus did show relevant changes (Fig. 4a). However, the theoretical shapes displayed a significant drop in the levels of integration compared to those of the realized taxa showing higher performance values (Fig. 5). In contrast to the shapes of actual species, the theoretical morphs showed an inverse trend between performance and integration levels, with the latter showing integration levels (measured within three different isoclines) that were also lower than the values displayed by the nonsubterranean taxa (Fig. 5).

Figure 5

Barplot indicating the values of performance (von Mises stress) and integration (global integration slope, positive values are reported for the sake of visualization) measured for each clade of the actual species and for each isocline including the theoretical morphologies.

Discussion

INCREASE IN FUNCTIONAL PERFORMANCE IS PAIRED WITH INCREASE IN MORPHOLOGICAL INTEGRATION

The humeri of moles show some of the most extreme morphological modifications known among vertebrates (Barnoski 1982), offering a remarkable example of how selection can reorganize the patterns of correlated changes of functionally relevant traits. Talpids occupy a wide range of environments, each requiring different sets of adaptations, with those of the subterranean taxa being the most demanding from a functional perspective (Nevo 1979; 1995).

Our analyses revealed that the increase in functional performance in the clade was paralleled by increased degree of integration throughout its evolutionary history (Fig. 5).

The higher degree of integration found in subterranean species supports the prediction that trait integration played a key role in maintaining a viable functional association within anatomical structures, whereas the lower degree of correlation detected among the nonsubterranean clades might imply a more independent evolution for the two humeral subunits, likely experiencing relaxed functional demands. In this context, our results support the view that natural selection can alter the pattern and magnitude of integration at the evolutionary level in parallel with functional performance (Martín‐Serra et al. 2015; Vidal‐Garcia and Keogh 2017; Sansalone et al. 2019).

The humeri of subterranean species are highly adapted for digging and high functional performance might have been acquired by coupling the expansion of the proximal region of the humerus, generating propulsion, with the expansion of the distal region to sustain the increased loadings (Piras et al. 2015; Sansalone et al. 2020). Moreover, the expansion and ossification of the two humeral subunits occur early during development (Bickelmann et al. 2014) and in a coordinated fashion, supporting the importance of integration in maintaining functional association between the parts. However, despite their functional relevance, the two humeral subunits evolved at different paces, with the muscular region characterized by higher evolutionary rates (Fig. 3b). These results support the conclusion that muscle attachment regions better reflect ecological differences, being the primary locus of wider morphological modifications from an evolutionary perspective (Young and Badayev 2006).

WHAT IS THE RELATIONSHIP BETWEEN MORPHOLOGICAL INTEGRATION AND FUNCTIONAL PERFORMANCE?

If performing complex tasks requires a tighter coordination between interacting parts, then it is reasonable to expect that a positive shift in functional performance should be followed by an increase in the degree of morphological integration at the evolutionary scale. In other words, higher correlation between traits may be required to allow a population to move upward along an adaptive landscape. However, contrary to this expectation, our experiment using the theoretical morphospace supported the opposing trend. Specifically, landscape regions characterized by higher performance values showed a marked drop in the integration level. Interestingly, the slopes measured within three different isoclines (see Figs. 4, 5) displayed values close to or lower than the critical value of 0.56 identified by Bookstein (2015) and Windhager et al. (2017) for integration at local scales (please note that we inverted the value to agree with the most recent proposition of this concept; see also Mitteroecker et al. 2020). To date, no biological samples have displayed regression values lower than 0.56 suggesting that biological structures should display some level of integration to perform adequately, and that complete traits independence in evolutionary terms is not feasible in the real world (Bookstein 2015; Windhager et al. 2017). The theoretical morphs reconstructed within the three isoclines characterized by extremely high‐performance values showed values of 0.58, 0.22, and 0.14 indicating that further increases in functional performance would generate highly nonintegrated geometries that possibly cannot be achieved by living organisms.

Here, different factors can be invoked to explain this pattern. At first, the theoretical humeral shapes showed an extreme expansion of the muscular area (proximal subunit) without correspondingly substantial changes in the articular area (distal subunit), likely generating the loss of correlation within the humerus. We can speculate that the enlargement of the muscular area could allow more room for the insertion site of the latissimus dorsi and teres major muscles, two of the main muscles involved in digging (Gambaryan et al. 2002). On the one hand, the expansion of the teres tubercle allows the attachment of larger muscles generating higher forces, but on the other, this could trigger an energetic imbalance. Foraging activity by subterranean moles is particularly demanding, necessitating that these animals daily consume a biomass equivalent to one third of their body weight (Gorman and Stone 1990; Bertolucci et al. 1999). Second, the high‐performance region may be inaccessible given the previously noted developmental pressures (Bickelmann et al. 2014). The early coordinated ossification of the two humeral subunits may result in reduced evolutionary potential. A similar pattern is shown by marsupials where intense functional demand for continuous suckling earlier in their development could drive high integration, resulting in higher covariation and limited variability of the marsupial oral apparatus (Goswami et al. 2016).

Finally, the high stability and predictability over time of the subterranean environment may have inhibited any ratchet‐like evolution of the phenotype (Tseng 2013; De Vos et al. 2015), meaning that the lack of substantial environmental changes over time may trap organisms in a local optimum rendering them unable to reposition on a new ascending slope.

In this context, the region marked by extremely high‐performance values may be considered “maladaptive” (sensu Simpson 1944) so that any upward transition over the functional landscape might be not possible for talpid moles. Thus, it can be hypothesized that subterranean moles are experiencing an evolutionary stasis around a stationary adaptive peak (Strathman 1978).

As noted previously, the proxy for performance used in this study may not represent a complete measure for fitness. However, among extant talpids, subterranean clades are by far the most specious and have a wider distribution than nonsubterranean taxa (many of these have now relict distributions; Ziegler 1999; Shinohara et al. 2003; Schwermann et al. 2019); this suggests that the increased functional performance necessary to exploit the subterranean environment could have played an important role in the evolution of the family.

Conclusions

Our findings demonstrate that combining biomechanical and integration analyses of theoretical models within a functional landscape can yield a more complete and accurate understanding of how natural selection could affect the patterns and degrees of integration at evolutionary scale and how these changes can be associated with higher functional performance.

More complex models and frameworks aimed at including multiple traits and functions linked by a trade‐off relationship show potential in expanding our knowledge of the mechanisms underlying the relationship between form and function (Dickson et al. 2021).

AUTHOR CONTRIBUTIONS

GS designed the study. GS, PR, and SW wrote the manuscript with contribution from all the authors. GS, PC, and RC assembled the sample. GS, SC, PR, and PC analyzed the data.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

DATA ARCHIVING

Averaged per species landmark data and phylogeny are available in the data repository Dryad (https://doi.org/10.5061/dryad.1c59zw3wh).

Associate Editor: Dr. T. Stayton

Handling Editor: Dr. M. L. Zelditch

Tables S1, S2, S3 and S4.

Click here for additional data file.