Introduction

Estuaries are dynamic environments with frequent disturbance and extreme conditions, which provide repeated opportunities for the decline and the establishment of new populations. Propagules can be recruited from nearby populations or be introduced from distant and different environments by migratory birds, or international ship traffic. One of the most common and productive aquatic plants inhabiting the world’s estuaries is Phragmites australis (Cav.) Trin. ex Steud., a cosmopolitan tall grass with high intraspecific variation and ecological amplitude1. Phragmites australis populations are extremely phylogenetically diverse in the Mississippi River Delta2,3 and in the Danube Delta4. Phylogenetic variation, including ploidy variation2,3,4, provides fitness in these dynamic environments due to the enlarged eco-physiological adaptability provided by a gene pool of genotypes of multiple distinct origins5,6,7. Due to their genetic diversity and the strong and variable selection pressure of the environment, estuarine marshlands are evolutionary hotspots of Phragmites diversity and can be dangerous gateways for invasions, especially cryptic invasions, in ranges where P. australis is native. For example, P. australis invasion in North America started at New York Harbor in marshes of the Hudson River estuary, where introduced European P. australis displaced the native populations of Phragmites australis ssp. americanus8,9.

Salt tolerance is a competitive trait in brackish estuarine environments. As a species, P. australis tolerates a wide salinity range10,11,12,13,14, but genotypes of different phylogeographic origin differ in salinity tolerance6,15 due to their different bioclimatic origins. For example, the invasive Mediterranean Phragmites lineage thrives in the brackish Mississippi River Delta thanks to its adaptation to draught evolved in the warm, dry climate of its native range in the Mediterranean and Middle East5,7. In the Yellow River Delta in China the genetic diversity of P. australis decreases with increasing soil salinity and only specific allelic phenotypes occur in the saltiest patches of the mosaic of saline habitats in the delta16. Such distinctive genotypes may have evolved locally from native populations as well as from introduced pre-adapted genotypes.

In this study we investigate phylogenetic diversity and its role in salt adaptation in three estuarine P. australis populations in China, in the Yangtze, Yellow River and Liaohe deltas. Unlike previous local studies of P. australis diversity in Chinese coastal marshes, we analysed the genetic diversity patterns of the Chinese populations within the global phylogeographic structure of the genus Phragmites17,18 in order to resolve their phylogeographic relationships and trace possible introductions. Several Phragmites lineages have recently been identified in China and Korea19,20,21 and classified according to the Phragmites classification system introduced by Saltonstall based on chloroplast DNA sequences9,22, however a macro-perspective of the evolutionary history of the exceptionally high diversity of Phragmites lineages in East Asia is lacking. There are at least three Phragmites species in East Asia: P. japonicus Steud., P. karka (Retz.) Trin. ex Steud.and P. australis, with ploidy variation from 2n = 4 × to 10x (with 2n = 8 × being dominant23,24), and several salt-tolerant ecotypes10,16,25. Four P. australis ecotypes have been phenotyped in the Yellow River Delta based on tissue Na+/K+ ratios as an indicator of salt stress10. Takahashi et al.25 showed that salt-tolerant reeds from the Yellow River Delta contained low Na+ and high K+ content compared to salt-sensitive plants from a freshwater river in Japan which contained high Na+ and low K+ when exposed to salt. This suggests that salt tolerant ecotypes in the Yellow River may be able to select K+ over Na+ in saline conditions and maintain ion homeostasis25. The HKT1 gene functions as a K+/Na+ co-transporter in P. australis and an insertion in an HKT1 intron causes alternative splicing of the mRNA into the two transcripts of the salt-tolerant and salt-sensitive reeds25. Variation in the genomic sequence of the HKT1 gene, as documented by25 and by20 in wild Korean populations, may therefore distinguish salt-adapted and non-adapted genotypes25. Therefore, in addition to (1) the phylogenetic relationships, we investigated also the (2) genetic structure of Phragmites populations (inferred by SSR—or microsatellites) in relation to soil salinity, (3) the HKT1 gene polymorphism and (4) the phylogenetic diversity of chloroplast DNA. Based on previous studies of estuarine Phragmites populations elsewhere, we expected high phylogenetic diversity in the study area and a salt-dependent distribution of the phylo-types of different origin.

Results

Phylogenetic relationships of the estuarine populations

Fragment size polymorphism in the TrnT-TrnL short region of the cpDNA revealed two haplotype profiles within the Chinese estuarine populations (Fig. 1, Supplementary Information S1). Following Phragmites classification system, the sequences matched P. australis Haplotype P and O and their intrahaplotypic microsatellite variants (P1, P2, P3, P4 and O2, O3) or differed from Haplotype P in one to two substitutions (new Haplotypes AS, AT, AU) (GenBank accessions no. KP994324 to KP994334). In total we found 9 different haplotypes of two distantly related phylogeographic groups (Fig. 2) previously defined as the Far East-Australian group (hereafter P-related genotypes) and “European” Phragmites, despite its almost cosmopolitan distribution (hereafter O-related genotypes). Compared to the previous phylogenies the inclusion of the new haplotypes from Asia and Australia moved the Asian/Australian group closer to the tropical species P. mauritianus and P. karka (note the haplotypes of P. karka marked blue in the lower part of the PCoA with 71% support, Fig. 2). The populations in the Yangtze River Delta (CML, JDSL, CXL and SJL) were entirely composed of Haplotype P-related genotypes, whereas the population in the Yellow River Delta consisted of O-related genotypes (with the exception of one single P-related genotype). The Tianjin population consisted of O-related genotypes and the Liaohe population was mixed with O- and P-related genotypes co-occurring in close proximity (Table S1).

Figure 1
figure 1

Map of the Chinese populations sampled in the study (ArcGIS ver. 10.6, Environmental Systems Research Institute, ESRI).

Full size image
Figure 2
figure 2

Phylogenetic position of the Chinese populations (purple) in the global phylogeographic structure of the genus Phragmites, inferred by cpDNA sequences (trnT-trnL and rbcL-psaI). The haplotypes previously identified in East Asia and Australia are in red. The coordinates account for 44% of the variation in the data (27% coord. 1 and 17% coord. 2). The circles with the numbers indicate the statistical support from the parsimony analiysis. Haplotypes IDs follow the names in GenBank. Haplotypes source: H28, H29, H30, H31, H32 and H33 (An et al. 2012). E4, S2 (corresponding to haplotype P1 in our study) and Phragmites japonicus (labeled Pj) (Chu et al., 2011). SLJ01 (Hurry et al. 2013). J, O, P, Q, X, Y (Saltonstall 2002). AN (Lambertini et al. unpublished). P1, P2, P3, P4, P5 and O2, O3 have been classified as microsatellites variants of haplotype P and O, as they differ from these in the numbers of repeats at the microsatellite loci. Statistical support is from the parsimony analysis. The PcoA is made with the program GenAlex ver. 6.445.

Full size image

The Bayesian analysis of SSR data inferred three ancestral populations for the Chinese estuarine populations (Fig. 3): one for the Haplotype-O related genotypes and P. japonicus (Haplotype AM) (orange), one for the Haplotype P-related genotypes (light blue), and one for the hybrids of Haplotype P and P. karka (dark blue). This analysis revealed several hybrids in the Yangtze Delta populations between P. australis P-related genotypes and a population of P. karka in the Mekong Delta in Viet Nam used as a P. karka reference in this study (Haplotypes I and U and their microsatellite variants), and in the Liaohe Delta population between P. australis Haplotypes P- and O- related genotypes and between P. australis Haplotype P and P. karka. Both hybrid types were Haplotype P. Only one hybrid in the Liaohe Delta population was Haplotype O and had admixed ancestry (higher than 20%) with both Haplotype P and P. karka (Fig. 3).

Figure 3
figure 3

Population structure of SSR data. Inferred ancestry probability is in Y axis and population ID in X axis. (a) Population structure for K = 3. (b) Evanno’s “deltaK” inferring the most likely number of ancestral populations (K). The graph is made with Structure ver. 2.3.447 and Evanno’s delta K graph is made with Structure Harvester48.

Full size image

Population genetic structure

The SSRs divided the sample set into two main groups, which reflected the haplotypic structure of the populations better than their geographic distribution (Fig. 4a,b), and the structure of HKT1 homo-and heterozygotes (Fig. 4c,d).

Figure 4
figure 4

Principal coordinate analysis of SSR pairwise genetic distances. The coordinates account for 37% of the variation (14% coord. 1 and 9% coord. 2). (a) the different colours indicate genotypes of populations of different geographic locations. (b) The different colours indicate genotypes related to Haplotype O and Haplotype P inferred by the trnT short fragment. (c) The different colours indicate homozygotes and heterozygotes at the HKT1-1 gene locus. (d) The different colours indicate homozygotes and heterozygotes at the HKT1-2 gene locus. The PcoAs are made with the program GenAlex ver. 6.445.

Full size image

Significant genetic divergence among the seven populations (Supplementary Information S1) was detected by the AMOVA of SSR data (Table 1a). The extent of divergence increased when the variance was tested between Haplotype O- and P- related genotypes (Table 1b) and disappeared among the populations related to Haplotype O when population structure was tested within haplotypes (Table 1c) and within population (Table 1d). It remained significant among the populations related to Haplotype P, both among populations of different deltas (Table 1e) and among the populations in the Yangtze Delta (CML, CXL and JDSL) (Table 1f). Differences in genetic variance due to SSR alleles were also detected between homozygotes and heterozygotes at the HKT1-1 locus of the HKT1 gene (Table 1g) (but not at the HKT1-2 locus, Table 1h), however such a distinction disappeared when the variance was tested between homo- and heterozygotes of Haplotype O (Table 1i,j). It remained significant between the homo- and heterozygotes of the Haplotype P-related genotypes both at the HKT1-1 and HKT1-2 loci (Table 1k,l).

Table 1 AMOVA of SSR data testing the structure within and among populations based on (1a) the geographic distribution of the populations, (2b–f) the haplotypic composition of the populations and (2g–l) the HKT1 composition of the populations in homo- and heterozygotes.
Full size table

Three to five alleles per locus were frequent in the SSR profiles of both the Liaohe and Yangtze Haplotype P-related genotypes, whereas Haplotype-O related genotypes showed a disomic inheritance pattern. However, two exceptions of genotypes with 3 alleles at locus Pagt4 were found also in the O- related populations in the Yellow River and Liaohe Deltas. Genetic diversity (I, h and uh) in SSRs was highest within the Haplotype-P related populations due to the higher number of alleles in P- related than in O-related genotypes (Table 2).

Table 2 Genetic diversity in haplotype O- and haplotype P-related populations.
Full size table

Genetic diversity distribution in relation to soil salinity

The factor that best explained the distribution of the genotypes in relation to salt was phylogenetic diversity, considering the four independent phylo-types: Haplotype O- and P-related genotypes and their hybrids (Haplotype P x P. karka and Haplotype P x Haplotype O) (P = 0.01). Genotypes related to Haplotype O were distributed in areas with higher salinity than those related to Haplotype P. Both hybrid types of Haplotype P had a higher variance in salt tolerance than Haplotype O- and P-related genotypes, and occurred in areas of intermediate salinity, compared to the two haplotypes (Fig. 5).

Figure 5
figure 5

Distribution of phylo-types (X axis) by the salinity classes of the SQ5 layer of excess salt (Y axis). Average salt excess classes values and Bonferroni’s intervals for each phylo-type. (A) (AB) and (B) refer to the significance of the Bonferroni multiple comparison. The graph is produced in R (R core Team).

Full size image

The distribution of HKT1-homozygotes and heterozygotes at both HKT1-1 and HKT1-2 gene loci was not affected by salinity (P = 0.87 for HKT1-1; P = 0.50 for HKT1-2). The HKT1 sequences of the homozygotes showed large variation, especially among the genotypes of Haplotype O. The sequences of the previously documented salt-tolerant genotype from Nanpi in the Yellow River Delta in GenBank matched the sequence of some of our Haplotype O genotypes in the Yellow River Delta. The sequences of the salt-sensitive genotype from Ustonia in Japan (also in GenBank) matched instead the sequences of our P. japonicus references in South Korea and Sakhalin Island (Russia), and was different from the sequences of Haplotype P-related genotypes in the Yangtze River Delta. The sequences of P. karka genotypes were different from those of P. australis Haplotype O and Haplotype P, as well as from those of P. japonicus. We found correlations neither between SSR allelic frequencies and the salinity classes of the Soil Quality layer 5 (SQ5) of the Harmonized World Soil Database, nor with our salinity measurements in the Yangtze Delta, nor between allelic and HKT1 homo/heterozygotes frequencies, nor between HKT1 homo/heterozygotes frequencies and salinity.

Discussion

Phylogenetic relationships of the estuarine populations

Haplotypes of two phylogeographic regions of P. australis co-occur in the coastal populations in East China. The genotypes related to Haplotype O are closely related to the European populations of P. australis and the East Asian endemic species P. japonicus, whereas Haplotypes P, AS, AT, AU are closely related to populations of P. australis in Australia. In this study they appeared closely related also to the tropical species P. karka in Asia and P. mauritianus in Africa. Haplotypes O and P are the dominant haplotypes in the studied populations and both have a wide distribution range. Halpotype O is found across temperate Europe9,18 and Asia19 and its cp-microsatellite variant O2 found in this study is shared by the two continents. Haplotype P has previously been found in East Asia and Australia9,10,15,18,19 and its cp-microsatellite P4 was also found in both continents. The other novel Haplotypes AS, AT and AU appeared to have evolved locally from Haplotype P- related genotypes, given their restricted ranges within the studied populations and the very few base substitutions compared to the sequences of Haplotype P-genotypes analysed in this study. Two populations were identified also by Gao et al.16 to explain the structure of the Yellow River population, and two ancestral populations were discovered for Korean Phragmites20. Gao et al.16 justified the two populations as salt-adapted and non-adapted genotypes, while Chu et al.20 explained the two ancestries as due to two Phragmites species, P. australis (sequenced as Haplotype P in that study) and P. japonicus in Korea. Two branches for Phragmites australis in East Asia were detected also by Yao et al.26 by ITS markers, that were identified as a suitable DNA barcode to capture the high genetic diversity of halo-tolerant Poaceae species in coastal areas.

In our study, Haplotype O-related genotypes occurred in areas with the highest salt impact and Haplotype P-related genotypes in areas with the lowest. As our samples from the Yellow River Delta were all Haplotype O-related genotypes, except for a single Haplotype P, and given our different geographic sampling from that of16, we cannot relate our results to those of their study directly. However, we also found a salt-dependent distribution of our samples that seems to be explained by genotypes of different phylogenetic origin which, in agreement with16, have distinct allelic patterns. In agreement with Chu et al.20, we found two distinct populations for P. japonicus and P. australis Haplotype P. As in our previous studies, P. japonicus is closely related to P. australis Haplotype O-related genotypes, both in cpDNA sequences and SSRs alleles18. In the present study Haplotype P appeared distantly related to its conspecific Haplotype O and more closely related to P. karka, with which it is also introgressed. Phragmites karka is known to occur and hybridize with P. australis haplotype P in southwestern China21. However, in contrast to21, P. australis Haplotype P is the seed donor of the hybrids that we found in the Yangtze and Liaohe deltas. Given the high degree of introgression and the lack of P. karka in our samples from the populations where the hybrids occur, P. karka introgression appears an inherited trait of Haplotype P rather than a recent local hybridization event.

The polysomic profiles of the SSRs alleles indicate a higher ploidy level for Haplotype P-related genotypes that consistently had more than two alleles, than that of Haplotype O-related genotypes that, instead, had a disomic pattern. Several studies have shown that Asian P. australis is dominated by octoploids and to a minor extent by hexaploids and decaploids, and that tetraploids are less frequent17,23,24. Based on our studies17,18 that sequenced the samples previously analysed for genome size by24, and our unpublished chromosome counts of Haplotype P and Haplotype O genotypes by “squash and stain”, Haplotype P-related genotypes are more frequently octoploids, whereas Haplotype O-related genotypes are more likely tetraploids, according to the sample sets that we studied.

The higher genetic diversity of Haplotype P octoploid populations (about the double of that of Haplotype O), and the high degree of introgression, suggests that Haplotype P is an allopolyploid that originated by hybridization between P. karka and another still unknown population different from Haplotype O and P. japonicus. Haplotype P shares in fact a large part of its genome with P. karka, but not with P. australis Haplotype O or P. japonicus, and the cpDNA is closely related to that of P. karka. Ishi and Kadono27 confirmed that octoploid P. australis is fertile in Asia. They found anomalies in neither pollen viability, nor availability, nor barriers to cross-pollination, nor reduced seed set rates in wild octoploid populations in Japan. The high genetic diversity that we found in our Haplotype P populations and the occurrence of hybrids between Haplotype P and Haplotype O in the Liaohe Delta further confirm that Haplotype P-related genotypes can reproduce sexually, and indicate that ploidy level is not a barrier to gene flow between the two Haplotypes. All hybrids between Haplotype P- and Haplotype O- related genotypes, except one, had the cpDNA of Haplotype P, indicating a higher ability of the octoploid genome to be recombined, but also a recent contact between Haplotypes P and O, as obvious hybrids were found only in the Liaohe populations where genotypes of the two lineages coexist. The Liaohe delta is one of the largest reed stands in the world and is managed for pulp and paper production28. Haplotype P could have been introduced here, and this could explain the co-dominance of the two distantly related haplotypes only in this delta.

The genetic diversity of Haplotype O-related populations was also relatively high and comparable to that of the European P. australis tetraploid population and higher than the diversity of the invasive European populations in North America and in the Gulf Coast3. Although this can rule out founder effects due to a recent introduction of Haplotype O to East Asia, the sympatry of the two distantly P. australis related lineages and in the same habitat, and their relationships to the two Asian species P. karka and P. japonicus, pose several questions concerning P. australis evolution in East Asia that cannot be explained by natural selection only. These results have therefore implications also for Phragmites taxonomy, systematics and ecology in East Asia and are central for future research.

Population genetic structure

The populations of the two haplotypes also differed in genetic structure. Insignificant population differentiation measures (PhiPTs) confirm gene flow and long-distance dispersal of pollen and seeds among the populations related to Haplotype O, whereas the significant PhiPTs among the populations related to Haplotype P indicate differentiation among populations even at the local scale, as seen among the populations in the Yangtze Delta. Nevertheless, the low PhiPT values (max 0.16) support long-distance gene flow even among the Haplotype P-related populations. This significant structure could be due to the high genetic diversity dispersed by the seeds of an octoploid allopolyploid. Lack of a regional geographic structure is a common feature of P. australis populations in Europe29,30,31, North America32 and South Africa33. Even in Australia, Haplotype P populations are not genetically distinct15. The significantly different genetic variation pattern of Haplotype P populations in East Asia is therefore unusual for P. australis. Previous studies of P. australis intraspecific variation in China have found similar patterns and attributed the spatial differentiation to edaphic factors in the Songnen Prairie34,35, swamp, salt and dune reed ecotypes in the Hexi Corridor in the Gansu province36, and salt tolerance in the Yellow River Delta11,16 . None of these studies, however, resolved or considered the phylogenetic diversity of the populations, or realized that population structure and ecology are different within lineages. Resolving phylogenetic relationships appears therefore crucial from this study to address ecology and adaptation in future Phragmites studies in East Asia and elsewhere.

Genetic diversity distribution in relation to salinity

The HKT1 gene encodes one of the main mechanisms of salt tolerance in vascular plants37. It varies in the genomic sequences of P. australis in East Asia20,25 and there is allelic variation associated with salt response in wild coastal populations of Arabidopsis thaliana in Northern Europe38. In our study salinity did not affect the distribution of homo- and heterozygotes at the HKT1-1 and HKT1-2 sites of the HKT1 gene, but the factor that better explained the relationship between salinity and genetic diversity, as well as the HKT1 variation pattern, was phylogenetic origin. Haplotype O-related genotypes occurred in areas with the highest salinity and Haplotype P-related genotypes in areas with the lowest. Both Haplotype P hybrid types (Haplotype P x P. karka and Haplotype P x Haplotype O) occurred in a wider range of salinity than Haplotypes O and P, and had higher average salt tolerance than Haplotype P. The P. karka population in the present study is from the Mekong Delta in Viet Nam, another estuarine area subject to increasing salinization39. Salt tolerance as a specific trait of P. karka remains poorly understood, but our results show that hybridization between estuarine populations of P. karka, and Haplotype O, increased salt tolerance in Haplotype P in East China.

Compared to25 sequences, the pattern of genetic variation of the complete HKT1 gene sequence of the homozygotes appeared more variable and complex than that of two distinct sequences for salt-tolerant and salt-sensitive genotypes. As could be expected, the coding sequence of salt-sensitive Haplotype P was different from that of salt-tolerant P. karka. Likewise, the coding sequence of salt-tolerant Haplotype O (matching the salt-tolerant coding sequence of25) was distinct from that of its close relative P. japonicus (matching the salt-sensitive sequence of25), suggesting that salt tolerance has evolved multiple times and independently in East China. The sequences of the hybrids could not be resolved with certainty because of heterozygosity at several loci.

Conclusions

As we have shown previously5,6,7, phylogenetic relationships play a crucial adaptive role for Phragmites populations. Two distantly related P. australis haplotypes dominate the estuarine populations in East China. One, Haplotype O, is related to European Phragmites sensu Lambertini et al.18, the other, Haplotype P, is highly introgressed with P. karka and has its distribution in East Asia and Australia. The two haplotypes have different salt tolerance, with Haplotype O occurring in the most saline areas in the deltas and Haplotype P in the lowest. Hybridization with haplotype O and introgression with estuarine populations of P. karka has increased salt tolerance in Haplotype P and expanded its distribution range to the most saline areas in the Yangtze and Liaohe deltas. Saline tolerance has therefore evolved multiple times and independently in Haplotype P and Haplotype O. Interestingly, salt-sensitive Haplotype P is closely related to salt-tolerant Asian P. karka, whereas salt-tolerant Haplotype O is closely related to Asian P. japonicus, whose HKT1 sequences matched those of the salt-sensitive genotypes of25. Future research may address the relationships of the two P. australis haplotypes with the two Asian species as well as the evolutionary and ecological significance of the two P. australis haplotypes in the range of the Asian species.

The present study has found no evidence of a cryptic invasion in East China, but cannot even resolve conclusively why two P. australis haplotypes, each with an intercontinental distribution range, do co-occur in the estuarine populations in East Asia and share the same habitat, i.e. whether East Asia is their center of diversity or they (one or both) have been introduced recently. The high genetic diversity in the populations of both haplotypes suggests an origin in East Asia, whereas the restricted distribution of their viable hybrids only to sympatry areas suggests a recent contact, or recent changes in the latitudinal distribution of the two haplotypes. The polysomic allelic pattern of Haplotype P further suggests a recent polyploidization. A larger sample of Haplotype O populations from across Europe and Asia and of Haplotype P from Australia and the Asian tropical region is needed to solve this enigma. Herbarium specimens may also help reconstruct past distribution ranges and can detect recent changes.

Methods

Sample set

We collected 174 samples of P. australis from 7 populations in East China in the Yangtze, Yellow River and Liaohe deltas and urban populations in the towns of Songjiang and Tianjin (Fig. 1). Accessions from outside of the area of investigation and of other Phragmites species were included to trace origins and relationships of the Chinese populations. In specific we used our sequence database3,18 and produced new sequences for P. karka and P. japonicus to cover the Asian Phragmites species variation. We used a P. karka population from the Mekong Delta in Viet Nam (N = 40) and samples of P. japonicus (N = 4) from populations in Korea and Sakhalin Island (Russia) (Supplementary Information S1).

DNA extraction

DNA was extracted with the E.Z.N.A. tissue DNA kit (Omega Bio-tek Inc.) from apical leaves conserved in silica gel. The samples were ground in a mortar with quartz sand in liquid nitrogen before adding the E.Z.N.A. extraction buffer, following the protocol for dry plant specimens. The samples were treated with RNAase, eluted with 100 µl elution buffer, and conserved at -20 °C until amplification. DNA quality and quantity were checked in a Nano Drop Spectrophotometer ND-1000 (Saveen Werner) at 280 nm wavelength. DNA concentration exceeded 50 ng µl−1 in all samples.

Chloroplast DNA markers

Two non-coding regions in the chloroplast DNA, the trnT-trnL and the rbcL-psaI, have previously been used to classify Phragmites diversity worldwide3,9,18,19,20,22. Of these two sequences, polymorphism in an indel of variable size in the trnT-trnL region can define most Phragmites haplotypes41. We used the primers of41 to amplify this variable region, and screened the polymorphism of our sample set. 2 µl DNA were added to 18 µl mastermix consisting of 10 µl 2xMastermix (VWR Amplicon), 10 pmol cy-labeled forward primer (5ʹ-CAT TAC AAA TGC GAT GCT CT – 3ʹ)9 , 10 pmol reverse primer (TrnT-R short: 5′- CGT CCG AGC CAT ATC AAA TT- 3ʹ) and sterile water to reach a final volume of 20 µl. Amplification was run in a Peltier Thermal Cycler PTC-200-MJ Research under the following conditions: 94 °C for 3 min, 40 cycles of 94 °C for 30 s, 52 °C for 40 s, 72 °C for 40 s, followed by 72 °C for 7 min. The amplified product (a fragment of about 350 bp) was diluted 20 × and loaded in a 7% acrylamide gel (Reprogel-Long REad) in the fragment analyser ALF Express II DNA Analysis System (Amersham Pharmacia Biotech). 5 µl of diluted PCR product were added to 3 µl loading dye (GE Healthcare) previously mixed with 1 µl each of 100 and 300 bp internal sizers (GE Healthcare). DNA was denaturated in the PTC-200 PCR at 94 °C for 5 min prior to loading in the gel. Electrophoresis conditions were 1500 V, 55 °C, 120 min. The first and the last slots of the gel were loaded with 50–500 bp external sizers (GE Healthcare).

All different trnT-trnTRshort profiles were subsequently sequenced with the trnT-trnL and rbcL-psaI primers developed by9. Annealing temperatures were 50 °C for both primer sets and the PCR conditions were as for the fragments amplification except that the extension time was increased to 1 min. The PCR product was diluted 20 × and sent to Macrogen Korea for Sanger sequencing with forward and reverse primers in an ABI system.

We downloaded all Phragmites trnT-trnL and rbcL-psaI sequences from GenBank22. We aligned our sequences with those downloaded from GenBank with the program BioEdit ver. 7.0 Sequence Alignment Editor40. The initial “clustal” alignment was completed manually. Repeated motifs (minisatellites and microsatellites) and indels were coded as multistate characters. Sequences differing in the number of repeats at microsatellite loci were classified as cp-microsatellite variants following18,22. The final matrix of 134 sequences, 33 of which were from this study, and 1625 base pairs plus 27 multistate characters, was analysed with PAUP ver. 4.0b1042. The parsimony analysis was performed with 37% jackknife deletion, 1000 replicates, jack resampling and stepwise addition43. Individual pairwise genetic distances were also calculated as "total character difference". The matrix of genetic distances was imported in GenAlex ver. 6.444,45 and analysed in a Principal Coordinate Analysis (PCoA). The PCoA obtained is comparable to that of Fig. 2 of18 and includes all haplotypes so far identified in East Asia.

We followed9 to define Phragmites Haplotypes and their intrahaplotypic microsatellites variants18,22 and17,18 for the names of the phylogeographic units.

HKT1 gene polymorphism and gene sequences

We downloaded all Phragmites genomic and mRNA sequences of the PhaHKT1 gene from GenBank20,25. We split the 2000 bp gene sequence into two fragments, HKT1-1 and HKT1-2, of about 1000 bp each. The HKT1-1 fragment contained an indel which was previously found to occur in the Yellow River population (Nanpi and Enchi) but not to have an obvious effect on the ability to tolerate salt, as both mRNA profiles occurred in salt-tolerant reeds25. The HKT1-2 fragment contained two introns, one of which was shown to splice and be translated into proteins of different size in salt-tolerant and salt-sensitive reeds25. We designed a set of primers around the indel in the HKT1-1 part of the gene and another set around the two introns in the HKT1-2 part of the gene, and analysed fragment size polymorphism.

The primers for the fragment containing the indel in the HKT1-1 part of the gene were HKT1 Full-F225 and HKT1-1R short (Supplementary Information S2) and annealing temperature was 55 °C. The primers for the fragment containing the two introns in the HKT1-2 part of the gene were HKT1-2F short and HKT1-2R End (Supplementary Information S2) and annealing temperature was 52 °C. Mastermix and PCR conditions were as described for the trnT-TrnL fragments. The amplified products (between 200 and 300 bp for the HKT1-1 fragments and around 400 bp for the HKT1-2 fragments) were run in a 1.5% agarose gel at 130 V for 2 to 5 h and band profiles scored manually.

We then sequenced all homozygote profiles. The HKT1-1 part of the gene was sequenced with the primers HKT1 Full-F225 and HKT1-1 1000R (Supplementary Information S2) and the annealing temperature was 60 °C. The HKT1-2 part of the gene was sequenced with the primers HTK1-2 1000F and HKT1R End (Supplementary Information S2) and the annealing temperature was 52 °C. Mastermix and PCR conditions were the same as described for the trnT-trnL and rbcL-psaI sequences. The amplified products were sequenced with the primers HKT1 Full-F2, HKT1-2 1000F and HKT1R End.

We calculated the genotypic frequencies of homozygotes and heterozygotes at the two loci as the number of individuals sharing the same profile divided by the total number of individuals in the population. HKT1-1 had two alleles and three genotypic classes, whereas HKT1-2 had three alleles and only two genotypic classes (homozygotes and heterozygotes with three alleles; Table S3). We used a two-loci codominant matrix to calculate allelic frequencies with the program GenAlex.

The three sequences obtained for each homozygote genotype (HKT1 Full-F2, HKT1-2 1000F and HKT1R End) were assembled with the program Genious ver. 6.1.7 (Biomatters) and the resulting consensus sequences of about 2000 bp were aligned with the same program. We compared our sequences, including those of P. karka and P. japonicus, with those of salt-tolerant and salt-sensitive ecotypes deposited in Gene Bank25 also included in the alignment.

SSRs

We studied the allelic variation at loci PAGT4, PAGT8, PAGT9 and PAGT13 previously reported to vary in the Yellow River Delta in relation to soil salinity16 with the primers developed by46. Mastermix, PCR and electrophoresis conditions were the same as described for the trnT-trnL fragment size polymorphism. Annealing temperatures were 50 °C for PAGT 8, PAGT9 and PAGT13 and 54 °C for PAGT 4. The amplified product was diluted 20 × before loading into the fragment analyser.

The chromatograms were aligned with the ALFwin fragment analyser 1.00.36 (Amersham Pharmacia Biotech). Given the polysomic nature of the sample set and the co-occurrence of several ploidy levels in East Asia24, we scored for the presence/absence of alleles of different size. The final binary matrix included 174 taxa, 39 characters and 90 missing data, corresponding to 1.13% of the data, due to fragments of uncertain size.

A PCoA based on the pairwise Euclidean distances of binary haploid profiles was obtained with GenAlex and the resulting structure was investigated in relation to the geographic distribution of the populations, the haplotypic diversity inferred by the TrnT short fragment, and the homo- and heterozygotes at the HKT1 gene loci. The significance of each factor (phylogenetic diversity, HKT1 polymorphism) on the SSR genetic structure was tested with AMOVA (using GenAlex) and with the Bayesian clustering software Structure (ver. 2.3.4)47. Statistical significance for the AMOVA was obtained with 9999 permutations. Concerning Structure, we used the admixture model with correlated allelic frequencies to infer the number of ancestral populations, 300,000 burn-in periods, 1,000,000 MCMC reps after burn-in, and 20 iterations. The initial admixture coefficient alpha was set at 1.047 and K was set from 1 to 10. Structure Harvester48 inferred the most likely K according to the Evanno´s method49.

Genetic diversity indices (samples size, number of alleles, number of effective alleles, Shannon Information Index, gene diversity and unbiased gene diversity) were calculated with GenAlex based on the frequency statistics of binary haploid profiles.

Soil salinity

We used the Soil Quality layer 5 (SQ5) of the Harmonized World Soil Database v 1.250 to extract salinity data at the coordinates of our samples. The SQ5 data file combines soil salinity with soil sodicity and soil phases influencing salt conditions and defines 7 categories of excess salt in the soil. The first five categories range from no to severe excess salt, while categories 6 and 7 include water bodies and permafrost areas. We compared the SQ5 layer with published soil salinity maps of the Yellow River Delta16,51,52,53 and of the Yangtze River Delta54 to ensure updated data and reproducibility of the SQ5 results. As some of our samples in the Yangtze River Delta were collected in areas that were still submersed in the SQ5 layer, we also mined salinity data from our own field work previously conducted at the sampling locations (Xiu Zhen-Li, unpublished data).

We tested differences in excess salt categories among haplotypes and among homo- and heterozygotes at the HKT1-1 and HKT1-2 loci in a multiple factor-ANOVA with Statgraphics Centurion, using General Linear Models and Bonferroni tests at the 95% confidence level. For the Yangtze River Delta populations, we repeated the ANOVA with salinity data measured in the field and also tested Pearson correlations between the SSR allelic frequencies and salinity (as per16), as well as between salinity and HKT1-1 and HKT1-2 allelic frequencies, and homo- and heterozygotes frequencies at the HKT1-1 and HKT1-2 loci, with Statgraphics Centurion.