Herbivore corridors sustain genetic footprint in plant populations: a case for Spanish drove roads

Studied species

Plantago lagopus L. (Plantaginaceae) is a diploid polycarpic annual or biennial forb, abundant in grazed annual-dominated pastures of Central Spain (Peco et al., 2005). This species presents a persistent soil seed bank (Peco et al., 2003), gynodioecous breeding system and its seeds are dispersed by cattle endozoochory (Malo & Suárez, 1995). By sheep and under field conditions, epizoochoric dispersal is anecdotal, but endozoochory figures range from three to seven seeds per g dry weight of manure in late spring when dispersal is highest and flocks are moving northward, to 0.1 seeds in autumn when flocks are moving southward (P Manzano, 2004–2006, personal observations). This translates into ca. two million P. lagopus seeds moved northward per transhumance day, or 50,000 seeds moved southward, for a typical flock of 1,000 sheep. Experiments have also shown seeds are largely viable after sheep ingestion (Peco, Lopez-Merino & Alvir, 2006) and remain attached for a long time on sheep coats (De Pablos & Peco, 2007), even under real transhumance conditions (Manzano & Malo, 2006). This species is self-compatible, with both wind and insects as major pollination agents (Sharma, Koul & Koul, 1993). Annual-dominated pastures in Central Spain are known to present particularly high diaspore availability related to summer drought but translate into massive dispersal by livestock, particularly at the time when transhumant herds are on the move (Manzano, 2015a). Such large dispersal, along with the faster dynamics in annual plants (overlapping generations in perennial plants may disturb some genetic descriptors due to the mixture of genotypes), make them particularly suited to test genetic signals of dispersal in fragmented scenarios (Ewers & Didham, 2006).

Landscape scenario and population sampling

The study area is located between Quintanar de la Orden, Tomelloso and Villarrobledo municipalities in Castilla-La Mancha, Spain (coordinates: 39°18′03.17″N and 2°49′54.50″W, Fig. 1). The area is a plateau (830–900 m a.s.l.) of sandstones, loams and clay materials under continental Mediterranean climate, with a mean annual rainfall of about 500 mm and severe summer droughts. The vegetation is a mosaic of agricultural dry cereal and sunflower croplands, abandoned fields and dry grasslands used by local flocks. The area is crossed by the CDR, one of the major road droves (ca. 410 km long) that are still used for transhumant sheep and cattle herds that move every year from the cooler and wetter mountains of Teruel, Cuenca and Guadalajara provinces to the wintering dehesas in Sierra Morena at lower altitude (Oteros-Rozas et al., 2012). The herds cross the study area twice a year (northward in June and southward in November), currently comprising about 9,000 heads of sheep and 1,200 of cattle. The size of transhumant herds, dominated by merino sheep, has shrunk from about half a million in the 16th and 17th centuries to 100,000–200,000 in the next two centuries and around 20,000 during the second half of the 20th century (Bacaicoa Salaverri, Elías Pastor & Grande Ibarra, 1993). Productivity of the area is low except for peaks in May and October, because of cold winters and hot dry summers. Consequently, the numbers of transiting nomadic livestock have been historically much higher than those of resident livestock (Manzano Baena & Casas, 2010), which only grew in the last centuries (Bacaicoa Salaverri, Elías Pastor & Grande Ibarra, 1993). Such resident livestock herds are, in turn, kept at stables in urban settlements at night and graze grassland patches that are spread across the landscape (Fig. 1), overwhelmingly under private tenure. They graze mostly during wetter times in spring and autumn, the time of the year when local vegetation is growing (cf. Manzanares town in Manzano Baena & Casas, 2010) and diaspore availability is very low—in spring because plants have not yet produced seeds, and in autumn because rains incorporate the dead plant matter and any remaining diaspores into the soil. At other times of the year they are supplemented with hay, concentrate and water. Drove roads tend to avoid urban settlements and run rather tangentially (Fig. 1), so resident livestock do little use of them because of the rather radial pattern of land use, with urban settlements in the center. As a result, management of local resident livestock is likely to have some effect on pollen flow because of available pastures and to cause some short-distance dispersal across them, but it is not likely to have effects on long-distance seed dispersal along the drove roads.

In winter 2013 we selected 13 populations of P. lagopus (six in the CDR and seven in isolated pastures interspersed in the agricultural matrix). Populations located in the agricultural matrix were located at least five km from the CDR. Each population consisted of a group of plants not further than 100 m from each other, and all populations were at least five km apart from each other (Fig. 1; Table 1). For each population, we randomly collected between 8 and 30 flowering individuals of P. lagopus, whose leaves were stored under dry conditions until DNA extraction. Sampling was done without prior permission, in accordance with article 32d of the regional regulating law on drove roads regarding common complementary uses, including educational and formative activities in the environmental and cultural fields (Junta de Comunidades de Castilla-La Mancha, 2003). Authorization is required for such activities when they threaten the protection of sensitive ecosystem, forests with high wildfire risk, and protected flora and fauna species, and infringements are foreseen in its article 41—but our sampling in the ecosystem was negligible, and P. lagopus is neither an endangered nor a protected species. We deposited five individuals at Real Jardín Botanico de Madrid—CSIC herbarium (reference numbers MA-01-00892341–MA-01-00892345).

To test for the potential effect of other landscape factors on genetic descriptors, we measured the Euclidean distance from the sampling points to the nearest human settlement using Google Maps. At the agricultural matrix points, we also checked the time since last plowing by inspecting historic aerial photographs from the area (Instituto Geográfico Nacional, Spain)—the oldest of which being 70-years-old. The CDR, as in the case of all other major drove roads in Spain, has been protected from plowing for the last centuries (Bacaicoa Salaverri, Elías Pastor & Grande Ibarra, 1993).

Microsatellite development and laboratory protocols

We developed specific microsatellite markers for P. lagopus in collaboration with the Unidad de Genómica, Parque Científico de Madrid (http://fpcm.es/en/servicios-cientificos/). Briefly, the library was developed from DNA of five individuals of P. lagopus extracted using QIAGEN plant Extraction Kit, using manufacturer’s protocol. DNA was fragmented with ultrasounds using Bioruptor (Diagenode, Liège, Belgium) and size fragment was selected for an approximated size of 1.2 kb. DNA sequencing was performed using a 454-FLX-Plus Roche sequencer. We obtained approximately 185,000 sequences. We used Pal_Finder to select those sequences with more than five repetition motifs and a size that ranged from 100 to 300 bp. We used Primer3 (Untergasser et al., 2012) to configure potential primers that amplify those microsatellite regions. Details of the primers used are shown in Table S1.

Genetic descriptors: diversity, inbreeding, differentiation and structure

We checked the six microsatellites used for genetic errors and misidentifications of alleles with Micro-checker (Van Oosterhout et al., 2004). We estimated four genetic diversity indices: observed and expected heterocigosity (H_O and H_E) using INEST 2.1 (Chybicki & Burczyk, 2009), the Shannon information index (I), as a measure of genetic diversity that considers the number of alleles (NA) and its dominance in the population and the number of private alleles (A_P) using GeneAlEx 6.5 (Peakall & Smouse, 2012). Because the number of sampled individuals is different in each population (ranging from 8 to 30, Table 1), we also calculated the genetic diversity using Hp-Rare (Kalinowski, 2005) to rarefact the allelic richness and the number of private alleles to the minimum number of individuals sampled in one population (i.e., eight individuals) for better comparisons between populations according to population size, hereafter A_r and PrivRare, respectively. Linkage disequilibrium among loci was assessed with GenAlEx 6.5 (Peakall & Smouse, 2012).

We calculated the inbreeding coefficient, F_IS, using INEST 2.1. (Chybicki & Burczyk, 2009) that corrects the effect of null alleles, inbreeding and genotyping failures. DIC values between models with and without inbreeding were compared to evaluate the effect of inbreeding over F_IS. Departure from Hardy–Weinberg (HW) equilibrium with GeneAlEx 6.5 (Peakall & Smouse, 2012) to evaluate the inbreeding status of P. lagopus populations. We also assessed the presence of genetic bottlenecks generated by population dynamics (e.g., colonization events, population size reduction) using Bottleneck 1.2.02 (Piry, Luikart & Cornuet, 1999), through a Wilcoxon test with 2,000 permutations under the three plausible mutations models, as recommended for microsatellite data. We applied Bonferroni correction for multiple testing.

We estimated genetic differentiation through the F_ST coefficient calculated with FreeNA (Chapuis & Estoup, 2006) to minimize potential biases caused by null alleles or allele dropout. We tested signification between each pair of populations with 999 permutations. We used HIERFSTAT 0.04 (Goudet, 2005) to evaluate the F-statistics grouping of the populations according to their location (either within the matrix or in the drove road) with 1,000 permutations. To assess the presence of isolation-by-distance (IBD) patterns, we performed a Mantel test with 999 permutations using genetic differentiation (i.e., F_ST coefficients) and Euclidean geographical distance between populations with Vegan package in R (Oksanen et al., 2007). Mantel tests were also developed considering the location of the populations (Drove Road or Matrix) with the same permutations. Genetic differentiation between patches was compared with the calculation of migration rates or number of migrants. BayesAss (Wilson & Rannala, 2003) and Migrate (Beerli, 2009) were used to estimate recent and past migration rates, respectively, using the default parameters recommended by the manual. This might show differences among ancestral or recent gene flow, given the progressive decrease of livestock movement along the drove road during the last decades. Mantel tests were also developed considering the location of the populations (Drove Road or Matrix) with the same permutations.

We assessed genetic structure using the Bayesian clustering method STRUCTURE v. 2.3.4 (Hubisz et al., 2009) with prior information on populations. We carried out 10 independent runs for each K value (i.e., number of groups or clusters), ranging from 1 to 15. Each run comprises a burning period of 100,000 permutations, followed by 1,000,000 Markov Chain Monte Carlo iterations. We assumed correlated allele frequencies and genetic admixture. To assess the most plausible number of clusters (K), we used the STRUCTURESelector module (Li & Liu, 2018) to obtain plausible values of K. We considered the Evanno method (Evanno, Regnaut & Goudet, 2005) with the ad hoc statistic ΔK and identified the maximum ΔK when plotted against K. MedMeaK (median of means), MedMedK (median of medians), MaxMeaK (maximum of means) and MaxMedK (maximum of medians) following Puechmaille (2016) were also estimated (threshold value = 0.6). These estimators are useful to discard spurious clusters and detect substructure in case of uneven sampling (Puechmaille, 2016). We also used Clump v.1.1.2 (Jakobsson & Rosenberg, 2007) to estimate the membership coefficient of each individual in each cluster and Distruct v 1.1 (Rosenberg, 2004).

Relationships of genetic descriptors with environmental factors

To assess whether movement of transhumant herds along the drove road leaves a footprint, we selected the genetic descriptors that have most commonly been used previously in similar studies, i.e., F_IS and A_r (Rico, Boehmer & Wagner, 2014a; Rico et al., 2014b; Rico & Wagner, 2016; DiLeo et al., 2017; Plue, Aavik & Cousins, 2019). We built multiple linear models fitted with generalized least squares (gls function, nlme R-package) to examine the relationship between genetic diversity metrics with the population location (agricultural matrix or drove roads), the distance to nearby human settlements (in km) and the interaction location × distance.

Due to the effect of pollinators activity in the crossing system (i.e., selfing, partial selfing or outcrossing, Turner, Stephens & Anderson, 1982; Dutech et al., 2005; Van Etten et al., 2015), for the analysis of F_IS the amount of grasslands in 500-m radius circles centered on sampling points (in percentage, log-transformed to achieve normality) was also included. The landscape descriptors (distance and percentage of grasslands) were standardized and centered to interpret model coefficients as standardized effects. Spatial autocorrelation was taken into account by building the GLS models with three alternative spatial covariance structures most commonly used to describe spatial correlation (Gaussian, Exponential and Spherical). The presence of nearby grasslands was included for the analysis of F_IS, for this genetic descriptor is considered to be particularly influenced by pollinator activity due to its effect over selfing and/or kinship crosses (Rico et al., 2014b). Long-distance seed dispersal will rather happen along drove roads and through transhumant livestock that bypass such nearby grasslands, because the resident livestock that use them tend to have short-distance movements around human settlements and do not use the drove roads that run rather tangentially to the settlements.

The differences in genetic diversity among populations may show a nested pattern, suggesting that former connections between them have been lost (Aggemyr & Cousins, 2012; Plue, Aavik & Cousins, 2019). Nestedness appears when the set of alleles in communities with lower genetic diversity are a subset of those found in communities with higher genetic diversity, which would imply a deterministic loss of alleles. Thus, we explored whether such nested patterns occurred among our study populations and if they could be linked to the presence of the drove road, the distance to settlements and the amount of neighboring grasslands. If nestedness among populations is caused by loss of connectivity, we expect sites within the agricultural matrix and far from settlements to have lower allelic richness, lacking some particular alleles present in the populations of the drove road (in other words, sites within the agricultural matrix would have a higher degree of nestedness). Neighboring grasslands could act as a local source of propagules (seed and pollen) deterring the loss of alleles, and thus the larger the availability of neighboring grasslands the less likely we would find nestedness patterns. To describe nestedness patterns we used the NODF index (Almeida-Neto et al., 2008). We tested if nestedness with a suitable null model based on 499 randomized matrices (null model “r2” with package vegan in R which constrains row totals –allelic richness—to equal the observed distribution and samples columns –alleles—with probability proportional to squared incidence frequencies). Then, we modeled the population rank order, from the highest rank of the population with most alleles to the lowest rank of the population with least alleles, with linear models built with GLS as above. The populations differed in the number of individuals analyzed. Then, to balance sampling effort among them we repeated 50 times the above procedure, by randomly selecting the same number of individuals per population. Final NODF, p-values and rank estimates were the mean of the above procedure for these 50 balanced replicates.

For the three genetic diversity metrics (F_IS and A_r) and for nestedness, we did not try to simplify the models but to identify the relevant predictors with p-values based on nested log-likelihood ratio tests. The genetic descriptors we used as response variables in the models were the inbreeding coefficient (F_IS), the genetic diversity statistic (A_r) and the nestedness estimates. We confirmed nearby grasslands to be only a relevant factor for F_IS by comparing models of the three metrics with and without the grassland factor via AICc. Additionally, we tested the correlation of the genetic descriptors with the plowing history (years since last plowing) with base R and Anova function, car package and just for the populations in the agricultural matrix—the drove road having been protected from plowing for the last centuries.

Genetic diversity and inbreeding

The values for genetic diversity (I, A_P, H_O, H_E, A_r and PrivRare) were similar between populations in the agricultural matrix and populations in the drove road, except for A_P, which was greater in the populations of the agricultural matrix, due to the large NA found in LMA08 (Table 1). The average (CV) expected heterozygosity H_E was 0.88 (0.02) in drove roads populations and 0.89 (0.03) in matrix populations. However, matrix populations showed greater variability in H_E, with values ranging from 0.92 in LMA13 to 0.84 in LMA08. The same pattern (similar average but different variability) was found in the other genetic diversity indexes (Table 1).

The inbreeding coefficient (F_IS) varied greatly among populations and between matrix and drove road locations. Drove road populations showed a mean value of 0.08 (ranging from 0.01 in LCR04 to 0.17 in LCR01), while matrix populations showed a mean value of 0.11 ranging from 0.06 in LMA05 to 0.17 in LMA14 (Table 1). Inbreeding had an important effect over Fis in LCR01, LMA11, LMA14 and LMA16 populations according to INEST model estimation. All populations had, at least, one locus that deviated from HW equilibrium (four loci in LCR01). Bottleneck signals were found in some populations under SMM model (LCR09, LMA08 and LMA11), but the results were not significant in all genetic models and when applying multi-comparisons corrections (Table 2). Linkage disequilibrium among loci showed only four significant comparisons, which were non consistent between loci pairs.

Genetic differentiation and structure

Pairwise genetic differentiation values (F_ST) between populations were low (Table S2). On average, LMA14 showed the greatest differentiation, with an average F_ST value of 0.08. Hierarchical analysis of the genetic differentiation showed low differentiation between matrix and drove road locations (F_CT = 0.032) and a moderate one between populations within each category (F_SC = 0.11). Both values are significant. Mantel test (r = 0.13; p = 0.22) did not detect the presence of an IBD pattern, but results were significant when only analyzing the six populations included in the drove road (r = 1; p < 0.01). The genetic structure analyzed with a Bayesian clustering approach (i.e., STRUCTURE, Fig. 2) showed different values depending of the approach consider to estimated K. The estimation of ΔK proposed values of K = 12 and K = 4 as more plausible approaches. MedMedK, MedMeanK and MaxMaxK also pinpoint that K = 4 is the most appropriate distribution of the genetic structure and MaxMedK suggested that K = 5 is the most accurate result. The K = 4 distribution showed an independent cluster (green) for LMA14 population, another cluster (purple) for LMA16, LMA13 and LCR10 populations and a third cluster (orange). The other populations are assigned to the last cluster, with admixture of the others. The K = 12 distribution showed an important population admixture and one cluster assigned to LMA14 population (Fig. 2). The results for gene flow (BayesAss and Migrate, Table 3; Fig. 3) showed differences for the migration rates estimated, with higher values for past gene flow (i.e., Migrate) than contemporary (i.e., BayesAss). The evaluation of the recent gene flow (i.e., BayesAss), focus the main differences in the population LMA08 that received a significant number of migrants and in lesser extent, LMA16 and LCR01. The remaining populations showed similar migrations rates between populations. The ancient gene flow analyzed (i.e., Migrate) showed similar mean values for drove roads and matrix populations but important differences between populations. Although LMA08 maintains high values of gene flow, LCR02 received the highest average migration rate meanwhile LCR10 showed the smallest average value.

Landscape factors governing genetic descriptors in P. lagopus populations

The distance to the nearest settlement had an impact on both genetic descriptors analyzed (Table 4). This relationship is independent of the frequency of plowing, because for the agricultural matrix plots the correlations of years since last plowing and the residuals of the regressions explaining the genetic descriptors were non-significant. Also, the values of genetic descriptors did not correlate significantly with the number of years since last plowing (Pearson correlations were 0.07, p = 0.89 for F_IS; and −0.13, p = 0.78 for A_r).

In addition to the distance to the nearest settlement, the interaction with population position showed a significant effect on inbreeding (F_IS; Table 4). Thus, F_IS increased with distance from human settlements in the agricultural matrix, but followed the opposite trend in drove road populations (Fig. 4A). The presence of nearby grasslands also proved to be a significant factor, showing higher rates of inbreeding at higher grassland cover (Table 4).

For the descriptor of genetic diversity (A_r), the distance to settlements had a consistent negative and significant effect, while the effect for population position and its interaction with distance was only marginally significant (Table 4). A_r values between matrix and drove road positions showed a pattern of convergence at short distances and divergence at high distances from human settlements, with the agricultural matrix then showing less diversity (Fig. 4B).

Significant nestedness was found among the P. lagopus populations (mean observed NODF = 34.7, mean p-value = 0.02, all averaged from 50 random selections of individuals to balance sample size among populations). None of the geographic parameters, however, significantly affected the rank (Table 4).

Implications for ecosystem function and conservation

Functional herbivore corridors would sustain diversity by reducing the risk of local extinctions due to increased inbreeding and loss of genetic variation in small populations (Caughley, 1994). The effect of drove roads on inbreeding observed in our study has relevant implications for conservation, as inbreeding depression is a major concern for the conservation of biodiversity in human-altered fragmented landscapes (Hedrick & Garcia-Dorado, 2016). We think that connectivity factors, related to herbivore type and capability of sustaining dispersal corridors for plants, as well as optimal matching of grazing times with the flowering and seeding phenology of plants, would therefore add to the complexity of factors that influence the relationship between grazing and plant diversity. Such factors include productivity and disturbance (Herrero-Jáuregui & Oesterheld, 2018; Olff & Ritchie, 1998), historical effects on communities (Cingolani, Noy-Meir & Díaz, 2005; Osem, Perevolotsky & Kigel, 2002) or herbivore species and sizes (Bakker et al., 2006; Liu et al., 2015). Functional herbivore corridors would sustain diversity by reducing the risk of local extinctions due to loss of genetic variation in small populations (Caughley, 1994).

The marginally significant, yet observable trend of drove roads in reducing genetic diversity loss of P. lagopus, in spite of the rapid transit of transhumance flocks—ca. 15 km/day or 1.5 km/h when moving (Hevia et al., 2013) is very noteworthy. Such a use, which is intense but nevertheless very restricted in time, would achieve extraordinary dispersal effects due to a “green wave effect” (Merkle et al., 2016), as transhumant flocks use drove road pastures when they are at their maximum productivity (Manzano Baena & Casas, 2010) and therefore at their maximum seed production stage—which would add to the reduced impact on pollinators we observed that is also related to the “green wave effect.” This is a radical contrast with other works studying rotational grazing systems (Rico, Boehmer & Wagner, 2014a; Rico et al., 2014b; Rico & Wagner, 2016; DiLeo et al., 2017; Plue, Aavik & Cousins, 2019) where the intensity of use per unit of area is much higher and the chances for consolidating dispersal between patches are also higher, and the effects on pollinator community are less clear. Indeed, a confirmation of such trend in future studies that are able to study more populations than we did would highlight the potential role of herbivore corridors to promote long-distance dispersal. The higher estimated migration rates in the past than in the present are also an alert on the negative effect that declining transhumant practices may have on ecosystems.

Mobile pastoralism is globally threatened by non-conductive policies (Manzano, 2017), weakening its role as a provider of important ecosystem functions such as habitat provision, nutrient cycling or bush encroachment/fire control (Leroy et al., 2018). The potential magnitude of the threat is great, for such land use covers more than half of the global land (Manzano, 2015b). A derived disruption of dispersal processes and the structural role for pollinators played by corridors themselves may further contribute to a generalized grassland biodiversity crisis motivated by land use intensification (Gossner et al., 2016) and climate change—and on the latter, a potential decoupling of herbivore “green wave” movements from the phenology and seed availability brings also concerns (Berg et al., 2010; Thackeray et al., 2016). It is urgent to quantify such dispersal beyond our study in order to orient adequate science-evidenced policy recommendations (Manzano Baena, 2012). Furthermore, corridors per se are important beyond livestock. We obtained our data in a drove road that has experienced a drastic decline in use in the recent centuries, from 600,000 head in the 16th century to less than 10,000 nowadays (Bacaicoa Salaverri, Elías Pastor & Grande Ibarra, 1993; Hevia et al., 2013). The effect on the genetics of plants must have been much more significant at that time, similar to what they must be in the systems presently steered by hundreds of millions of mobile pastoralists (Manzano & Agarwal, 2015), whose capacity to cover wide territories is shown by the fractal structure of their drove roads in different countries (Manzano-Baena & Salguero-Herrera, 2018). This effect should also be greater in the world’s greatest wild herbivore migrations (Berger, 2004).