Field surveys on the QTP were performed from 2010 to 2015, and in total, 254 sites were examined (Supplementary Figure 1). Based on morphological identification, M. spicatum and M. sibiricum were found at 46 and 28 sites, S. pectinata and S. filiformis were found at 40 and 92 sites, and R. trichophyllus and R. subrigidus were found at 35 and 54 sites (Figure 1). In the summer of 2015, our field work was focused on mixed populations and plant material collection for molecular analysis. We did not detect new mixed Myriophyllum populations, and a total of 25 sites (Myriophyllum: 10 sites that have been included in our previous work, Wu et al., 2015; Stuckenia: 8 sites; Ranunculus: 7 sites) contained individuals of both parental species or individuals with morphological ambiguity (Figure 1 and Table 1). The plant materials were also sampled and preserved in silica gel for the molecular identification of hybrids.

The variables climate, soil, and landcover were used as environmental factors for the niche analysis (Table 2). Environmental data for the entire studied region (25–40°N, 75–105°E) were obtained from open sources with a resolution of 1 km (climate: www.worldclim.org; soil: Harmonized World Soil Database, v 1.2; landcover: www.resdc.cn; Table 2). Because autocorrelated variables may bias the niche comparison analyses (Broennimann et al., 2012), we kept one of any set of highly correlated variables (Pearson’s r > 0.75). Thirteen variables were selected as evaluator variables, including five bioclimatic variables (Bio1, Bio2, Bio3, Bio4, and Bio12), altitude, solar radiation, average wind speed in the growing season, three soil parameters (texture type, pH, and salinity), and landcover (Table 2).

The ecological niche of each plant was quantified through both a species distribution model (SDM) and an ordination approach. 1. The most widely used SDM algorithm, maximum entropy (Maxent, Phillips et al., 2006), was chosen to model the potential distribution of each species based on the locations of field-observed occurrence. Through initial runs, Maxent performed better on our dataset than other generally used algorithms, such as generalized additive models, boosted regression trees and random forest, based on the AUC (area under the receiver operator curve) and TSS (maximizing true positive rate/sensitivity and true negative rate/specificity) values. One-fifth of the presence was used as testing data for cross-validation, and the threshold of occurrence was generated based on the sum of the sensitivity and specificity. The algorithm either runs 1,000 iterations or ends at convergence (threshold = 0.00001). 2. An ordination approach, such as principal components analysis (PCA), was used to identify the niches of species observations, with the available ecological data of the entire study area as the background (Broennimann et al., 2012). The observed occurrences were converted into smooth densities and represented using the first two PC axes, which was named the “pca-env” technique by Broennimann et al. (2012). It was proven to perform well in the subsequent niche overlap analysis among different ordination and niche model techniques (Broennimann et al., 2012). When executing the analysis, we used the proportions of sand and clay instead of the categorical variable soil texture.

The niche overlap of each species pair was measured with D metrics (Warren et al., 2008; Broennimann et al., 2012). Statistical tests of niche equivalency were performed for the null hypothesis that the niches of the parental species were identical in the study area. The overlaps of the pseudoreplicate model were estimated by randomly reassigning the presence sites to a pair of species and then compared to the observed niche overlap (Warren et al., 2008). We also performed background/niche similarity tests for the null hypothesis that the empirical overlap between two species is greater than the simulated overlap between the niche of one species and the backgrounds of the other (Warren et al., 2008). Contributor variables of PCA were used for canonical discriminants analysis (CDA), to compare the niche requirements among parent species and their hybrids. Wilks’ test was employed to evaluate the goodness of fit statistic and measure the model performance.

All the above analyses were conducted in R version 3.6.2. The SDM-based analyses were performed using the packages “dismo” (Hijmans et al., 2017) and “ENMTools” (Warren et al., 2010). The Java program Maxent v 3.3.3k was also implemented in “dismo” (Hijmans et al., 2017). The PCA-based analyses were performed in the package “ecospat” (Di Cola et al., 2017). The “candisc” package was used for CDA analyses (Friendly and Fox, 2021).

Our previous work identified hybrids in mixed Myriophyllum populations using microsatellite and cpDNA markers (Wu et al., 2015). Molecular studies have also been suggested for the identification of Stuckenia and Ranunculus hybrids because they are morphologically variable and difficult to distinguish from the parental species (Kaplan, 2008; Zalewska-Gałosz et al., 2015). We conducted a molecular analysis of a total of 303 individuals from all 15 indeterminate populations for the two species pairs (Stuckenia: 91 individuals of 8 populations; Ranunculus: 112 individuals of 7 populations) (Figure 1 and Table 1). We also included parental populations of each species pair in the molecular analysis as reference data (Figure 1). A total of 21 pure Stuckenia populations (75 individuals of S. pectinata and 131 individuals of S. filiformis) and 20 pure Ranunculus populations (110 individuals of R. trichophyllus and 215 individuals of R. subrigidus) were involved.

Total genomic DNA was extracted using the DNAsecure Plant Kit (Tiangen Biotech, Beijing, China). Microsatellite markers (7 loci for Stuckenia, Nies and Reusch, 2004; Wu et al., 2020; 17 loci for Ranunculus, Wu et al., 2017; Wu et al., 2019) and one cpDNA fragment (rpl32-trnL) were chosen to identify kinships. The primers “rpL32-F” and “trnL^(UAG)” (Shaw et al., 2007) were used to amplify and sequence the region. The PCR conditions were performed according to the respective protocols.

Microsatellite genotyping was performed using GeneMarker v1.5 (SoftGenetics, State College, PA, United States). A Bayesian clustering method approach was used for individual genetic assignment, implemented in STRUCTURE v2.3.4 (Pritchard et al., 2000; Falush et al., 2003). Ten independent runs were performed at K = 2 with a burn-in period of 500,000 iterations and 2,000,000 MCMC iterations under the admixture model. The average Q value of all STRUCTURE runs was used to indicate individuals as parental species or hybrids, using Q < 0.10 or Q > 0.90 for pure individuals and 0.10 < Q < 0.90 for hybrids (Burgarella et al., 2009). We also assessed the optimal number of clusters for further genetic divergence (Evanno et al., 2005), when K = 2 was suggested for all species pairs.

To mutually verify hybrid identification and evaluate hybrid class assignment, a model-based Bayesian approach was applied using NewHybrids to further assign all individuals to one of six possible classes through evaluation of the Bayesian posterior probability of membership: two parental species, F1s, F2s, and backcrosses to both parents (Anderson and Thompson, 2002). In NewHybrids program, pure samples of the two species need not be specified for hybrids identification, because the allele frequencies of the two species are not required as a priori. The runs were performed using Jeffreys’ prior distribution, with a burn-in of 200,000 iterations and a run length of 2,000,000 iterations. The threshold of assignment probability was set to 75%. When Myriophyllum populations were analyzed, the same dataset was used as that in Wu et al. (2015). All cpDNA fragment sequences were aligned using the program Mafft v6.7 (Katoh and Toh, 2008). The cpDNA haplotypes were extracted using DNASP v5.10 (Librado and Rozas, 2009). The published rpl32-trnL sequences for each genus (for Ranunculus, R. section Batrachium species were included) were used to confirm that no other closely related species involved. We implemented a median-joining network for the genealogical relationships among the genotypes and haplotypes in NETWORK v10.1 (Bandelt et al., 1999).

The potential distribution of each species was quantified using Maxent. The AUC of the training values ranged from 0.935 to 0.957 for the SDM of the six parental species (Supplementary Table 1), indicating that the models of all species performed well and generated excellent evaluations. The thresholds of presence/absence for each species ranged from 0.211 to 0.276, as suggested by the TSS (Supplementary Table 1). We defined and represented the species potential habitats as “high potential” (>0.6), “moderate potential” (0.4–0.6), “low potential” (threshold-0.4) and “no potential” (<threshold), based on the probability of presence suggested by the models (Supplementary Figure 2).

Temperature, precipitation, and elevation were the most decisive environmental factors for the distribution of M. spicatum on the QTP (Table 2). Its highly suitable region on the QTP was found to include the edge of the plateau and the river valley areas (Supplementary Figure 2A). The distribution of M. sibiricum can be mainly attributed to the conditions of temperature and soil texture. The northern slope of the Himalayas and border areas among Qinghai, Sichuan, and Tibet were the main regions that were suitable for M. sibiricum (Supplementary Figure 2B). Elevation, temperature isothermality, and soil pH were the environmental variables that constrained the distribution of S. pectinata (Table 2). Its suitable areas were mainly at the northeast and southeast edges of the QTP (Supplementary Figure 2C). Elevation, annual temperature, and precipitation mainly influenced the distribution of S. filiformis (Table 2). The species had broad potential distribution areas over the plateau (Supplementary Figure 2D). Elevation, temperature isothermality, and soil texture were the main factors influencing the distribution of R. trichophyllus (Table 2), which was most likely to occur in central Tibet and western Sichuan (Supplementary Figure 2E). Annual temperature, landcover, and solar radiation greatly affected the distribution of S. subrigidus (Table 2). It was the second most common submerged plant on the QTP and had extensive suitable distribution areas (Supplementary Figure 2F).

By multiplying the presence distributions (>threshold) of related species, we obtained the potential areas for the coexistence of each species pair (Figure 2). For Myriophyllum, the two species could be separated by the conditions of temperature, seasonality, elevation, and solar radiation (Supplementary Figure 3). The potential sympatric areas for the two species were river valleys in southern Tibet and western Sichuan (Figure 2A). For Stuckenia, the temperature diurnal range and seasonality, elevation, wind speed, and soil texture were decisive for the coexistence of the two species (Supplementary Figure 4). Their potential sympatric regions were medium-elevation areas in north Qinghai and mosaic areas along the Yarlung Tsangpo River (Figure 2B). The annual temperature and elevation mainly separated the two Ranunculus species on the QTP (Supplementary Figure 5). They might cooccur in central and southern Tibet (Figure 2C).

Obvious niche overlap was suggested in all species pairs (Myriophyllum D = 0.594; Stuckenia D = 0.424; Ranunculus D = 0.545, Table 3). The null hypothesis of niche equivalency was rejected in all species pairs, indicating that the niches of ecologically isolated parental species did not shift to equivalence on the QTP (Table 3 and Supplementary Figure 6). The results of the background test in Myriophyllum and Ranunculus were significant, indicating that the niche overlap between the two parental species pairs was more obvious than that between either parent and random background points (Table 3 and Supplementary Figure 7).

The ecological niche was also measured by an ordination approach. The first two axes in PCA explained 76.61% of the variation in environmental factors (Supplementary Figure 8). Axis 1 was associated with the temperature variables, precipitation, soil texture, and wind speed, while axis 2 was mainly associated with the elevation and solar radiation (Supplementary Figure 8).

Obvious niche overlap was also identified when quantified using the ordination approach (Figure 3). Limited unique niches for each parental species and few changes in niche centers were found for Myriophyllum and Ranunculus (Figures 3A,C). The niche overlap between S. pectinata and S. filiformis was obvious as well, while the high-occurrence density areas of the two species were different (Figure 3B). The values of niche overlap between each species pair were 0.625 (Myriophyllum), 0.446 (Stuckenia), and 0.516 (Ranunculus) (Table 3 and Figure 3).

Similar to the results based on SDM, no niches of species pairs were identical according to the niche equivalency tests (Table 3 and Supplementary Figure 9). The background tests showed that the niche similarity between Myriophyllum species and between Ranunculus species was higher than expected at random on the QTP (Table 3 and Supplementary Figure 10). Consistent with the results of background tests, CDA analyses yielded one discriminant function that explained 100% of the variance in PC contributors, showing ecological niche differences were significant between parent Stuckenia (Wilks’ λ = 0.669, p < 0.001), but not significant between parental Myriophyllum species (Wilks’ λ = 0.705, p = 0.054) or Ranunculus species (Wilks’ λ = 0.774, p = 0.053) (Supplementary Figure 11). While moderate Wilks’ λ values (close to 1) suggested all species pairs occupied some regions that were ecologically similar.

Totals of 82 and 126 multilocus genotypes were revealed in the ambiguous populations of Stuckenia and Ranunculus using microsatellite markers. Under the assumption that K = 2, the two genetic clusters suggested by STRUCTURE analysis corresponded to the parental species (Supplementary Figure 12). All the individuals from parental populations had high posterior probabilities close to 0 and 1, except for DXS of S. pectinata and MQR of R. subrigidus (Supplementary Figure 12). Of the eight ambiguous Stuckenia populations, most contained individuals who were genetically admixed, with a probability range from 0.200 to 0.860 (Table 1, Figure 1B, and Supplementary Figure 12A), while QS was composed of two pure parental species. Of the seven ambiguous Ranunculus populations, all contained hybrids, with a probability range from 0.106 to 0.803 (Table 1, Figure 1C, and Supplementary Figure 12B). Two populations (DXS and MQR) assigned to S. pectinata and R. subrigidus by morphology included hybrid individuals (Table 1, Figure 1C, and Supplementary Figure 12B). The PL of S. filiformis contained one individual of S. pectinata.

The hybrid individuals suggested by NewHybrids were consistent with the results of STRUCTURE with few exceptions (Table 1 and Figure 4). All genotype classes were found in the Myriophyllum hybrids. In MDa, more hybrids were suggested compared to the STRUCTURE results (Table 1 and Figure 4A). Stuckenia hybrids were rarely F1s. The hybrid individuals in GEM suggested by STRUCTURE were assigned to S. filiformis (Table 1 and Figure 4B). We assessed the STRUCTURE results under K = 3 in Stuckenia, which did not suggest an assignment of hybrids to independent cluster, but a split of S. filiformis consistent to genetic structure in Wu et al. (2020). For Ranunculus, most hybrids were F1s, while the hybrid individual in MQR was assigned to backcross-R. subrigidus. Additionally, three individuals with DJ were suggested to be F1s (Table 1 and Figure 4C).

One cpDNA fragment (rpl32-trnL) was sequenced for kinship identification. The lengths of the rpl32-trnL alignment region were 731 bp for Stuckenia and 887 bp for Ranunculus. For Stuckenia, 26 polymorphic sites were found between the parental species (Supplementary Figure 13). Twenty-one haplotypes (H1–H21) were found in S. filiformis populations, and three haplotypes (H22–H24) corresponded to S. pectinata (Supplementary Figure 13A). Among the hybrids suggested by SSR data, only individuals from GD had the haplotype of S. pectinata, whereas cpDNA haplotypes of other hybrids were assigned to S. filiformis (Table 1). For Ranunculus, 21 polymorphic sites were found between the parental species (Supplementary Figure 13B). Two and four haplotypes were found in R. subrigidus (S1–S2) and R. trichophyllus (T1–T4), respectively (Supplementary Figure 13B). The hybrids identified by microsatellite loci from DTR and MQR had haplotypes of R. subrigidus, while other hybrids contained the cpDNA of R. trichophyllus (Table 1).

Overall, combining the microsatellite and cpDNA data revealed the distributions of hybrids. Twenty-three individuals from eight sites were identified as hybrids of Stuckenia (Table 1). The hybrids were suggested to be post-F1s, and S. filiformis tended to be the maternal species (Table 1). Half of the mixed populations contained moderate proportions of hybrid individuals (Table 1). Thirty-one hybrids from nine sites were found in Ranunculus, and they were likewise not dominant in mixed populations except NM (Table 1). Most hybrid individuals of Ranunculus were F1s, while R. trichophyllus was more frequent, and reciprocal hybridization occurred in Qinghai, where R. trichophyllus was rarely distributed (Table 1 and Figure 1).

We assessed the distribution suitability of parents at the locations of their hybrids based on SDM results. We found that most sites of hybrids were suitable for at least one parental species, except that DL of Myriophyllum, HX and MQS of Stuckenia, NM and RT of Ranunculus were rarely suitable for either parent (Supplementary Table 2). According to the PC scores of the hybrid populations in “pca-env” analysis, many were situated in the overlapping niches of parents, while the DL of Myriophyllum, HX, GEM and MQS of Stuckenia, DTR, MQR, RT, and NM of Ranunculus were distributed in the marginal niches of both parents (Supplementary Figure 14 and Supplementary Table 2).