Limited effect of a highway barrier on the genetic structure of a gypsum soil specialist

Study species

Lepidium subulatum L. (Brassicaceae) in an endemic chamaephyte of the Iberian Peninsula with a few populations in North Africa, specialist of gypsum soils (Hernández Bermejo & Clemente, 1993). It is a perennial shrub (up to 26 years; Eugenio et al., 2012), 20–60 cm in height (Palacio & Montserrat-Martí, 2005; Romão & Escudero, 2005). It is a generally entomophilous and primarily self-incompatible species (Matesanz et al., 2015). The species is diploid with eight chromosome pairs (2n = 16; Hernández Bermejo & Clemente, 1993). Outcrossing and atelechorous seed dispersal with seeds covered by a mucilage able to anchor the seeds in the vicinity of mother plants surely affects its ability to connect habitat remnants (Gómez-Fernández, Alcocer & Matesanz, 2016).

Study area

The area is a homogeneous landscape at the south-east corner of Madrid, central Spain. It is characterised by the presence of fragments of gypsum vegetation, surrounded by an intensively managed matrix of dry, extensive crops (coordinates in UTM ETRS89 30N: the centroid of the most north-westerly population (Pop 13) 4,438,624, 489,512; the most south-easterly population (Pop 12) 4,432,335, 495,277). The climate is arid, with an annual average temperature and rainfall of 13.8 °C and 440 mm, respectively and almost no rainfalls during the summer (Meteorological Station of Arganda del Rey, 459,310, 4,462,676; 38–48 km away from the study area). The area is dissected northwest to south east for more than 15 kms by the A-3 highway, which has operated for more 35 years in its present configuration (Madrid-Valencia, a motorway with four traffic lanes and a width of 25 m; Boletín Oficial del Estado, 2016), dividing the area into two similar halves.

We chose 24 populations or patches distributed at both sides of the A-3 motorway, selecting the same number of populations on each side of the area in a strio of 3 km from the highway (Fig. 1) and maximising the range of fragment sizes and connectivity. Populations were selected from the kilometre point 65 to 75 along the highway. In each gypsum patch we collected samples from 8–20 individuals (average 16.2 ± 4.2 individuals) for DNA extraction (Table S1). We calculated L. subulatum population density as the number of intersected individuals in a transect of 250 m across the patch. Remnants with insufficient number of individuals were discarded to avoid confounding effects on flow arising from small population sizes. For each remnant patch, the following variables were also estimated fragment area, patch perimeter and distance to the highway from the patch’s edge and from the centroid using ArcGis 10.8 (ESRI, 2020). In addition, we calculated the connectivity among patches with Tremlová and Münzbergová’s connectivity index (Tremlová & Munzbergová, 2007). This was done taking into consideration all the remnants present in the area (see Arenas et al., 2017 for details).

In order to connect genetic parameters with the physical-chemical characterisation of each patch (remnant’s quality), we collected two soils samplings per patch with a core ten centimetres deep and six centimetres in diameter. Samples were analysed for total amount of nitrogen (N), phosphorous (P), potassium (K), organic carbon, pH, conductivity and glucosidase and phosphatase’s activity in Nutrilab/URJC (https://nutrilab-urjc.es/). We subjected the samples to a sulphuric acid digest and then used an SKALAR San⁺⁺ Analyser to estimate N and P (Skalar, Breda, Nertherland; Maestre & Puche, 2009). For K, samples were shaken with distilled water in the ratio 1:5 prior to digestion (Maestre & Puche, 2009). We used the Walkley–Black method (Gelman, Binstock & Halicz, 2012; Sato et al., 2014) for estimating organic C. We measured pH and conductivity using a solution 1:2.5 (weight:volume). We used the protocols described in Maestre & Puche (2009) to measure glucosidase and phosphatase activity.

DNA extraction, amplification and sequencing

Material collection was approved by the Consejería de Educación e Investigación de la Comunidad de Madrid (REMEDINAL-TE CM, P2018/EMT-4338). The genetic study was performed with eight specific microsatellites markers: Lsub01, Lsub02, Lsub03, Lsub04, Lsub05, Lsub07, Lsub08, Lsub12 (Martínez-Nieto et al., 2012). DNA was extracted from 60 mg of dried young leaf tissues using the SpeedTools Plant Extraction kit (Biotools, Madrid, Spain). We amplified individually each microsatellite using Polymerase Chain Reactions, mixing 14.2 μl of milli-Q water, 2 μl of T10X buffer with MgCl₂, 0.8 μl of dNTP mix [10 mM], 0.5 pmol of each primer (forward labelled with fluorescence dye and reverse) 1 U of Taq-polymerase (Biotools, Madrid, Spain) and 1.2 μl of DNA (DNA had a concentration range of 10–30 ng/ul). The thermocycler’s protocol was: 4 min at 95 °C, followed by 30 cycles of 45 s at 95 °C, 45 s at 51 °C (annealing temperature) and 1 min at 72 °C, followed by a final step of 7 min at 72 °C. PCRs were checked using electrophoresis in agarose gel pre-stained with Gel Red (Biotium, Hayward, CA, USA). Lastly, the amplified fragments were analysed with an ABI 3730XL sequencer (Applied Biosystems, Foster City, CA, USA) of the Unidad de Genómica (Universidad Complutense, Madrid, Spain).

Statistical analysis

We tested the Hardy-Weinberg equilibrium with Genepop 4.7 (Raymond & Rousset, 1995; Rousset, 2008) and the presence of null alleles by MicroChecker 2.2.3 (Van Oosterhour et al., 2004). We also calculated the allelic richness per population/patch, the total number of private alleles, the expected and observed heterozygosity (H_E and H_O), the inbreeding coefficient (F), and the number of migratory individuals (Nm), using GenAlEx 6.5 software (Peakall & Smouse, 2006). Complementarily, migration rates or number of migrants were estimated with BayesAss 3.04 (Wilson & Rannala, 2003) and Migrate 4.4.3 (Beerli, 2009), to estimate recent or past gene flow rates respectively. We calculated the F_ST with a null alleles correction with FreeNA (Chapuis & Estoup, 2007; Chapuis et al., 2008). Due to the different numbers of individuals sampled in each population, the allelic richness was calculated by rarefaction analysis using HP-Rare v 6 June 2006 programs (Kalinowski, 2005). To evaluate the existence of any genetic bottleneck, we used the Wilcoxon test in Bottleneck v.1.2.02 software (Cornuet & Luikart, 1996a). This measures deviations in heterozygosity and deficiencies in allelic richness. It takes into consideration that, when a bottleneck is produced, the effective population size drops, triggering a decrease in the allelic number that exceeds the effects of heterozygosity (Cornuet & Luikart, 1996b). We run three different mutation models: infinite allele model (I.A.M.), step wise-mutation model (S.M.M.) and intermediate two-phase mutation model (T.P.M.). We used the parameters by default for T.P.M. (70% of proportion of S.M.M. and 30% of variance). For each model, we performed 2,000 iterations. We also evaluated whether any of the genetic parameters (i.e., allelic richness, private alleles, H_O, H_E, F, F_ST average or Nm average) were affected by some soil variables (i.e., N, P, K organic carbon, pH, conductivity, phosphatase and glucosidase) which are surrogates of habitat quality, connectivity (C_i), other landscape variables (fragment size, perimeter of the remnant and minimum distance to highway) and population density. For this, we utilised General Linear Models (GLMs) using R v.3.3.0 (R Core Team, 2016); in which we tested a model with all independent variables but without interactions among predictors. P-values were corrected using Bonferroni’s adjust to avoid the multiple comparison effect.

We studied the genetic structure of the 24 populations using a Bayesian clustering method with STRUCTURE v.2.3.4 with prior information on population membership. Ten independent runs were carried out for each K value (i.e., number of potential clusters, ranging from 1 to 27), each one with a burn-in period of 10⁵ iterations followed by 10⁶ MCMC iterations and 10⁴ thinning. Analyses were developed with allelic correlation frequencies and genetic admixture. According to the different population sampled, we used the STRUCTURESelector (Li & Liu, 2018) module to obtain K values. We considered the Evanno method (Evanno, Regnaut & Goudet, 2005), together with MedMeaK (median of means), MedMedK (median of medians), MaxMeaK (maximum of means) and MaxMedK (Maximum of medians) proposed by Puechmaille (2016) with a threshold value = 0.6, to discard spurious clusters and detect potential substructure. Also, to complement the STRUCTURE analysis (Puechmaille, 2016; Janes et al., 2017), we performed a spatial Principal Components Analysys by adegenet v. 2.2.1 package in R (Jombart, 2008; Jombart & Ahmed, 2011) to evaluate the global and local spatial structure.

To evaluate the effect of the highway as a modifier agent of the genetic landscape, genetic diversity and differentiation was estimated, considering the two sides of the motorway (NE vs. SW population groups) separately. In addition, an AMOVA (Analysis of MOlecular VAriance), performed by GenAlEx 6.5 software (Peakall & Smouse, 2006), was used to estimate and evaluate the distribution of the percentage of intra- and interpopulation variability, including both highway sides as different regions. These results also complemented the genetic structure analyses (see below). To test if the highway triggered a barrier effect between both sides, we calculated the average of F_ST of one side (northeast or southwest) with each population of the other highway side to observe if the value was significantly different with the F_ST average of the first side. We also analysed the effect of highway and geographic distance between pairs of populations using Generalised Least Squares (GLS) models with the nlme and corMLPE packages (Pinheiro et al., 2018; Pope, 2019). The corMLPE package constructs regressions using distance matrices with non-linearity and random effects. In this case, we selected two populations and used as independent variables the side of highway (binary variable: 0 = same side of the motorway; 1 = different side) and the distance between populations. We selected the best model using AIC, using F_ST as the dependent variable and geographic distance, highway side and their interactions as the independent variable.

Patch genetic characterisation

None of populations in remnants were in H-W equilibrium, except population 13 and 14 (Table S1). Some microsatellites could present some null alleles (Table S1). We found an average allelic richness per patch/population of 7.3 ± 1.3 alleles. Population 21 presented the highest value (n = 10.4 ± 1.8) and population 4, the lowest (n = 4.6 ± 1.2) (Table 1; Table S1). When the average number of alleles was rarefied, the average richness decreased to 4.9 ± 1.1 (Table 1). In this case, population 4 continued as the poorest (n = 3.8 ± 0.9), whereas population 12 was the richest (n = 8.3 ± 0.6) (Table 1; Table S1). The total number of private alleles was 35, but several populations did not have any (populations 1, 8, 12, 15, 16, 17, 18, 19, 20) (Table 1; Table S1). In contrast to allelic richness, population 4 had the highest number of private alleles (n = 5) (Table S1). Expected heterozygosity was higher than observed in all populations (average H_E = 0.7588 ± 0.0 and average H_O = 0.6507 ± 0.1, respectively). Population 21 showed the highest expected heterozygosity, while population 4 had the lowest (H_E = 0.8±0.0 and H_E = 0.7±0.1, respectively) (Table S1). Nevertheless, the population with the highest observed heterozygosity was population 5; and the lowest, population 12 (H_O = 0.8 ± 0.2 and H_O = 0.5 ± 0.3, respectively) (Table S1). Finally, the average inbreeding coefficient (F) was 0.1 ± 0.1, with great variation between populations. Population 5 had the lowest inbreeding coefficient, showing a patent heterozygosity excess (F = 0.0 ± 0.2); meanwhile, population 12 had the highest coefficient (F = 0.3 ± 0.3) (Table S1). The presence of bottlenecks was detected for all populations. (Table S3).

Genetic differentiation and migrants: highway effect

The average number of migratory individuals between population/patches (Nm) was 4.02 (Table 1). Population 22 had the highest number of migratory individuals (Nm = 5.6); and population 12, the lowest (Nm = 1.4) (Table S1). Genetic differentiation based on F_ST coefficients, showed an average value between population pairs of 0.1, suggesting moderate genetic differentiation across the study area (Table S1) suggesting certain degree of isolation.

When comparing the populations located on either side of the highway (see Fig. 1), there were significant differences in the number of private alleles and number of migratory individuals. In the case of private alleles, the populations of the SW side presented a higher number (n = 22 and n = 13, respectively) (Table 1), while the NE half had a higher number of migratory individuals (Nm = 5.1and Nm = 2.8, respectively). The NE population with the highest number of migratory individuals was population 22 and the population with the lowest was population 24 (Nm = 6.6 and Nm = 3.2, respectively). In contrast, in the SW side, the population with the highest number of migratory individuals was population 1 whereas the lowest number was found in population 12 (Nm = 3.8 and Nm = 1.2, respectively). A similar pattern was found in the evaluation of recent gene flow (i.e., BayesAss). All migration rates were smaller than one migrant per generation, although higher rates were found in populations at both sides of the highway (i.e., population 15 at NE side and population 7 at SW side), forming small clusters (Table S4). Ancient gene flow showed a similar pattern (i.e., Migrate), with populations with high and small migration values at both sides of the highway (Table S5). Only two populations (pop 15 and 3) showed high values of gene flow with both approaches.

The best GLS model or explaining differences in F_ST values contained as predictors the Euclidean distance between patches or the highway side (Table 2). Therefore, F_ST was higher when geographic distance among populations was higher or when populations belonged to different highway side. However, these differences should be considered with caution, due to the small reduction in AIC values when this last variable was included in the model (Table 2).

Genetic structure and association with environmental variables

Bayesian analyses suggested different number of genetic clusters (K). Namely, the Evanno method suggest K = 4 as the most plausible value. MedMedK and MedMeaK support a K = 7 and MaxMedK and MaxMeaK indcated that K = 8 (Fig. 2). These differences suggest the existence of some internal substructure within the main 4 clusters. This heterogeneity in the number of clusters is the result of internal divisions suggesting the existence of nested structures. For instance, some of the clusters proposed by K = 4 split into different clusters with K = 7 or K = 8. The cluster that groups populations 16–24 are divided into three different clusters in K = 7 or 8 (populations 17,18 and 19 in one cluster, population 24 in by itself, and the other populations in another cluster). In contrast, sPCA showed a global structure (p-value = 0.01) but no a fine scale pattern (p-value > 0.05; Figs. 3C and 3D). This was surveyed in the map of genetic clines (Fig. 3A), where individual scores were similar in all the study area, except in populations 3, 15–19 and 11–12, and in the eigenvalues (Fig. 3B), where only it is remarked one among the positive values (global structure). In the case of the linear models, we found that some genetic descriptors, such as allelic richness (estimate = 0.001; p-value < 0.01) and expected heterozygosis (estimate = 0.01; p-value = 0.04), were positively related to the population density of the patch (Fig. 4). We did not find any other relation with the rest of genetic parameters and environmental variables.

Path genetic differentiation

Populations of L. subulatum presented significant variation in genetic diversity, gene flow and structure (e.g., allelic richness and heterozygosity). This variation might be associated with selective drivers that constrain the populations, but also to neutral processes that may foster different levels of genetic diversity in the remnant populations, including biological, demographic and ecological features. Although selective and neutral processes are acting, L. subulatum should sustain levels of genetic diversity that are enough to maintain plastic responses to different environmental conditions (Matesanz et al., 2020a, 2020b). L subulatum is a generalist entomophilous and primarily self-incompatible species (Loveless & Hamrick, 1984; Gómez-Fernández, Alcocer & Matesanz, 2016), with an important and long-standing presence in the soil seed bank (Caballero et al., 2005) and significant plant recruitment under field conditions (Escudero et al., 2000). These factors may support the relatively high levels of genetic diversity, even with the presence of selective forces and demographic or stochastic events which could tend to minimise this variation.

We detected evidences suggesting recent bottlenecks or founder events. Perhaps, due to these, some patches presented a particular concern due to the existence of local moderate F inbreeding coefficients. This probably is due to crossing with close relatives within patches (dispersal is very inefficient, see Escudero et al., 2000), causing an increase of kinship and, therefore, inbreeding (Fischer, Hock & Paschke, 2003). It is well known that inbreeding (or high kinship similarity) has a negative effect on population dynamics because it reduces the reproductive success of the species (Jolivet, Rogge & Degen, 2013; Hermansen et al., 2015) and the frequency of heterozygosity, fixing more easily any deleterious alleles (Mattey & Smiseth, 2015) which potentially decrease the adaptive ability of the species (Porcher & Lande, 2016). However, as the loci number is low, these F inbreeding coefficient values could be due to presence of null alleles.

Genetic differentiation and migrants: highway barrier

Our most remarkable result is the weak effects of the highway on the genetic structure (allelic richness, heterozygosity, genetic differentiation, gene flow parameters) of this species. This type of infrastructure could act as a barrier but could also improve the gene flow between the both sides due to the movement of vehicles or animals, with principal or auxiliary roads acting as dispersal pathways (Schmidt, 1989; Tikka, Högmander & Koski, 2001). Other potential causes that might explain the reduced effect of the highway could be related to the relatively small scale of the study (along 10 km of the motorway), the molecular tool employed (i.e., neutral markers, Balkenhol, Waits & Dezanni, 2009) and the fact that it would be necessary additional generations to verify any genomic signal in the populations (see Eugenio et al., 2012). It is worthwhile to note that some populations on the SW side of the study site might be suffering a potential barrier effect, as showed by the relatively high F_ST values between populations pairs (Table 1; Table S1). Therefore, it would be interesting to perform a detailed study to separate the possible existence of genetic differentiation derived from selection or demographic dynamics that could suggest potential local extinction risks.

Remnant populations retained most of the genetic diversity of the species in the area, which has been demonstrated as essential to maintain the capabilities of the species to cope with the uniqueness of this habitat (Matesanz et al., 2020a). The presence of certain levels of gene flow (present and ancient) and the moderate F_ST values could indicate that there is not a deficit in heterozygotes in the remnants. The exchange of genetic material could be mainly favoured by some pollination features (Aguilar et al., 2008; Matesanz et al., 2015; Gómez-Fernández, Alcocer & Matesanz, 2016) and by the cattle herds and their movements. The existence of dense soil seed banks (Eugenio et al., 2012) can also ameliorate the genetic differentiation among populations. On the other hand, in some patchy populations, there is a high number of private alleles. This genetic parameter is related to the amount of genetic flow among patches (Slatkin & Barton, 1989), and could indicate a reduction of genetic flow among populations, or at least, suggests some particular population dynamics in those patches. The number of private alleles is higher than those obtained in the other study of the species (Gómez-Fernández, Alcocer & Matesanz, 2016). This discrepancy could be due to differences in the study areas that favour or restrict genetic flow, for instance, the shape or the size of the gypsum habitat remnants.

Genetic structure and association with environmental variables

Estimates of the genetic structure differed slightly depending on the followed approach. We obtained K = 4 when using Bayesian analyses (i.e., STRUCTURE), while 20 neighbours were necessary to connect the network in the sPCA, although the global structure was also the most important. Nonetheless, both cases showed an important admixture of population of both sides of the highway. Consequently, the highway did not exert a significant barrier effect. This high level of admixture in the genetic clustering might be associated with the presence of different factors that foster the dispersal of the species and maintain the high levels of genetic diversity. Efficient seed dispersal among patches could be associated with the active presence of cattle or mainly stationary (or partially stationary) sheep herds within the area. Highways could work as a barrier but also improve plant dispersal, since vertebrates used roads to move, creating corridors (Suárez-Esteban, Delibes & Frediani, 2013). In this sense, several studies have showed the role that herbivores plays in seed dispersal in gypsum specialists (Pueyo et al., 2008). Sheep seem to have a remarkable effect in the short and long-distance seed dispersal (Manzano & Malo, 2006; García-Fernández et al., 2019), even when seeds do not have the appropriate surface structures (e.g., seed hooks; Bakker et al., 1996; Kuiters & Huiskes, 2010). Foraging dynamics could profoundly affect the structure of the whole landscape since sheep can graze remnants in an unpredictable pattern and occurs variably as a function of the season and palatable biomass.

Surprisingly, we found only marginally significant relationships between genetic parameters and the environmental or landscape variables. As it occurs with other edaphic specialists, L. subulatum only grows in gypsum soils and such specialisation probably acts as an almost binomial filter for the species; the species is excluded of non-gypsum soils but is very abundant in any adequate gypsum island (Palacio et al., 2007; Escudero et al., 2015). Nevertheless, only population density had a positive effect in the allelic richness and the expected heterozygosis, suggesting that patches with a higher number of individuals, independently of other soil or landscape conditions, show higher levels of genetic diversity. This, in accordance with controlled studies in common garden experiments, induces high levels of phenotypic plasticity, a necessity for facing stressful conditions (Matesanz et al., 2020a).