Within-species genetic diversity is shaped by multiple evolutionary forces within the confines of geography, and has cascading effects on the biodiversity of other taxa and levels. Invasive species are often initially limited in genetic diversity but still respond rapidly to their new range, possibly through ‘pre-adapted’ genotypes or multiple sources of genetic diversity, but little is known about how their genetic structure differs from that of native species and how it alters the genetic-species diversity relationship. Here, we selected a widespread native species (Phragmites australis) and its co-occurring invasive competitor (Spartina alterniflora) as our model plant species. We investigated the genetic structure of P. australis using two chloroplast fragments and ten nuclear microsatellites in 13 populations along the Chinese coastal wetlands. We discovered a distinct geographical differentiation, showing that the northern and southern populations harbored unique genotypes. We also found a significant increase in genetic diversity (allelic richness and expected heterozygosity) from south to north. Combined with previous studies of S. alterniflora, the Mantel tests revealed a significant correlation of genetic distances between P. australis and S. alterniflora even when controlling for geographic distance, suggesting that the invasive species S. alterniflora might exhibit a phylogeographic pattern similar to that of the native species to some extent. Furthermore, our results suggest that the S. alterniflora invasion has altered the relationship between the genetic diversity of the dominant native plant and the associated species richness of soil nematodes. The reason for the alteration of genetic-species diversity relationship might be that the biological invasion weakens the environmental impact on both levels of biodiversity. Our findings contribute to understanding the latitudinal patterns of intraspecific genetic diversity in widespread species. This work on the genetic diversity analysis of native species also provides significant implications for the invasion stage and ecological consequences of biological invasions.

Sampling, genotyping and sequencing of P. australis

We collected 194 individuals of P. australis from 13 sites along the Chinese coast. To avoid collecting the same clone and to get enough genetic variation, we selected five P. australis reed populations within each site, with a distance of at least 1 km between stands. We collected three individuals from each population. All individuals were transplanted with rhizomes in a common garden at the Jiangwan Campus of Fudan University in Shanghai (31.28°N, 121.48°E). After the plants regrew, we collected and dried young leaves from all individuals, and stored them in zip-lock plastic bags with silica gel at room temperature until DNA isolation. We extracted total DNA from the dried leaves according to a modified cetyltrimethylammonium bromide (CTAB) method. We examined the quality and quantity of extracted DNA with 1% agarose gels and a microscope spectrophotometer, and stored DNA at -20℃ until later genotyping and sequencing.

To measure genetic variation, we used 10 microsatellite primer pairs previously designed for P. australis (Saltonstall 2003, Yu et al. 2013). Forward primers were labeled at the 5’ end with the fluorescent dyes FAM, HEX or TAMRA. We performed polymerase chain reaction (PCR) as described by Liu et al. (2022), and separated the PCR products by capillary electrophoresis using an ABI 3730XL DNA capillary sequencer (Applied Biosystems, Foster City, California, USA) after confirming the PCR product on a 2% agarose gel. We scored fragment profiles and carefully check the stutter peaks and the low-frequency alleles with GeneMarker 2.2.0 to reduce the potential effect of null allele. We did not discover the null alleles with the Hardy-Weinberg equilibrium-based method, since there is no reliable approach to elimination of allele dosage for our polyploid data. The same clones were detected by the function assignClones in R package polysat (Clark and Jasieniuk 2011). The duplicated genotypes were removed for further genetic estimates.

To determine the haplotype, we amplified two non-coding chloroplast regions by PCR in one sample of each stand, using the primer pairs [trnT (UGU) “a”-trnL (UAA) “b” and rbcL-psaI] as described previously (Saltonstall 2002). We sequenced the PCR products in both directions on an ABI 3730XL DNA sequencer (Applied Biosystems). We assembled and checked the sequencing with SeqMan 7.7.0 (Lasergene, Santa Clara, USA) and identified haplotypes to the naming scheme of P. australis described by Saltonstall (Saltonstall 2016).

Data analysis of genetic diversity and structure of P. australis

To estimate the genetic diversity level of P. australis, we calculated the number of alleles per locus or allelic richness (N_a), and the expected heterozygosity (H_e) with R package polysat (Clark and Jasieniuk 2011). We assessed the relationship between genetic diversity and latitude using linear regression.

To assess the genetic structure of P. australis, we calculated genetic differentiation (F_st) with the R package polysat. We also calculated Pairwise Bruvo distances based on microsatellite variation, and used the genetic distance matrix for principal coordinates analysis (PCoA) and hierarchical cluster analysis using the unweighted pair-group method with arithmetic means (UPGMA). We applied Bayesian clustering with Structure 2.3.4 (Pritchard et al. 2000) to detect the genetic structure of P. australis. We performed 20 replicates of the clustering analysis at each value of K from 1 to 10 under the admixture model with 50,000 burn-in steps and 500,000 Markov Chain Monte Carlo repeats. We calculated Delta K using the online program Structure Harvest (Earl and vonHoldt 2012) to determine the most likely cluster number (K value) for our genetic data, grouped replicates in CLUMPP 1.1.2b (Jakobsson and Rosenberg 2007) and visualized in DISTRUCT 1.1 (Ramasamy et al. 2014).

Correlation analysis between geographical and genetic distances of P. australis and S. alterniflora

We used the previously published nuclear microsatellites and chloroplast sequences data of S. alterniflora in China (Qiao et al. 2019, Shang et al. 2019, Xia et al. 2020). We extracted the geographical coordinates and the diversity indices (i.e., Allele number, N_a; Expected heterozygosity, H_e) of surveyed populations of S. alterniflora from Shang et al. (2019) and Xia et al. (2020) for further comparisons. For genetic analysis of S. alterniflora, we used the raw data of 11 nuclear microsatellites from Qiao et al. (2019). We removed three loci from the raw dataset because there were many missing values or null alleles in loci 5, 7 and 9. We calculated geographic distances using the function distm in R package geosphere and used pairwise Fst for genetic distance.

We used Mantel test and multiple matrix regression with randomization (MMRR) (Wang 2013) to examine relationships between geographic and genetic distance matrices at the site level. We ran correlation analyses between geographical and genetic distance matrices using the function mantel in R package vegan, and regression analyses using the function MMRR written by Wang (2013) with genetic distance as the dependent matrix and geographical distances as the independent (predictor) matrices with 9,999 permutations. The correlation of genetic distances between the two species were also performed with partial Mantel test while controlling the geographical distance for seven common sites.

Correlation between genetic variation and nematode community

We used the geographic records of nematode genera from a published work (Zhang et al. 2019). These nematode data were investigated to reveal the biotic homogenization of nematode communities by exotic S. alterniflora in China. This study found a clear latitudinal cline (nematode diversity increased with increasing latitude) and a strong correlation of nematode diversity to environmental variables in soils for P. australis, but weak for S. alterniflora (Zhang et al. 2019). Because the TJ and TS sites were located within a very short distance (approximately 54 km) around the same bay, we considered them as one site when comparing the variation in geography, genetics, and community. Thus, we had seven common sites with both genetic and nematode information. We estimated the Jaccard distances between nematode communities using the function dist with a method binary parameter. We used these Jaccard distances to perform principal coordinates analysis (PCoA) of nematodes. Matrix correlation analyses between geographic, genetic and nematode distance matrices for seven common sites were performed using both Mantel test using the function mantel in R package vegan and MMRR using the function MMRR written by Wang (2013) with nematode distance matrices as the dependent variable using 9,999 permutations.