The effect of landscape on functional connectivity and shell shape in the land snail Humboldtiana durangoensis

Resistance surfaces

The geographic centroids of each one of the seven genetic groups of H. durangoensis previously defined by microsatellite loci by López, Gómez & Mejía (2017) were used to determine the effect of the landscape on functional connectivity (Fig. 1). Whereas the landscape can include a large number of variables, in the present work, we followed two approximations to evaluate the functional connectivity between snail populations. The first was to use an approximation of the Grinnellian niche defined from a set of bioclimatic variables (Bell et al., 2010; Ortego et al., 2012; Poelchau & Hamrick, 2012); the second was to analyze the effect of vegetation cover, because it is known that it affects the dispersion of terrestrial snails (Labaune & Magnin, 2002; Armbruster, Hofer & Baur, 2007; Ström, Hylander & Dynesius, 2009; Edworthy et al., 2012; Kappes et al., 2009), especially in mountain populations where periods of glaciation and deglaciation promoted the contraction and expansion of vegetation cover (Armbruster, Hofer & Baur, 2007). In both cases, the different models were generated for three different time frames, including the current period and two time periods representing the extreme conditions experienced during the late Quaternary: the middle Holocene (6,000 years bp), which was warmer and wetter than the present, and the last glacier maximum (LGM), which was characterized by dry and colder climates (21,000 years bp).

To reduce the error in the parameterization, validation and comparison of the models (Barve et al., 2011), the available geographic space for the taxon (M) was defined as the Ecoregion Madrense Centro (González-Elizondo et al., 2013). Grinnellian niche models were constructed with the 19 climatic variables available in WorldClim (Hijmans et al., 2005) and the 18 climatic and topographic variables available in ENVIREM (Title & Bemmels, 2018). The models were made at a resolution of 30 arc-seconds, but in the case of the LGM, the variables were used at their native resolution of 2.5 min, and a bilinear interpolation was performed to decrease the resolution to 30 arc-seconds with the disaggregate function of the raster library ver. 2.6-7 in R (Hijmans, 2017). The atmospheric circulation model used was the MPI-ESM-P, since it has shown better performance with respect to other models of circulation (Tang et al., 2017). The bioclimatic variables were clipped to the geographic space with the crop and mask functions of the raster library ver. 2.6-7 in R (Hijmans, 2017).

Species niche model

The environmental suitability areas were defined by a maximum entropy algorithm (MAXENT v. 3.2.19, Phillips, Anderson & Schapire, 2006) from 28 records of H. durangoensis available in museums and our own collections. We selected this algorithm because it produces reliable results even with a small quantity of data (Elith et al., 2006; Heikkinen et al., 2006; Hernandez et al., 2006). In a preliminary analysis, the 19 WorldClim and 18 ENVIREM variables were included with the default parameters and log output to minimize the correlation and maximize their contributions to the model. The relative importance of each variable was determined from its percentage of contribution and for the loss of predictive power when each variable was excluded using a jackknife test. In addition, to select those variables with correlation coefficients lower than 0.6, environmental information was extracted from each geographic point, and a Pearson correlation test was performed with the function correlation test in the psych library of R (Revelle, 2018). Thus, the geographic distribution model was obtained with the selected variables and assumed 10,000 pseudoabsence points separated by 1 km from the presence records (Barbet-Massin et al., 2012). The statistical evaluation of the model was carried out in 10 repetitions and the data were partitioned into 75% for training and 25% for evaluation with a logistical output. The predictive power of the model was evaluated using a partial ROC test with 100 bootstrap replicates (Barve, 2008). Finally, the suitability area available for the species in each temporal frame was estimated with the DEM surface tools in ArcGIS 10.

Vegetation models

The random forest (RF) classification algorithm was used to obtain the modeled vegetation cover (Breiman, 2001). This method categorizes a set of data based on the classification and regression of the trees from a bootstrap analysis (Breiman, 2001). The INE-INEGI (1997) vegetation cover map was used as an input file. Because this classification contains many vegetation types for the Madrense Centro ecoregion, prior to the analysis, the vegetation types were reclassified into five categories based on the ombrothermal horizons of the Sierra Madre Occidental (Macías-Rodríguez, Giménez de Azcárate-Cornide & Gopar-Merino, 2017): (1) temperate forests, (2) cold temperate forests, (3) grasslands, (4) tropical forests and (5) drylands. The model was trained to take into account the variables sets of BIOCLIM and ENVIREM and to select only those that explained more than 50% of the variation based on the mean decrease accuracy criterion; these models were made with the RF 4.6-14 library in R (Breiman, 2001).

Isolation by resistance

To evaluate the resistance of the landscape between the genetic groups, the resistance isolation model (IBR) was implemented in CIRCUITSCAPE 3.4.2 (McRae, 2006). This method produces a resistance/conductance matrix between the pairs of sites that are obtained by assigning an arbitrary resistance/conductance value per pixel corresponding to the relative resistance of the landscape to the genetic flow. The result was a resistance value that depended on the distance between the localities, the number of possible pathways and the heterogeneity of the landscape (McRae, 2006). The following resistance values were assigned to the forest structure: 60 (cold temperate forests), 110 (temperate forests), 200 (grasslands), 300 (tropical forests) and 360 (drylands).

For the surface derived from the niche modeling, resistance values were assigned considering five symmetrical categories defined by the range between the minimum training presence and the highest suitability value obtained by the Maxent model. The values were assigned with the ifelse and raster R libraries (Hijmans et al., 2005), and three different approaches were used to evaluate the relationship between the paired F_ST values among the seven genetic groups (López, Gómez & Mejía, 2017) and resistance values. Three matrices were considered in this analysis: the genetic paired distances, the Log10-transformed Euclidean geographical distances, and the paired resistance distances obtained from CIRCUITSCAPE for the two evaluated resistance surfaces (climate and vegetation). In the first approximation, the Mantel partial test was used to evaluate the effects of the two variables while controlling for the effect of a third. The significance of the partial correlation of the Mantel test was obtained by 1,000 random permutations using the partial mantel test function of the NFC library (Bjørnstad, 2013). In the second approach, a distance-based redundancy analysis (dbRDA) was used in the vegan 2.5 library (Oksanen et al., 2013) considering the genetic distances, geographic distances, and the effect of vegetation cover, as well as the effect of the climatic distances on the mean of the resistance values (Noguerales, Cordero & Ortego, 2016). The characterization of the environmental space was performed with the random Points function in R that generated 1,000 random geographic points and with the extract function to obtain the climatic point value per site. Then, the main function in R was used to perform the PCA, and later, the dist function in R was used to obtain the eigenvalues of the environmental distances for the first three components considering only the loadings of the geographic points corresponding to the genetic groups; finally, the significance of the dbRDA was evaluated with the anova.cca function in R. Lastly, in the third approach, given that in the two previous analyses the climate component was not significant (see Tables 1 and 2), we only evaluate the effects of geography (G), vegetation cover (E) and both (G + E) on functional connectivity through a Bayesian framework implemented in the SUNDER 0.0.4 library (Botta et al., 2015). The algorithm implemented in SUNDER assumed that the covariance of the allelic frequencies among the populations would decrease as a function of the geographical and environmental distances (Botta et al., 2015). Thus, to estimate the effect of the set of G, E, and G + E variables, 10 independent chains with 10⁷ iterations and sampling every 1,000 steps were used with uniform priors with large upper bounds (Botta et al., 2015).

Shell morphometrics

A total of 129 shells of H. durangoensis adults from the seven genetic groups used by López, Gómez & Mejía (2017) were analyzed: Las Peñas (20 shells), El Salto (8 shells), Progreso (3 shells), Topia (7 shells), Potrero (46 shells), Los Herreras (25 shells) and Guanaceví (20 shells). The shape of the shell was obtained using two approaches: a classical approach that assumed four linear shell measurements (height, SH; width, SW; aperture height, AH; maximum aperture width, AW) obtained with a digital micrometer with an accuracy of 0.01 mm; in addition, globosity (G = SH/WD), spiral height (SP = SH−AH) and shell volume (V) were calculated (Fig. S1). These variables were Log10 transformed to remove the size effect following the method described by Mosimann (1970). Finally, the eigenvalues of the mean and the centroid values for each one of the genetic groups were recovered from a principal component analysis (PCA) for posterior analysis (Hartigan & Wong, 1979). On the other hand, the shell shape was evaluated from 11 landmarks according to Mumladze, Tarkhnishvili & Murtskhvaladze (2013) (Fig. S1). Following the method proposed by Kistner & Dybdahl (2013), a total of five photos were taken per individual to eliminate the error associated with the orientation. The X/Y coordinates were digitized in TPSDIG ver 2.12 (Rohlf, 2008). The average shape per genetic group and the deformation grids were obtained from a generalized Procrustes analysis in order to visualize changes in the shape of the shell with the gpagen function implemented in geomorph 3.0.7 (Adams, Collyer & Kaliontzopoulou, 2019). To eliminate the phylogenetic effect on the variation in shell shape, a phylogenetic principal component analysis was performed considering the shell shape of both, classical and geometric morphometrics approaches and a tree based on distances generated from F_ST values with the function phyl.pca in phytools (Revell, 2012); additionally, the deformation grids of the average shell shape of the genetic groups were ploted in the phylomorphospace with the plotGMPhyloMorphoSpace function implemented in Geomorph. Lastly, to determine whether there was a relationship between the shell shape and environmental conditions, a redundancy analysis (RDA) was performed considering the three matrices generated (means and centroid values for each one of the seven genetic groups as well as the average shape obtained from the geometric morphometrics analysis) with the rda function following the method proposed by Borcard, Gillet & Legendre (2018) in the vegan library ver 2.5 (Oksanen et al., 2013).

Environmental suitability and vegetation models

The results obtained for the modeling of the distribution area of H. durangoensis in the Madrense Centro region showed that the models constructed for the three temporal frames were satisfactory (P = 0). In general, our findings suggested that the areas of environmental suitability had decreased considerably in the last 21,000 years (38,197 km² or 28.5% of the total area in the LGM, 32,945 km² or 24.5% in the mid Holocene and 23,620 km² or 17.6% in the current). Our findings show that at present, the areas with high probability of occurrence are restricted to the northern portion of the distribution area (Fig. 1). The model of vegetation cover generated from the current vegetation map with RF showed that the estimated success rate was 76.77% for the LGM, 77.48% for the middle Holocene and 75.57% for the current period (Table S1). In addition, in the last 21,000 years, a variation in the coverage area of each plant community was estimated, and the temperate forests increased the most, while the grasslands decreased the most (Table S1; Fig. S2).

Resistance and functional connectivity

The maps generated by CIRCUITSCAPE considering the structure of the vegetation cover suggested that the connectivity routes between the H. durangoensis genetic groups in the Central Madrense region have changed little in the last 21,000 years, although in the actual period, the areas of high resistance are larger compared to those in the LGM (Table 1; Fig. S2). The resistance surface from the environmental suitability models for the Mantel test and Mantel partial tests were not significant (Table 2; Fig. S3). On the other hand, when considering the effects of vegetation cover, the Mantel test between the values of F_ST and vegetation cover was once again not significant in any of the three time frames; however, the Mantel partial tests yielded significant correlations when controlling for the effects of geography and vegetation cover in the three time periods (Table 2). In the case of the RDA, the marginal tests for the three time frames showed a significant association between the genetic differentiation and geographic distance, explaining 24.26% of the variance, but were not significant when the resistance distances generated from the vegetation cover or from the climatic variables were considered (Table 3). In contrast, in the conditional tests as in the Mantel partial test, a relationship was again observed with the structure of the vegetation cover but not with that of the climate (Table 3). With respect to the results generated by SUNDER, when the climatic component was no longer considered, it was observed that during the LGM, it was the geographic component that best explained the variation, while for the middle Holocene and the actual period, both the geographic component and the vegetation cover were important (Table 4).

Variation in shell size and shape

The values estimated from the morphometrics classical approach allow us to establish that the populations located in the north of the distribution area (Topia, Potrero, Los Herreras and Guanaceví) had larger sizes and higher spires in comparison with the populations in the center (Progreso) and south (Las Peñas and El Salto) of the distribution area (Table 5). The percentage of variance explained by the first three phylogenetic components was 99.47% for the means of the linear variables, 99.98% for the centroid size and 91.84% for the average shape obtained from the analysis of the geometric morphometrics. Finally, the RDA obtained from the analysis of the first three phylogenetic components was statistically significant (P < 0.05). The bioclimatic variables associated with each dataset were different, although in all cases, they were exclusively temperature variables, with the temperature annual range (Bio 7) being the only common variable (Fig. 2). Although it was difficult to establish a pattern, the data retrieved from geometric morphometric analysis allowed us to suggest that larger shells with higher spirals are related to the max temperature of the warmest month (Bio 5), while smaller shells with the lower spirals were related to the temperature annual range (Bio 7) and mean temperature of the wettest quarter (Bio 8) (Fig. 2). Additional support for the aforementioned results proceed from the analysis of the deformation grids in the phylomorphospace, where the populations located to the North (Guanaceví, Los Herreras, Topia and Potrero) tend to have higher spires and higher values of whorl expansion ratio that lead to squared shells, in contrast, the lower spires and the lower values of whorl expansion ratio in the Center (Progreso) and South (El Salto and Las Peñas) populations lead to wider an more rounded shells (Fig 3).

Effects of the landscape on functional connectivity

The functional connectivity in terrestrial snails was determined by the availability of microhabitats suitable for dispersal. Our findings showed that the variables related to the humidity and relative aridity of the terrain, as well as the precipitation of the driest and the wettest month, had a greater contribution to the potential distribution model generated by Maxent. These variables were related to the apparent rupture of the estivation period in May and to the period of activity and dispersion between July and September, as has been suggested for other members of the group (Baur, 1986; Aubry et al., 2006). However, the climate component defined through the environmental suitability analysis with the MAXENT maximum entropy algorithm and by the method proposed by Noguerales, Cordero & Ortego (2016) did not contribute significantly to explaining the functional connectivity of H. durangoensis populations.

A possible explanation for this phenomenon might be related to the spatial resolution provided by the bioclimatic layers. It has been demonstrated that the geographic patterns of the areas of environmental suitability in the terrestrial mollusks were particularly dependent on the resolution of the grid, since this increases or diminishes the heterogeneity of the geographic space (Kadmon & Heller, 1998). However, the models generated for land snails at a resolution of 30 arc-seconds (1 km²), as used in this study, have been shown to be efficient in explaining the historical demographic reductions that are the consequence of contractions in the areas of environmental suitability (Horsák et al., 2010; Pfenninger et al., 2014; Mumladze, 2014; Patrão et al., 2015). In this sense, the areas of environmental suitability for H. durangoensis have decreased from 28.5% in the LGM to 17.6% at the present, a result congruent with the population reductions recovered for this species with microsatellite markers and DNA sequences (López, Gómez & Mejía, 2017; López, Zúñiga & Mejía, 2019). Therefore, although the climate component apparently did not make a significant contribution to functional connectivity, its influence on the taxon cannot be denied because H. durangoensis likely experienced environmental tracking as a consequence of climate change, as has been demonstrated in alpine populations of Arianta arbustorum (Baur & Baur, 2013).

On the other hand, the RF algorithm has been shown to perform well in predicting the current vegetation types in heterogeneous geographic areas, as it was very robust in relation to the number of classes in which plant communities were clustered, as has been verified by paleopalinological records for models generated for LGM (Waske & Braun, 2009; Rodriguez-Galiano et al., 2012; Vanselow & Samimi, 2014; Hais et al., 2015). Thus, the results of efficiency in the assignment to plant categories with the RF algorithm (Table S1) fall within the values obtained in other works (Waske & Braun, 2009; Hais et al., 2015), suggesting that predictions of vegetation cover in this study are correct. Although our paleovegetation maps apparently did not show significant changes in vegetation cover (Fig. S2), the resistance results from Circuitscape suggested that these changes have occurred and that resistance values have increased from the LGM to the present (Table 1). One of the main limitations of analyses based on resistance surfaces is that the values assigned to each of the categories are arbitrary; however, it has been shown that the assigned resistance values have no effect on the habitat categories in a fragmented landscape (Schweiger, Frenzel & Durka, 2004; Wang et al., 2008). Consequently, as has been reported for other mountain snails (Scheel & Hausdorf, 2012; Hugall et al., 2002; Sherpa et al., 2018), the altitudinal displacement of plant communities in mountainous regions during Quaternary climate changes could explain the dynamics of functional connectivity in H. durangoensis as has been postulated for other species distributed in the SMOc (Metcalfe et al., 2000; Anducho-Reyes et al., 2008; Bryson et al., 2011; López-González, Correa-Ramírez & García-Mendoza, 2014).

Based on these findings, we hypothesized that the functional connectivity of H. durangoensis on different temporal scales has been promoted by the presence of both temperate and cold temperate forests and that two patterns can be distinguished as has been suggested in Helix aspersa and Cepaea nemoralis (Arnaud, 2003; Schweiger, Frenzel & Durka, 2004; Barahona-Segovia et al., 2019). The first is a model of isolation by distance on a larger geographic scale (Pfenninger & Posada, 2002; Arnaud, 2003; Schweiger, Frenzel & Durka, 2004), and the second is possible dynamic metapopulation promoted both by environmental and landscape heterogeneity on a fine geographic scale, as has been documented for other land snails (Arnaud et al., 2001; Baur & Baur, 2013).

Variation of the shell in H. durangoensis

The relationship between shell size and shape in land snails with climatic variables of temperature and precipitation has been widely studied and is well known (see review in Goodfriend (1986)). However, while the effect of the genetic component on shell shape variation has been studied (Goodacre, 2001; Dowle et al., 2015; Sherpa et al., 2018), few studies have attempted to control this effect (Webster, Van Dooren & Schilthuizen, 2012; Kotsakiozi et al., 2013), and none so far have evaluated this effect at the intraspecific level. Our findings showed, after controlling for the genetic effects, that the shell size and shape were determined by climatic variables of temperature and precipitation (Fig. 2). However, whereas these variables were not significant to explain the genetic relationships among the groups, they suggested that both the phenotype and genotype were the results of independent processes (Haase & Misof, 2009); that is, the microhabitat conditions had a great effect on the shell despite the existence of gene flow (Chiba & Davison, 2007; Fiorentino et al., 2013; Stankowski, 2013; Proćków, Kuźnik-Kowalska & Mackiewicz, 2017). Thus, whereas it has been suggested that the use of comparative phylogenetic methods at intrapopulation levels may generate poor informative results (Niewiarowski, Angilletta & Leaché, 2004), the power of resolution of these methods may depend on the taxon and the assessed trait (Martins & Housworth, 2002), as has been found in this study.

In addition, our results suggested that populations with larger shells and apertures are distributed to the north, while populations with smaller shells and apertures were distributed to the south. The altitudinal interval of the sampled localities in the northern region (1,702–2,400 m asl) was lower than the altitudinal interval in which the populations in the southern region were collected (2,587–2,759 m asl), which was consistent with the results previously found in intrapopulation studies of the species of the genera Arianta, Vestia and Trochulus (Buria & Stahel, 1983; Baur & Raboud, 1988; Sulikowska-Drozd, 2001; Proćków, Kuźnik-Kowalska & Mackiewicz, 2017), where the populations from colder climates had smaller shells. This could be related to a greater probability of survival of organisms with small shells in unfavorable climatic conditions (Baur et al., 2014) and the greater resistance to crystallization temperatures (Ansart et al., 2014). At the same time, at higher altitudes, the duration of individual growth time is shorter (Anderson, Weaver & Guralnick, 2007; Proćków, Kuźnik-Kowalska & Mackiewicz, 2017). However, there were also differences in the sizes of the aperture and the heights of the spires between the north and south regions. These shells attributes could reflect microclimatic conditions, where small apertures tended to occur in the drier and higher altitude regions, meanwhile large apertures and higher spires occured at lower altitudes as has been reported in other species (Anderson, Weaver & Guralnick, 2007; Haase & Misof, 2009; Dowle et al., 2015).

How many ESUs?

In the literature, only two published works that addressed the definition of the ESUs of land snails have been published (Holland & Hadfield, 2002; Ursenbacher et al., 2010); however, they did not remove the phylogenetic effects, which impacted their results. In the first study, a fragment of the mtCOI DNA was used and only the phylogenetic trees, genetic distances and AMOVA analysis in the 12 populations of the tree snail Achatinella mustelina were recognized as six ESUs that were reproductively isolated and distributed throughout a longitudinal transect of 24 km (Holland & Hadfield, 2002). In the second study, which used microsatellite loci and performed a genetic structure analysis, two main clusters were found in Trochulus aureatus, although the authors decided to define each one of the nine sampled populations as different ESUs, even though they were separated by less than 200 m (Ursenbacher et al., 2010). In opposition, our results suggest that each of the seven genetic groups previously identified by the analysis performed by López, Gómez & Mejía (2017) must be considered an ESU, not only because of their genetic distinctiveness but also due to the phenotypical differences. The removal of the phylogenetic effect shows that temperature and precipitation variables were strong determinants of the shell size and shape of the species, which explained the morphological differentiation (Fig. 2).