Non-linear genetic diversity and notable population differentiation caused by low gene flow of bermudagrass [Cynodon dactylon (L.) Pers.] along longitude gradients

Plant materials and DNA extraction

Plant materials consisting of 249 C. dactylon individuals were sampled at 13 different geographic origins locating 105°57′34″–119°27′06″E, China (Fig. 1). Mean annual precipitation and temperature at the collection sites were provided by the China Meteorological Administration (Table 1). Total genomic DNA was extracted using a Zymo Research Plant DNA Extraction Kit (Xinmaijie, Zhengzhou, China) and the quality and concentration of DNA were checked by electrophoresis on 1% agarose gels.

EST-SSR amplification

The 20 EST-SSR primer pairs from previous research (Zhang et al., 2019) produced clear and polymorphic bands and were subsequently used to evaluate genetic diversity for this study. 25 µL PCR amplifications reaction solution contains 1 µL of template DNA, 2.5 µL 10 × PCR buffer, 2 µL of 25 mM MgCl₂, 0.5 µL of 10 mM deoxy nucleoside triphos phate (dNTP), 0.5 µL each of forward and reverse primers, 0.2 µL of Taq DNA polymerase, and 17.8 µL of nuclease-free water.

The amplifications were performed with the following PCR protocol used in the previous research (Zhang et al., 2019): first 1 cycle of 3 min at 95 °C, followed by 10 cycles at 95 °C, 60 °C and 72 °C respectively for 30 s, finally terminating with 1 cycle of 6 min at 72 °C. Next, the samples were subjected to 1 cycle at 95 °C for 3 min, then 20 cycles of 30 s at 95 °C, 55 °C, and 72 °C respectively, and finally 1 cycle at 72 °C for 6 min. PCR amplified fragments labelled with one fluorescent tag (6-FAM (6-carboxyfluorescein) (blue) and HEX (green)) were detected on a 3730XL Genetic Analyzer (ABI, USA) with the size standard GS500 LIZ. Capillary electrophoresis analysis was made by using GeneMapper Software Version 4.0 (Applied Biosystems, Foster City, CA, USA).

Analysis of genetic diversity

Genetic diversity for all populations and ploidy levels was calculated using POPGENE V. 1.32 Software (Yeh et al., 1997). These indicators included Nei’s gene diversity index (He), Shannon diversity index (I) and polymorphic information content (PIC) (Smith et al., 1997; Yeh, Young & Boyle, 1999). He was generated according to genetic distance to analysis genetic diversity. I was computed based on the obtained allele frequencies, and PIC reflected the degree of variation of microsatellite DNA. Ploidy level of each sampled individual was determined by using a flow cytometer (Cube8, Partec, Germany) based on the method outlined by Zhang et al. (2020). Linear association between ploidy levels and longitude was investigated under a regression analysis using the software package SPSS13.0 for Windows (SPSS Inc. Chicago, IL, USA) considering correlation coefficient values and significance of variables in the model. All individuals of the same ploidy level were treated as a population to compare the genetic diversity among different ploidy levels. Significant differences in genetic diversity parameters among different populations and ploidy levels were evaluated under a one-way analysis by Duncan’s test using SAS 9.3 software (Institute Inc., Cary, NC, USA) (Duncan, 1955). We divided all the individuals into four groups by longitude from east to west to calculate the genetic diversity parameters for the four groups with four populations within one group as described above.

Analysis of genetic relationship among all populations of C. dactylon

Distance matrices of intra-population genetic diversity indicators, environmental distance matrices, genetic distance matrices and geographical distance matrices based on the dataset for 13 populations were calculated using GenAlEx 6.1 (Peakall & Smouse, 2006). Mantel correlation tests (Mantel, 1967) between intra-population genetic diversity indexes and environmental factors (temperature and precipitation) were conducted to examine the effects of environments on the intra-population genetic diversity. Mantel correlation coefficient between genetic distances and longitude distances for C. dactylon at different spatial scales was also investigated using NTSYS pc version 2.0 (Khan et al., 2002) and the probabilities for the significance were assessed using Bonferroni inequality adjustment. Partial Mantel test was conducted in R (ver. 4.0.5; R Core Development Team, 2011) to estimate correlation between genetic distances and geographic distance, while controlling for mean annual precipitation Analyses of molecular variance (AMOVA) and population differentiation index (FST) were estimated using ARLEQUIN v3.5 (Excoffier & Lischer, 2010). The significance of FST was determined by random resampling of the genotypic data through 1000 permutations. Principal coordinates (PCoA) analysis of genetic variations at different longitudes was performed using GenAlEx Version 6.1. With the help of clustering in NTSYSpc-v.2.1 software, UPGMA (unweighted pair-group method with arithmetic means) was used for cluster analysis to generate a dendrogram for genetic variation along longitude gradients (Rohlf, 2000).

Population structure analysis

The population structure of the 249 individuals was constructed using a Bayesian clustering method implemented in STRUCTURE 2.3 (Hubisz et al., 2009) as previously described in Zhang et al. (2019). The model-based approach implemented was used to subdivide the individuals into different subgroups. Specifically, twenty independent runs were performed for each simulated value of K ranging from 2 to 13. For each run, 10,000 iterations were set as the burn-in time following 100,000 iterations of Markov Chain Monte Carlo (MCMC). Outputs were processed with CLUMPP V1.1.2 (Jakobsson & Rosenberg, 2007) and STRUCTURE barplots were displayed using DISTRUCT 1.1 (Rosenberg, 2004). The best K value was determined from the modal value of the run with the delta K. Gene flow was calculated using POPGENE V. 1.32 Software (Yeh et al., 1997) based on the formula of Nm = 0.5(1- Gst)/Gst.

Mismatch distribution analysis (MDA; Rogers & Harpending, 1992) was used to investigate whether C. dactylon had undergone recent population expansion. Goodness-of-fit was tested with the sum of squared deviations (SSD) between observed and expected mismatch distributions (Librado & Rozas, 2009) and Harpending’s raggedness index (Harpending, 1994) (HRag), using 1,000 parametric bootstrap replicates. The formula t =τ/2u was used to estimate the time since expansion began for the expanding population. The value u was calculated as u = µkg, where µ is the substitution rate (6.5 × 10⁻⁹ synonymous substitutions per synonymous site per year (Gaut et al., 1996)), k is the average sequence length used for analysis (105 bp) and g is the generation time in years (1 year was used as an approximation for g according to age of first reproduction of C. dactylon). Fu’s FS (Fu, 1997) and Tajima’s D (Tajima, 1989) neutrality test were used to infer demographic history of each population. Statistical significance was estimated by performing 1,000 random permutations. Demographic analyses were performed in ARLEQUIN.

Genetic diversity and differentiation along longitude gradient

The genetic diversity parameters including the PIC value, He and I were calculated for EST-SSR markers in 13 populations (Table 2). The PIC values ranged from a high of 0.68 and a low of 0.40. Highest He value of 0.19 and lowest of 0.06 were calculated. The value of I ranged from 0.11 to 0.34 respectively for EST-SSR markers. Non-significant difference in genetic diversity parameters among groups (group 1 consisted of populations 1–3; group 2 consisted of populations 4-6; group 3 consisted of populations 7–9; group 4 consisted of populations 10–13) was observed along longitude gradients (He among populations: F = 0.027, P = 0.974). However, populations 4, 5, 6 and 8 at mid-longitude and populations 1, 12, 13 at lowest and highest longitude had higher genetic diversity than other populations according to I and He (Fig. 2). In general, there was also no significant correlation (P > 0.05) between intra-population genetic diversity and climate (mean annual temperature and precipitation), suggesting that temperature and precipitation showed little influence on intra-population genetic diversity along longitudes.

AMOVA results indicated that 38.23% of the genetic variation was attributed to difference among populations. About 61.77% of the genetic variation showed within population, which accounted for a greater part of the whole genetic variation (Table 3). Mantel’s test revealed significantly positive correlation between genetic distance and longitudinal distance (r = 0.276, P = 0.022). Environmental measures along longitude gradients (geographical distance (r = 0.438, P = 0.001); mean annual precipitation (r = 0.305, P = 0.033)) were correlated with genetic distance, while mean annual temperature was not (r = 0.006, P = 0.418). However, a partial Mantel test accounting for mean annual precipitation was used to estimate that genetic distance was not independently correlated with geographic distance (P = 0.146). Gene flow (Nm) among these C. dactylon populations along longitude gradients was 0.7701 (<1), which mainly reflected low communication frequency of gene among populations.

Genetic separation by UPGMA dendrogram and principal coordinate analysis

A UPGMA dendrogram dependent on Nei’s genetic distance showed that 13 different populations were grouped into three distinct clusters (Fig. 3). The first cluster consisted of populations from Lianyungang and Tancheng that located at high-longitude in the east of China near the Yellow Ocean. Most populations along longitude gradient assigned to the second cluster, so a certain degree of kinship was identified among populations of C. dactylon located at different regions along longitude gradient. Populations 7, 9, 10 clustered together and population 3 was close to population 5. The third mixture cluster contained plants coming exclusively from interval locations along the longitude, such as Shanxian (4), Zhengzhou (6), Sanmenxia (8), and Baoji (12). In addition, each of four populations in the mixture cluster had higher intra-population genetic diversity. The Principal coordinate analysis (PCoA) (Fig. 4) revealed a similar pattern as inferred with UPGMA dendrogram. The first, second and third principal coordinates respectively explained 31.57%, 17.78% and 11.90%, altogether 61.25% of the variability. At axis one, 13 populations could be divided into two parts with one of cluster 2 and others of clusters 1, 3. At axis two, cluster 1 consisting of populations 1 and 2 could be clearly separated from cluster 3. Additionally, populations of cluster 1 were clearly far from the other populations. Populations 3, 5, 7, 9 and 10 clustered together in the cluster 2. Populations 4, 6, 8 and 12 with relatively higher genetic diversity formed the cluster 3.

Population genetic structure of C. dactylon

Bayesian analysis of the population structure of C. dactylon was used to conduct separate calculation from K = 2 to K = 20. The results clearly demonstrated that the appropriate K values were obtained when the 13 populations were subdivided into 4 or 8 groups (Fig. S1). The results of 13 populations from K = 2–13 strongly indicated different clusters in populations of C. dactylon (Fig. 5; Fig. S2). When K = 4, population 1 and 2 clustered together (blue) and the cluster with green consisted of populations 4, 6, 8, 12. Populations 5, 7, 10, 13 were differentiated from populations 3, 9, 11 compared to K = 3. And when K = 8, one group included populations 3, 5, 7, 9, 10, and another group consisted of the populations 8 and 12. Each of the other six groups consisted of only one population respectively, which indicated a rich genetic differentiation along longitude gradients. Admixture among populations of medium longitude (populations 5, 6, 8, 9) could be observed from the structure of C. dactylon populations (K = 2 to 13), the clustering of populations at low and high longitude were relatively stable except for the populations 1, 12, 13.

Mismatch analysis and neutral test of C. dactylon

The mismatch distribution for each population of C. dactylon (Table 4) suggested populations 1–11, 13 possibly underwent recent population expansion. The result was also proved by significantly negative values of Fu’s FS among the 13 populations of C. dactylon. Significant SSD values and non-significantly positive Tajima’s D values of the population 12 might suggest relatively stable population size in Baoji. In addition, the demographic expansion time of each population was calculated based on the value of τ.

Genetic variation and differentiation among different ploidy levels

Genetic diversity (He and I) results of C. dactylon with five ploidy levels showed that triploid had significantly lower genetic diversity than pentaploid (P = 0.019 for He; P = 0.008 for I) and hexaploid (P = 0.049 for He; P = 0.020 for I) (Table S1). The genetic diversity of the C. dactylon individuals with higher ploidy level was significantly higher than that of individuals with lower ploidy level except for diploid. Significant genetic diversity exists among different ploidy level, but there is no significant difference in genetic diversity between pentaploid and hexaploid based on trend lines of PIC (P > 0.05). Diploid and triploid respectively had significantly lower genetic diversity than other polyploidy. All the individuals of C. dactylon at different ploidy levels fell into three major clusters based on UPGMA. The first and second clusters contained diploid and triploid individuals, respectively. Tetraploid, pentaploid, and hexaploid individuals formed the third cluster (Fig. 6). The genetic distance between two groups among different ploidy levels increased with the increase of ploidy level. In addition, there was no obvious growth trend of the distribution of ploidy levels, and regression analysis indicated that there was no significant linear association between ploidy level and longitude (r = 0.071; P > 0.05 –Fig. 7).

Genetic diversity and structure among populations along longitude gradients

Identifying relative roles of spatial and environmental factors in shaping patterns of population genetic differentiation is of interest for understanding landscape genetics. Each of 13 populations of C. dactylon along longitude gradients had a remarkable genetic variation within population based on the study. Genetic diversity can be improved by sexual propagation, and clonal propagation by rhizome and stolon can promote stable inheritance not being easily affected by genetic drift (Loveless & Hamrick, 1984). Intra-population genetic diversity had no obvious linear relationship with longitude in this study. In addition, mean annual temperature and precipitation at different longitudes had no significant influence on the intra-population genetic diversity. However, according to our previously published paper, higher within-population genetic diversity appeared at lower latitude with higher temperature (>25 °C), while relatively low temperature at higher latitude might have little effect on the intra-population genetic diversity (Zhang et al., 2019). So higher temperature (>25 °C) may be associated with higher genetic variation of C. dactylon. From sampling regions, there was a little higher genetic diversity of populations at mid-longitudes (4, 5, 6, 8) than other populations at low- and high-longitudes, because outcrossing and self-incompatibility of C. dactylon could produce low communication among different populations of the east–west transitional locations. Genetic differences among distant geographical sites with unique alleles can increase genetic diversity (Brochmann et al., 2003; Westergaard et al., 2011; Parducci et al., 2012; Eidesen et al., 2013).

Landscape heterogeneity may lead to the formation of local adaptive populations, and different genotypes may be selected to adapt to different environments (Felsenstein, 1976; Rieseberg & Willis, 2007; Anderson, Willis & Mitchell-Olds, 2011). Environmental distance matrices were also calculated to estimate the influences of environmental distance matrices (IBE) on the genetic distance. No significant correlation was observed between genetic distance and geographic distance after controlling the effect of mean annual precipitation, which might be due to a correlation between geographic distance and mean annual precipitation. P. arundinacea (Phalaris) also exhibits a certain broad-scale population genetic structure in east-to-west gradients (Perdereau et al., 2017). The adaptive differentiation among populations is influenced by the natural selection of heterogeneous landscape and gene flow (Dionne et al., 2008; Poelchau & Hamrick, 2012). Low east-to-west inter-population gene flow (Nm = 0.7701) indicated that there was little adaptive selective pressure mediated evolutionary dynamics among different populations along longitude gradients (Casler et al., 2007). In addition, these factors to hinder further communication among populations include the transmission mechanism of pollen and seeds of species and the barrier caused by environmental conditions. In our study, restricted seed dispersal may exist among populations at different longitudes exhibiting significant IBD, and mountains in the Henan–Shanxi region may hinder further genetic communication among populations as natural barriers. Alpine plants from several mountain ranges in Europe usually have a strong patterns of longitudinal differentiation (Ronikier, Cieślak & Korbecka, 2008; Ronikier, Schneeweiss & Schönswetter, 2012). Geographical isolation along an ecological gradient or environmental isolation can lead to strong genetic differentiation for widespread plants (Hoskin et al., 2005).

The result of AMOVA indicated 61.77% of the genetic variation was attributed to differences within populations and FST was 0.3823. Higher FST along with increasing distance along longitude than latitude gradients may indicate local genetic differentiation among populations along longitude gradients (Edillo et al., 2009). Both spatial and environmental factors interweave in the formation of population differentiation pattern. In our case, low gene flow among neighboring populations may show the isolation-by-distance pattern among two neighboring populations. Several reproductive barriers in different environments within the same geographical area can block gene flow (Lowry, 2012; Fishman et al., 2013). Interestingly, some populations separated by 200 km in distance are close in genetic distance. A probable explanation is that human-mediated introduction allow contact among these populations. The mismatch distribution test and neutrality test both indicated expansion within population along longitude gradients, with an estimated time around 12.5 Ma (under the demographic expansion model), which may be related to forests shrink and grasslands spread during Miocene time. The dramatic expansion of C4 grasslands in India is a part of increasing tropical C4 savannah vegetation worldwide in the late Miocene, which has been variously ascribed to uplift of the Tibetan plateau (Polissar et al., 2015). In contrast, the population in Baoji with non-significant SSD and HRag values was interpreted as reflecting past demographic stability.

The low number of admixed genotypes at mid-longitude can be explained by the geographical distribution of the collection sites. Natural or human-mediated transmission may have increased genetic diversity and affected the genetic structure of several grass species (Balfourier, Imbert & Charmet, 2000; McGrath, Hodkinson & Barth, 2007; Petit et al., 2003; Taberlet et al., 1998). Hybridization at population boundaries and gene migration appear to promote homogenization. Populations at low and high-longitude may evolve much earlier than populations from sites of mid-longitude based on the result. Therefore, populations of C. dactylon at mid-longitude may evolve from complex introgressions between populations at low and high-longitude with close genetic affiliation. Historic changes in landscape configuration, spatial isolation and related gene flow changes may gradually lead to intraspecific genetic variation between eastern and western populations of C. dactylon in different geographical areas. The admixture among populations determines their potential to adapt to new environmental conditions (Lee, 2002).

Relationship between ploidy levels and genetic diversity

The higher values of He, I and PIC were observed in populations with pentaploid and hexaploid than triploid. Polyploid retains most of the population’s gene pool in the form of fixed heterozygosity (Brochmann et al., 2004). Hybridization of polyploid progeny from different parents of the same species can increase the genetic diversity. Multiple introductions from different origins and recombination between different genomes led to genetic mixture. Moreover, groups of individuals with the same ploidy level were separated by genetic distance along with the increasing ploidy levels. The higher genetic distance between two different ploidy levels along with the increasing ploidy levels was observed in the study, which indicated that each chromosome replication process might gradually increase the genetic dissimilarity. Diploid plants can participate in the formation of many polyploid (Ramsey & Schemske, 2002; Otto, 2007), and polyploid are often separated from their progenitors to produce new genotypes (Čertner et al., 2017; Alsayied et al., 2015; Sweigart, Martin & Willis, 2008). Different chromosome numbers could lead to reproductive segregation or sterility, which limited gene flow between hybrids and their parents (Tate, Soltis & Soltis, 2005). Polyploidy not only produce direct genetic redundancy, but also promote emergence of new gene functions and drive diversification and evolution (Alix et al., 2017). The hybrid genome composition of two related heterotetraploid species was analyzed using molecular cytogenetics to elucidate the role of heteropolyploid in the realization of Hordeum diversity (Cuadrado, Bustos & Jouve, 2017). From the facts above, polyploidization can lead to rich genetic variation and diversification for probable speciation of C. dactylon.

Comparison of genetic pattern between longitude and latitude

No significant difference exists in genetic diversity parameters within-population along longitude gradients, while higher within-population genetic diversity appeares at low latitude with higher temperature along latitudinal gradients (>25 °C) (Zhang et al., 2019). The result showed higher temperature might promote higher within-population genetic variation of C. dactylon. The barrier of hybridization of rice mainly focused on the environmental variation, just as temperature influenced cellularization of the endosperm to promote hybrid (Folsom et al., 2014). Higher temperatures may have affected the genetic similarity of Culex pipiens populations sampled early and later in the breeding season (Edillo et al., 2009).

A relatively low gene flow (Nm = 0.7701) means a rich genetic differentiation along longitude gradients (Nm < 1) (Wright, 1965). In our previous study on the genetic diversity and genetic structure of C. dactylon along a latitudinal gradient, low genetic differentiation and notable admixture structure existed among populations along latitudes (Zhang et al., 2019). Low genetic differentiation caused by different latitudes may be attributed to the wind-pollination outcrossing system and vegetative materials through human activities, which indicated that the latitude interval could not be an obstacle to genetic exchange. In addition, stronger admixture could be observed along with the latitude in another study (Zhang et al., 2019) than that along with longitude in this study, which indicated relatively lower gene flow among populations of C. dactylon occurred along the longitude gradients. North-south direction gene flow has existed between the eastern Lemus chinense populations and the western ones (Long et al., 2019). The eastern Eurasian populations of some species with global distribution pattern are a little different from-the western Eurasian populations, which occurs between late Miocene and Pleistocene (Tabata et al., 2018; Hardouin et al., 2015; Rodrigues et al., 2014). East Asia was less extensively glaciated and had a relatively mild glacial climate. The Qinghai-Tibetan Plateau uplifted rapidly in Pleistocene and the Asian summer monsoon intensified, so the formation of the inland arid belt of Eurasia resulted in long-term ecological isolation (Clark & Mix, 2002; López-Pujol et al., 2011). A probable explanation is that the genetic differentiation along longitude gradients was little and the gene flow was frequent before the uplift of Tibet Plateau, but the situation reversed after that because populations were ecologically isolated. Therefore making a comparison between latitude and longitude, a stronger gene flow from the north-south than east–west direction is characteristic of genetic structure of C. dactylon in China. Distinct genetic divergence of C. dactylon presented between the lower and the higher ploidy level along both longitude and latitude gradients. In addition, widespread polyploidization may cause significantly increasing genetic diversity with the ploidy level.