1 Introduction

Anthropogenic climate warming is driving rapid shifts in species distributions and ecosystem composition globally, fundamentally reorganizing biodiversity (Brodie et al., 2025; Lawlor et al., 2024; Zhang et al., 2024). While species distribution models (SDMs) remain indispensable for forecasting these ecological responses, they overwhelmingly treat species as monolithic, genetically homogeneous entities (Chardon et al., 2020; Qiu et al., 2024). This assumption critically overlooks intraspecific genetic diversity, which is increasingly recognized as a primary determinant of a species adaptive capacity and vulnerability to environmental change (Cheng et al., 2021; Exposito-Alonso et al., 2022; López-Jurado et al., 2019). This knowledge gap is particularly concerning for foundation species, which disproportionately modulate community structure, habitat provision, and biogeochemical cycles (DuBois et al., 2022; Qiao et al., 2021). Consequently, ignoring their intraspecific dynamics risks missing a cryptic but profound mechanism of ecosystem reorganization: climate-driven range shifts and lineage replacement within a species native range.

The common reed, Phragmites australis, provides a powerful model to investigate how intraspecific variation shapes responses to global change (Eller et al., 2017; Meyerson et al., 2016). Phragmites australis is generally considered an allotetraploid, and hexaploids are of allopolyploid origin, whereas octoploids in Asia are most likely autopolyploid (Liu et al., 2022; Wang et al., 2024). As a nearly ubiquitous wetland foundation species, it harbors remarkable genetic and phenotypic diversity, deeply shaped by its complex evolutionary history and multiple polyploidization events (Meyerson et al., 2025; Saltonstall, 2002; Tanaka et al., 2017; Wang et al., 2024). The introduced populations of P. australis from Europe have shown an aggressive expansion capacity in North America (Guo et al., 2013; Saltonstall, 2002), due to smaller genome size (Pyšek et al., 2018) and greater B chromosome expansion (Wang et al., 2024; Wang et al., 2026). Similar lineage-specific dynamics may be unfolding within its native range in China, where multiple ploidy levels coexist across distinct geographic regions.

In China, which is a recognized hotspot for P. australis genetic diversity, three major genetic lineages occupying distinct geographic regions have been identified: the cold-adapted CN lineage (possessing O or M haplotypes) across northern China, the FEAU lineage (possessing the P haplotype, predominantly octoploid) in eastern and southern China, and the SW lineage (possessing U or I haplotypes) restricted to the subtropical southwest (Liu et al., 2022). This correspondence between genetic identity and ploidy level, initially suggested by the geographical overlap of chloroplast haplotypes (An et al., 2012) and chromosome counting-based ploidy level (C. Chen et al., 1993), has since been well established by flow cytometry (Lambertini et al., 2020), allelic number comparisons from microsatellites (Liu et al., 2022), and the alternative allele frequencies of reads mapped to the reference genome (Wang et al., 2024).

The evolutionary legacy of whole-genome duplication, resulting in polyploidy, like in the FEAU lineage, represents a fundamental source of its divergent adaptive potential (Cheng et al., 2021; Kolář et al., 2017). Polyploidy can generate genetic novelty, alter gene expression, and significantly enhance physiological stress tolerance (Bureš et al., 2024; Van de Peer et al., 2017), potentially pre-equipping specific lineages to occupy new geographical range and endure environmental shifts (Cheng et al., 2021; López-Jurado et al., 2019). Specifically, transcriptomic and ecophysiological evidence suggests that the octoploid FEAU lineage possesses large morphological traits (K. Chen et al., 1993; Guo et al., 2025; Liu et al., 2021b, 2026; Yin et al., 2024), stronger salt tolerance and higher thermal tolerance in transcriptome (Wang et al., 2021) than its tetraploid relatives. Over the past decade, FEAU lineage of P. australis has expanded to Caohai Lake, a plateau lake historically not occupied by this species in China (Li et al., 2026; Ran et al., 2025). Previous work suggested asymmetric introgression between the CN and FEAU lineages, but the limited number of microsatellite markers, uneven geographic sampling, and lack of dosage correction for polyploids left its direction, magnitude, and role in climate-driven expansion insufficiently resolved (Liu et al., 2022).

However, a critical question remains: how will the intraspecific differences translate into spatial biogeographic dynamics under rapid climate change? Addressing this question is essential for revealing a subtle yet potent ecological process, climate-driven lineage turnover. This process of intraspecific niche replacement represents a dynamic analogous to biological invasions but operates entirely within the native range (Paudel et al., 2025; Zhao et al., 2024).

Here, we integrate population genomics, multi-site common garden experiments, and species distribution modeling to systematically investigate how historical evolutionary legacies shape climate responses among the major genetic lineages of P. australis in China. We explicitly test the hypothesis that the octoploid FEAU lineage possesses superior thermal tolerance and growth performance, enabling significant range expansion under future warming, whereas the other lineages will exhibit limited capacity to track climatic shifts. By linking genomic identity, physiological tolerance, gene flow, and spatial forecasting, this study provides a mechanistic test of how polyploidy-mediated intraspecific variation governs future biogeographic patterns and the structural integrity of foundation species under global change.

2 Materials and Methods

2.1 Microsatellite genotyping and population genetics

We collected 495 leaf samples of P. australis from 22 provinces across China (Table S1). Approximately 30 mg of dried leaf tissue from each sample was placed in a 2.0-ml tube with a 5-mm glass bead, flash-frozen in liquid nitrogen, and ground using a Tissuelyser II (Qiagen) at 30 Hz. Genomic DNA was isolated using the DNAsecure Plant Kit (Tiangen). DNA integrity was checked on agarose gels, and quality and quantity were measured with a NanoDrop 2000 spectrophotometer (Thermo Scientific). Based on the established nomenclature (Saltonstall, 2002), Phragmites australis haplotypes were identified using sequences from the cpDNA trnT-trnL and rbcL-psaI regions; we determined the haplotypes for 202 out of 495 samples in this and previous studies.

Phragmites australis is an allopolyploid lineage with a base allotetraploid genome, which justifies the common practice in previous microsatellite studies of selecting markers that amplify at most two alleles per individual for most tetraploids (Saltonstall, 2003). In Asia, the prevalent octoploids are most likely autopolyploid derivatives of the tetraploid, allowing ploidy level to be inferred from the number of alleles per locus observed in an individual (Liu et al., 2022).

Hexaploid individuals, on the other hand, appear to be rare and mostly occur in contact zones between different lineages, where they likely originate from inter-lineage hybridization events (Wang et al., 2024). In the present study, the chosen microsatellite markers very rarely produced more than four alleles in a single individual (Table S2); therefore all samples were consistently treated as tetraploid (i.e., four homologous copies). This empirical observation supports the expectation that tetraploids carry two allele copies per locus, while octoploids carry four copies, consistent with a likely autopolyploid origin of the Asian octoploids.

Microsatellite markers were designed with SSRMMD (Gou et al., 2020) based on the transcriptome sequencing of EU and FEAU lineages of P. australis (Wang et al., 2021). The minimum number of repeats was 6 and 5 for trinucleotide and tetranucleotide repeats, respectively. The 64 newly developed and 5 previous (PaGT9, PaGT11, PAGT12, PAGT12, PaGT14) (Saltonstall, 2003) microsatellite markers were tested with single PCR amplification in 16 representative samples. Five newly developed markers failed to get clear bands, and the remained 64 markers were tested in the multiplex PCR. In the multiplex PCR, six marked failed to amplify, and thirteen might have variation in primer regions among samples. Therefore, a total of 58 markers were used in the microsatellite genotyping for all samples (Table S2). All 58 markers were aligned to the reference genome of P. australis (Wang et al., 2024) using minimap2 (Li, 2021), which revealed that five markers mapped to more than one chromosome. During genotyping, eleven markers, including four of the five multi-mapping markers, were removed because more than ten samples exhibited more than four alleles per sample at these loci. An additional five markers were discarded due to amplification failure in more than ten samples. One multi-mapping marker was retained owing to its ideal performance in downstream applications. Consequently, a total of 42 markers (Table S2) were selected for population genetic analysis. These markers are distributed across 18 out of the 25 chromosomes, including one B chromosome.

PCR amplicons (<300 bp) from each sample were pooled and tagged with 8 bp sample-specific barcodes via index primers. The barcoded products were combined into one pool for subsequent library construction. Libraries were prepared following the Illumina standard protocol and subjected to paired-end sequencing (2×150 bp) on an Illumina HiSeq 2500 platform, achieving a mean coverage of >5000× per SSR locus per sample. Raw reads were quality-checked with FastQC. An in-house Perl script named “SSRSeq count” (str_count.pl, available at https://github.com/ccoo22/SSRseq_count) was employed to process high-quality reads and generate an SSR read count table (Cui et al., 2022). Briefly, paired-end reads were merged using FLASH, and the resulting reads were aligned to P. australis sequences where the SSR markers were located via Blastn. The script finally produced the microsatellite read count table summarizing read counts corresponding to alleles with varying repeat numbers for each locus and sample.

We performed microsatellite genotyping following the methodology established by the previous study, using the Perl script str_type.pl and its online implementation SSRSeq V1.1 (Cui et al., 2022). Genotyping was carried out based on read count tables, with allele calling using a threshold set according to ploidy to minimize allelic dropout. The approach incorporated corrections for stutter peaks based on empirical slip ratios and adjusted for amplification bias among alleles with varying repeat numbers. The final genotype outputs containing corrected allele dosages were formatted for subsequent population genetic analyses.

To delineate major evolutionary lineages, we performed principal coordinate analysis (PCoA) based on Bruvo’s genetic distance, which is robust for polyploid data and captures the primary axes of genetic differentiation without imposing a population model. To further investigate admixture and introgression among lineages, we applied Bayesian clustering in STRUCTURE, which assumes Hardy–Weinberg equilibrium within clusters and estimates individual ancestry coefficients.

We performed Principal Coordinate Analysis (PCoA) based on genetic distance using the polysat package (Clark & Jasieniuk, 2011) in R, which is specifically designed for analyzing polyploid microsatellite data. Since allele dosages had already been estimated based on read counts in prior steps, we did not utilize the genotype estimation functions within polysat. Instead, genetic distances between individuals were directly calculated using meandistance.matrix() under the default Bruvo’s distance metric. The resulting distance matrix was then subjected to classical multidimensional scaling via the cmdscale() function.

Meanwhile, we employed a Bayesian clustering approach implemented in STRUCTURE v2.3.4 (Pritchard et al., 2000) to infer individual genetic ancestry. The number of clusters (K) was tested from 1 to 10, with five independent runs per K. Each run consisted of a burn-in period of 1,000,000 steps followed by 1,000,000 MCMC iterations. The optimal K value was determined using Structure Harvester (Earl & vonHoldt, 2012). For the selected K, we aligned the replicate runs with CLUMPP (Jakobsson & Rosenberg, 2007) to obtain a consensus population structure.

Following individual ancestry inference by STRUCTURE, we classified samples into three admixture groups (pure CN, pure FEAU, and mixed) using thresholds of 80% ancestry from a single component. To visualize spatial patterns, individual data were aggregated at the provincial level and displayed on a map of China using proportional pie charts. Regional distributions of CN ancestry were examined through histograms across three geographic regions. We quantified admixture levels using an index derived from ancestry coefficients (1 - 2 × |CN - 0.5|), which was log-transformed for analysis. Linear mixed-effects models with Province as a random effect were used to assess relationships between admixture level and geographic coordinates (latitude and longitude).

A binomial test and a one-sample t-test were used to evaluate whether the proportion and mean of FEAU ancestry in admixed individuals were greater than 0.5, respectively. To infer the direction of introgression between the CN and FEAU lineages, we identified private alleles using reference groups of geographically distant pure individuals. Pure CN here reference individuals were selected from northwestern China (longitude < 100°) with STRUCTURE ancestry coefficient Q_CN > 0.9 (n = 70). Pure FEAU reference individuals were selected from southeastern China (longitude > 115°) with Q_FEAU > 0.9 (n = 44). For each of the 42 microsatellite loci, allele frequencies were calculated separately for the two reference groups. A private allele was defined as an allele with frequency > 0.05 in one reference group and < 0.01 in the opposite group. For every individual, we scored the presence or absence of private alleles from each lineage. The proportion of pure individuals (Q > 0.8 for the respective lineage) carrying private alleles from the opposite lineage was calculated, and binomial tests were used to assess whether these proportions were significantly greater than 0.

2.2 Bioclimatic niche comparisons

We compiled the geographic distribution information of major published P. australis haplotypes in China (O, M, P, U, I). Based on our previous microsatellite and genomic studies as well as conventional nomenclature, we designated haplotypes O and M as the CN lineage (China), P as the FEAU lineage (East Asia and Australia), and U and I as the SW lineage (southwest China, or subtropical lineage). For samples whose haplotypes were not experimentally determined in this study, lineage assignments were inferred based on results from microsatellite-based PCoA. In total, we collected 837 records, comprising 365 CN, 431 FEAU, and 41 SW lineage individuals.

All bioclimatic factors for contemporary climate conditions (1979–2013) were downloaded from WorldClim database (Fick & Hijmans, 2017) for the niche analysis. To address multicollinearity among the bioclimatic variables available, pairwise Pearson correlations were calculated (Figure S1). One variable from each pair exhibiting a correlation > 0.80 was randomly excluded. The initial selection of bioclimatic factors was based on the environmental requirements of P. australis and previous species distribution model study of this species (Guo et al., 2013). The filtering yielded a final subset of eight key predictors for the niche comparison and further MaxEnt modeling: annual mean temperature (bio1), mean diurnal range (bio2), isothermality (bio3), temperature seasonality (bio4), mean temperature of wettest quarter (bio8), precipitation of driest month (bio14), precipitation seasonality (bio15), and precipitation of warmest quarter (bio18). Statistical differences in bioclimatic variables among lineages were assessed using Kruskal–Wallis tests with Benjamini–Hochberg (BH) p-value adjustment.

The environmental space of the CN and FEAU lineages was characterized using the “Characterize Environmental Space” component in the Wallace 2 platform (Kass et al., 2023). This module was employed to visualize and compare the Hutchinsonian niches of the two lineages. Specifically, a Principal Component Analysis (PCA) was performed using the Environmental Ordination module to reduce the dimensionality of the environmental data, followed by a niche overlap analysis conducted via the Niche Overlap module.

2.3 Common garden experiments

To investigate the growth differences and adaptive strategies between the CN and FEAU lineages, we collected growth data, including total biomass, shoot height, density, and specific leaf area (SLA), from two parallel common garden experiments: (1) A previously published common garden experiment (Song et al., 2021) conducted in 2017 across Jinnan (36.43°N, 117.45°E) and Panjin (41.20°N, 122.02°E), using CN (n = 11) and FEAU (n = 9) lineages, for which we determined the haplotype information of all samples; (2) A new common garden experiment established in 2021 across Qingdao (36.36°N, 120.69°E) and Shanghai (30.20°N, 121.29°E), with CN (n = 9) and FEAU (n = 8) lineages (Table S3).

In both experiments, rhizomes were propagated hydroponically and segmented to produce 5–6 segments per population (each with 2–3 buds), which were then transplanted individually into 20 L pots containing a standardized soil substrate. Identical populations were planted simultaneously at both sites under rigorously controlled conditions, with uniform substrate, container size, fertilization, and pesticide application, thereby isolating ambient climate as the sole environmental variable. Garden management was synchronized across sites throughout the growing season, including bi-weekly irrigation (adjusted for ambient precipitation), monthly application of balanced fertilizer from June to September, and standardized pest and disease control.

Measurement protocols were consistent across both gardens. Plant height was defined as the distance from the soil surface to the tallest inflorescence or, in its absence, the highest leaf tip. For each pot, the heights of the three tallest shoots were measured and averaged. Five fully expanded leaves (from nodes 3–5) per pot were randomly selected, placed in labeled plastic bags, and rehydrated in darkness for 12 hours prior to measuring saturated fresh mass. Leaf areas were determined from scanned images using ImageJ software calibrated with a millimeter scale. All leaf samples were then oven-dried at 80 ℃ for 48 hours and weighed. Aboveground biomass was harvested by cutting shoots at soil level, placed in paper bags, and dried at 80 ℃ until constant mass was achieved. Total belowground dry mass was estimated by extrapolating the fresh-to-dry mass ratio obtained from a subsample to the entire root system.

For each trait separately, we calculated the mean and standard error (SE) for each lineage within each garden. To assess the statistical significance of differences between the two lineages within each garden, pairwise Student’s t-tests were performed independently for each garden-trait combination. No correction for multiple testing was applied, as each common garden was considered an independent experimental environment.

2.4 Heat tolerance measurement

Heat tolerance was assessed on a sunny morning after approximately three months of cultivation. For each of the 12 P. australis genotypes (6 FEAU and 6 CN) grown in a common garden (36.38°N, 118.90°E) located in Weifang, Shandong Province, China in 2024, the third or fourth fully expanded leaves were selected. After excision, the leaves were wrapped in moist tissue paper, placed in ziplock bags with ice packs, and transported to the laboratory under dark conditions to minimize transpiration. In the laboratory, the leaves were cut into segments approximately 1 cm in length. These segments were then wrapped in wet filter paper, sealed in plastic bags, and subjected to heat treatment in a water bath. The temperature regimen consisted of sequential exposure to 25 ℃ (ambient control), 38 ℃, 40 ℃, 42 ℃, 44 ℃, 46 ℃, 48 ℃, 50 ℃, 52 ℃, 54 ℃, 57 ℃, and 60 ℃. Following heat treatment, all samples were kept in moist, dark conditions for 24 hours. After dark adaptation, the maximum quantum yield of photosystem II (Fv/Fm) was measured using a pulse-amplitude modulation (PAM) fluorometer (PAM-2500, Walz, Effeltrich, Germany). Each temperature treatment included five biological replicates.

Parameters of heat tolerance (T_crit, T₅₀, and T₉₅) were derived from a logistic decay model that describes the relationship between Fv/Fm and temperature. Parameter estimates were obtained through nonlinear least-squares regression using the nls routine in R. Mean values for each parameter were calculated using bootstrap resampling, implemented within a custom R function (psiiht), which also provided optional diagnostic plotting (Perez et al., 2021). The critical temperature (T_crit) was identified as the temperature at which the slope of the Fv/Fm decline reached 15% of the maximum negative slope observed during heating, with the additional constraint that this point must occur below the temperature causing 50% inhibition. The temperature of 50% inhibition (T₅₀) was defined as the point at which Fv/Fm decreased to half of the value measured at the control temperature. Similarly, the temperature of 95% inhibition (T₉₅) corresponds to a 95% reduction in Fv/Fm relative to the control level. All three parameters were compared between CN and FEAU lineages using a t-test.

2.5 Species distribution models

To assess the impact of climatic changes on the distribution of P. australis lineages, species distribution models were developed using eight selected bioclimatic variables. Potential geographic distributions were projected using the maximum entropy algorithm (MaxEnt) within Wallace (Kass et al., 2023). MaxEnt models species distributions by comparing environmental conditions at occurrence sites with those across the background environment. The model was calibrated with 10,000 randomly generated background points from the study region. Regularization multipliers were tested ranging from 0.5 to 5.0 at intervals of 0.5, and all feature class combinations (Linear, Quadratic, Hinge, Product) were examined. Model selection followed a random k-fold cross-validation approach (k = 4), with the best model chosen based on the lowest corrected Akaike Information Criterion (ΔAICc).

For projecting future distributions across both native and introduced ranges, bioclimatic variables at 5-arcmin resolution were obtained from WorldClim. Future projections for the period 2061–2080 incorporated two representative climate scenarios: the low-emission pathway SSP1-2.6 and the high-emission pathway SSP5-8.5 (IPCC, 2021). Future climate projections were derived from the Centre National de Recherches Météorologiques Earth System Model 2-1 (CNRM-ESM2-1) (Séférian et al., 2019).

3 Results

3.1 Genetic structure and lineage delimitation

Given that STRUCTURE showed optimal K = 2 (Figure S2) and failed to fully resolve the SW lineage (Figure 1A), we used PCoA for lineage delineation (Figure 1B) and STRUCTURE for admixture quantification (Figure 2). This complementary approach allows us to distinguish between evolutionary lineage identity and contemporary gene flow.

Figure 1

Figure 2

The PCoA clearly resolved three major genetic clusters, along with a few introgressive and nuclear-chloroplast discordant samples (Figure 1B). Coordinate 1, explaining 22.05% of the variance, separated the FEAU lineage (haplotype P and related haplotypes such as AS) from the CN lineage (haplotypes O and M). Coordinate 2, which explained 9.36% of the variance, distinguished the SW lineage (haplotypes U and I).

For the Bayesian clustering analysis in STRUCTURE, at K = 2, two genetic clusters were identified, distributed across northern and southern China, respectively (Figure 1A), corresponding to the previously reported CN and FEAU lineage. When K was increased to 3, a southwestern genetic component was further differentiated (Figure S3).

Admixture between CN and FEAU lineages was widespread across China, with a general predominance of FEAU ancestry (Figure 2A). The spatial distribution of admixture groups, however, exhibited distinct regional patterns (Figure 2B–E). Southern China was predominantly occupied by pure FEAU and mixed individuals, northwestern China was dominated by pure CN with a limited presence of mixed types, while northern China was primarily characterized by admixed individuals alongside both pure FEAU and CN lineages (Figure 2B, C). These geographical gradients resulted in a significant longitudinal effect on admixture levels (p = 0.013), whereas no significant latitudinal pattern was observed (p = 0.873) (Figure 2D, E).

Among the 131 individuals classified as admixed (Q_CN and Q_FEAU both between 0.2 and 0.8), 80 individuals (61.1%) had a FEAU ancestry proportion greater than 0.5 (0.569 ± 0.16). Both the proportion (binomial test, p = 0.007) and the mean (one-sample t-test, t = 4.89, df = 130, p < 0.001) were significantly greater than 0.5, indicating a systematic bias toward FEAU ancestry in admixed individuals. Using reference groups of 70 pure CN individuals and 44 pure FEAU individuals, we identified 16 CN-private alleles and 19 FEAU-private alleles across the 42 microsatellite loci. Among 194 pure CN individuals (Q_CN > 0.8), 56 individuals (28.9%) carried at least one FEAU-private allele (binomial test, p < 0.001). Among 168 pure FEAU individuals (Q_FEAU > 0.8), 77 individuals (45.8%) carried at least one CN-private allele (p < 0.001). These results demonstrate bidirectional gene flow between the two lineages, with CN → FEAU introgression being more prevalent at the level of private alleles.

3.2 Differentiation of bioclimatic niche among lineages

Significant differences were observed among lineages across all eight bioclimatic variables analyzed (Figure 1C). For bio1, bio3, bio8, and bio18, the FEAU lineage exhibited values significantly higher than those of CN but lower than those of SW. In contrast, for bio2 and bio4, values in FEAU were significantly lower than in CN but higher than in SW. Meanwhile, CN showed significantly lower values than both FEAU and SW in bio14, but higher values in bio15 (Figure 1D).

Following environmental PCA (Figure S4), niche similarity tests revealed a statistically significant moderate overlap between the CN and FEAU lineages (Figure S5) (Schoener’s D = 0.48, p = 0.01). Specifically, 16% of the environmental space was occupied exclusively by the CN lineage, while 40% was unique to the FEAU lineage. Approximately 60% of the environmental space was shared by both lineages, indicating considerable common habitat suitability despite their distinct ecological characteristics.

3.3 Differentiation of growth performance between two lineages

The four common gardens are all located in eastern China where the CN and FEAU lineages co-occur, with mean annual temperature (Bio1) ranging from 9.3 ℃ (Panjin) to 16.8 ℃ (Shanghai) (Figure 3A). The total biomass of the FEAU lineage was significantly greater than that of the CN lineage in Jinan, Panjin, and Qingdao, but no significant difference was observed in Shanghai (Figure 3B). Similarly, plant height was significantly greater in the FEAU lineage than in the CN lineage in Jinan, but not in the other common gardens (Figure 3C). In contrast, the CN lineage exhibited significantly higher density than the FEAU lineage in Jinan and Panjin, while no significant differences were detected in Qingdao or Shanghai (Figure 3D). For specific leaf area (SLA), no significant differences were found between the two lineages in any of the common gardens (Figure 3E).

Figure 3

3.4 Heat tolerance of two lineages

The heat tolerance parameters T_crit, T₅₀, and T₉₅ were compared between CN and FEAU lineages (Figure 3F; Figure S6; Table S4). The mean T_crit for the CN group was 42.9 ± 1.04 ℃, which was significantly lower than that of the FEAU group (44.2 ± 0.96 ℃; p < 0.05). Similarly, a significant difference was observed in T₅₀, with the CN group exhibiting a mean of 47.7 ± 0.56 ℃ compared to 48.5 ± 0.78 ℃ in the FEAU group (p < 0.05). In contrast, no significant difference was detected in T₉₅ between the two groups; the mean values were 51.6 ± 0.83 ℃ for CN and 52.1 ± 1.36 ℃ for FEAU.

3.5. Potential distribution projection of three lineages

The optimal model configurations and performance metrics varied among the three lineages (Table S5). The CN lineage achieved the best fit with the LQHP feature class and a regularization multiplier of 0.5, resulting in high predictive accuracy (AUCtrain = 0.947) and excellent model stability (CBItrain = 0.973). This model was strongly supported (AICc weight ≈ 1.000). The FEAU lineage performed best under the LQH feature class and RM = 0.5, also showing high training AUC (0.946) and good validation performance (AUCval = 0.898). In contrast, the SW lineage was best modelled using a simpler LQ feature class with RM = 1.0, yielding the highest training AUC (0.965) but lower complexity (7 coefficients). All models demonstrated high predictive performance, with training AUC values exceeding 0.94.

Model projections indicated notable shifts in habitat suitability under future climate scenarios (Figure 4; Figure S7; Table S6). For the CN lineage, the high-emission scenario (SSP5-8.5) resulted in a marked expansion of high-suitability areas to 16.5%, compared to 9.4% under current conditions. In contrast, the FEAU lineage showed the most pronounced response under SSP5-8.5, with high-suitability habitat increasing substantially to 25.2%, while low-suitability areas contracted to 56.7%. The SW lineage remained relatively stable across projections, with the majority of habitats consistently classified as low suitability, though a slight increase in high and medium suitability was observed under SSP5-8.5. These results suggest divergent responses to climate change among lineages, with FEAU exhibiting the greatest potential for range expansion under high-emission scenarios.

Figure 4

4 Discussion

4.1 Ecophysiological basis of thermal tolerance and range expansion

Our findings demonstrate that the octoploid FEAU lineage of P. australis possesses greater heat tolerance and biomass production than the tetraploid CN lineage. Under a high emission scenario (SSP5 8.5), the projected suitable habitat for the FEAU lineage expands by 18.6%, while the CN lineage exhibits a much smaller relative increase. These results support the hypothesis that polyploidy can confer physiological advantages that translate into asymmetric range dynamics within a species’ native range (Bureš et al., 2024; Van de Peer et al., 2017).

The superior heat tolerance of the FEAU lineage, characterized by an approximately 1.3 ℃ higher T_crit and 0.8 ℃ higher T₅₀, aligns with evidence that whole genome duplication can enhance stress tolerance through gene dosage effects, subfunctionalization, or epigenetic modifications (Van de Peer et al., 2017). In P. australis, transcriptomic comparisons have revealed that octoploids upregulate heat shock proteins and maintain photosystem II integrity under elevated temperatures (Wang et al., 2021). Similar patterns occur in other polyploid systems, such as Dianthus broteri (López-Jurado et al., 2019) and Fragaria (Wei et al., 2020), although this advantage can be context dependent (Kolář et al., 2017).

A 1.3 ℃ difference in T_crit could be critical when summer temperatures approach or exceed a lineage’s photosynthetic thermal thresholds. Mean annual air temperature over China has risen by more than 1 ℃, with more frequent and hotter summer days (Wei et al., 2023). In eastern China, extreme heat events are becoming more frequent and intense, with peak summer temperatures increasingly approaching the critical thresholds (42.9 ℃ for CN, 44.2 ℃ for FEAU) identified here (Zhou et al., 2024). Direct evidence linking PSII heat tolerance to survival in wetland plants remains limited (Perez et al., 2021), yet the functional consequence is clear: a lineage that maintains photosynthesis at 44 ℃ rather than ceasing at 43 ℃ will sustain carbon assimilation longer during heatwaves. Unlike woody plants in drylands, where heatwave mortality often results from hydraulic failure (Breshears et al., 2021), P. australis in wetlands is unlikely to face severe water limitation; direct thermal damage to the photosynthetic apparatus is the primary constraint. Even marginal differences in thermal tolerance can be amplified through positive feedback: prolonged carbon assimilation leads to greater biomass, which in turn enhances competitive ability and accelerates range expansion (Cheng et al., 2021). Therefore, the consistent physiological differences documented here are likely to translate into meaningful demographic advantages under the intensifying heat regimes expected across eastern China.

4.2 Niche overlap, admixture, and the trajectory of lineage replacement

Bayesian clustering and principal coordinate analysis resolved three distinct genetic lineages, revealing a broad zone of admixture between the FEAU and CN lineages in northern China. A substantial proportion of individuals exhibited mixed ancestry, with admixture levels increasing significantly with longitude but not latitude, consistent with a previous study using limited microsatellite markers (Liu et al., 2022). This east–west gradient suggests ongoing gene flow, which may be facilitated by both natural and human-mediated dispersal. Eastward water dispersal via the Yellow River (Liu et al., 2021a, 2021b) and historical introductions of P. australis germplasm for reed production in areas such as the Panjin region (Brix et al., 2014) likely contribute to this pattern. However, the persistence of pure CN and pure FEAU individuals within admixed populations implies that reproductive barriers or selection against hybrids exist, potentially driven by ploidy differences (tetraploid versus octoploid) or local adaptation.

The moderate environmental niche overlap between CN and FEAU (Schoener’s D = 0.48) indicates that approximately 60% of their environmental space is shared. Under current climatic conditions, parapatric coexistence and the observed admixture patterns may be maintained by dispersal limitations or competitive trade-offs (Guo et al., 2024; Liu et al., 2026). For example, the CN lineage’s higher shoot density at northern sites may contribute to its persistence in cooler regions, potentially offsetting the demographic disadvantages of its lower heat tolerance and morphological size. However, coexistence theory predicts that when two lineages share a large fraction of environmental space and one possesses a consistent performance advantage, competitive exclusion is likely to occur in a changing environment (Pastore et al., 2021).

Furthermore, the octoploid FEAU lineage likely possesses greater genomic compatibility than the tetraploid CN lineage, leading to asymmetric introgression: hybrids are expected to backcross more frequently with the high-ploidy lineage (Bartolić et al., 2024; Han et al., 2015; Liu et al., 2022; Zohren et al., 2016), gradually shifting the genetic composition of admixed populations toward the FEAU genome. Supporting this hypothesis, private allele analysis revealed bidirectional gene flow, with 28.9% of pure CN individuals carrying FEAU-private alleles and 45.8% of pure FEAU individuals carrying CN-private alleles. Crucially, admixed individuals showed a significant bias toward FEAU ancestry (61.1%), consistent with asymmetric backcrossing favoring the octoploid genome.

Our common garden results corroborate this trajectory: the FEAU lineage produced significantly greater biomass than CN in three gardens (Jinan, Panjin, and Qingdao), but not in the warmest (Shanghai). This exception may reflect other limiting factors in Shanghai or indicate that FEAU’s thermal advantage diminishes near its upper physiological limits. Because the two experiments differed in starting material and growth duration, we focus only on the direction and significance of lineage differences within each site, not on absolute biomass across gardens.

As warming progresses, the competitive balance is expected to shift decisively toward FEAU. The current admixture zone may act as a moving hybrid front, with FEAU gradually displacing CN in warmer regions while potentially maintaining a stable hybrid zone in transitional climates. This projected lineage replacement differs fundamentally from the non-native invasion of P. australis in North America (Saltonstall, 2002), where an introduced European lineage spread into novel territory; in China, the expansion involves only native genotypes. Nevertheless, both processes threaten to erode unique genetic diversity and reduce the species’ adaptive capacity.

4.3 Vulnerability of the SW lineage and potential ecosystem consequences

The SW lineage, with only 41 occurrence records and a restricted distribution, occupies a distinct climatic niche and is projected to remain spatially limited under future scenarios. This stability likely reflects a narrow fundamental niche rather than resilience. Field observations of FEAU encroachment into Caohai Lake, a site historically not occupied by P. australis (Li et al., 2026; Ran et al., 2025), echo patterns seen in other systems where warm adapted lineages displace cold adapted ones (Cheng et al., 2021; Paudel et al., 2025). Although not yet quantified, this trend merits urgent monitoring. The SW lineage should be considered a high priority conservation target, and its habitats must be carefully managed to minimize the inadvertent introduction of FEAU propagules (Tanaka et al., 2017).

The potential replacement of the CN lineage by the FEAU lineage, if realized, could have ecosystem-level consequences beyond genetic composition. As a foundation species, P. australis modulates wetland structure and biogeochemical cycles. The FEAU lineage’s tendency toward greater biomass but lower shoot density may alter light availability, sediment dynamics, and habitat complexity for invertebrates (Yan et al., 2021). However, these ecosystem consequences remain speculative. Direct measurements of decomposition rates, greenhouse gas fluxes, and faunal communities across monospecific stands of each lineage are required to test these hypotheses.

4.4 Limitations, future directions, and conservation implications

While our SDM projections provide a robust framework for assessing vulnerability, they inherently assume unlimited dispersal and an absence of novel biotic interactions (Di Cola et al., 2017; Phillips et al., 2017). The natural expansion of the FEAU lineage will ultimately depend on propagule dispersal distances (seeds and rhizomes) and landscape barriers (Liu et al., 2021a). The current admixture zone in northern China confirms that gene flow is active, but the actual rate of spread may lag model projections. Future work should integrate mechanistic dispersal models and long-term field monitoring of contact zones to validate these spatial predictions. Moreover, precipitation changes and land use conversion were not fully accounted for, both of which could alter habitat suitability across regions (Guo et al., 2013).

Additionally, although polyploidy is strongly implicated in underpinning the FEAU lineage’s thermal advantage, ploidy was not experimentally manipulated in this study. Future common garden experiments utilizing near isogenic lines of different ploidies, or transcriptomic analyses that explicitly partition ploidy effects from lineage effects, are needed to strengthen causal inference (Wei et al., 2020). The potential role of admixture in facilitating the adaptive introgression of heat tolerance alleles also warrants deeper investigation (Suarez-Gonzalez et al., 2018).

The projected replacement of a cold adapted native lineage by a warm adapted native lineage presents a nuanced conservation challenge. Unlike classical biological invasions, no non-native species is involved; nevertheless, the loss of unique regional gene pools (CN and SW) may erode the species’ overall adaptive capacity to future environmental volatility (Exposito-Alonso et al., 2022). To safeguard evolutionary resilience, we recommend three practical measures: using strictly locally sourced germplasm in wetland restoration to prevent the human mediated spread of the FEAU lineage; protecting established hybrid zones as reservoirs of genetic diversity; and establishing long-term monitoring networks in geographic contact zones to detect lineage turnover early (Brodie et al., 2025; Lawlor et al., 2024; Vranken et al., 2021).

Conclusions

The integration of population genomics, common garden experiments, and species distribution modeling reveals three interacting mechanisms that jointly drive lineage turnover under climate warming. First, the octoploid FEAU lineage possesses superior thermal tolerance (1.3 ℃ higher T_crit) compared to the tetraploid CN lineage. Second, the octoploid lineage exhibits intrinsically greater biomass production, a trait that is likely an inherent consequence of polyploidy and further reinforced by its thermal tolerance, together enhancing its competitive ability. Third, asymmetric introgression, resulting from preferential backcrossing of hybrids with the octoploid FEAU lineage, progressively shifts the genomic composition of admixed populations toward FEAU. These mechanisms, coupled with moderate niche overlap and ongoing gene flow, strongly suggest a trajectory of gradual lineage replacement across eastern China. This prediction is supported by initial field observations of FEAU expansion in a plateau lake (Li et al., 2026), but decadal monitoring and dispersal constraints are needed for validation. Incorporating evolutionary history and functional genetic identity into ecological forecasting remains essential as climate change reshapes ecosystems.

Data availability

The raw sequence read data were submitted to the NCBI Sequence Read Archive with the BioProject ID of PRJNA1126092. Data and code are available on Zenodo at https://doi.org/10.5281/zenodo.20362090.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (No. 32470388; U22A20558; 32301317), Shandong Provincial Natural Science Foundation (ZR2024QC197), and Tang Scholar Program. The authors thank Professor Jianquan Liu from Lanzhou University and Dr. Kristin Saltonstall from the Smithsonian Tropical Research Institute for providing the geographic information of the known haplotypes.

Additional information

Author Contributions

L.L. (Lele Liu) led the conceptualization, investigation, methodology, validation, visualization, and writing of the original draft. L.L. also contributed equally to formal analysis and acquired funding. W.S. contributed equally to formal analysis and writing – review & editing, with supporting efforts in investigation. Y.W. contributed equally to conceptualization and funding acquisition, and led the writing – review & editing, with supporting roles in formal analysis and investigation. L.L. (Lele Lin) contributed to conceptualization, funding acquisition, investigation, and writing – review & editing. C.W. contributed to conceptualization, methodology, visualization, and writing – review & editing. H.S. contributed to conceptualization, funding acquisition, investigation, and writing – review & editing. Y.G. contributed equally to writing – review & editing, with supporting roles in investigation, methodology, and visualization. W.G. contributed equally to conceptualization, funding acquisition, and writing – review & editing, and provided overall supervision for the project. All authors reviewed and approved the final manuscript.

Funding

National Natural Science Foundation of China (32470388)

• Lele Liu

MOST | National Natural Science Foundation of China (NSFC) (32301317)

• Huijiaq Song

MOST | National Natural Science Foundation of China (NSFC) (U22A20558)

• Weihua Guo

Shandong Provincial Natural Science Foundation (ZR2024QC197)

• Lele Liu

Additional files

Supplementary file 1.