Abstract
The impact of ecological factors on natural hybridization is of widespread interest. Here, we asked whether climate niche influences hybridization between the two closely related plant species Myriophyllum sibiricum and M. spicatum. Eight microsatellite loci and two chloroplast fragments were used to investigate the occurrence of hybridization between these two species in two co‐occurring regions: north‐east China (NEC) and the Qinghai‐Tibetan Plateau (QTP). The climate niches of the species were quantified by principal component analysis with bioclimatic data, and niche comparisons were performed between the two species in each region. Reciprocal hybridization was observed, and M. sibiricum was favoured as the maternal species. Furthermore, hybrids were rare in NEC but common in the QTP. Accordingly, in NEC, the two species were climatically distinct, and hybrids only occurred in the narrow geographical or ecological transition zone, whereas in the QTP, obvious niche overlaps were found for the two species, and hybrids occurred in multiple contact zones. This association between hybridization pattern and climate niche similarity suggests that the level of hybridization was promoted by niche overlap. Compared with the parental species, similar climate niches were found for the hybrid populations in the QTP, indicating that other environmental factors rather than climate were important for hybrid persistence. Our findings highlight the significance of climate niche with respect to hybridization patterns in plants.
Introduction
Natural hybridization is common and considered to be an important driving force for evolution and speciation (Anderson & Stebbins, 1954; Lewontin & Birch, 1966; Barton & Hewitt, 1985; Hewitt, 1988; Arnold, 1997; Ellstrand & Schierenbeck, 2000; Burke & Arnold, 2001). In many cases, hybrids do not persist and reproduce due to genetic incompatibility, ecological selection and competition from parents, live in a tension zone that maintains a balance between dispersal and selection or adapt to distinct habitats (Templeton, 1981; Barton & Hewitt, 1985; Nosil et al., 2009; Abbott et al., 2013). However, the extent and position of hybrid zones are variable (Barton & Hewitt, 1985), and hybridization patterns have been shown to vary among multiple co‐occurring areas of parental species (e.g. Gaskin & Schaal, 2002; Zeng et al., 2011). Ecological factors can also promote or limit interspecific gene flow (Anderson, 1948), and the success and maintenance of hybrids may be more common in novel habitats than in the parental environment, depending on environment‐dependent superiority (Arnold & Hodges, 1995; Arnold, 1997; Barton, 2001; Burke & Arnold, 2001; Seehausen, 2004). Recently, integrating environmental data and genetic analyses has allowed us to seek ecological explanations for evolutionary patterns (Kozak et al., 2008), including gene flow patterns between distinct genetic lineages promoted by niche similarity (Arteaga et al., 2011) or suitability (Ortego et al., 2014).
In this study, environmental data were combined with population genetic data to explore hybridization between two watermilfoil sister species: Myriophyllum sibiricum Komarov and Myriophyllum spicatum L. (Moody & Les, 2010). Both sexual and vegetative reproduction can be observed in the two species, and vegetative propagules include stolons, fragments or turions (Aiken et al., 1979). The two species have different geographical ranges and climatic niche breadths. In particular, M. sibiricum (northern watermilfoil) has an affinity for colder climates and is rarely found south of the mean January isotherm of 0 °C (Aiken & Walz, 1979) with a circumboreal distribution in the northern part of Eurasia and North America (Aiken et al., 1979; Aiken & McNeill, 1980; Faegri, 1982; Ceska & Ceska, 1986; Brown et al., 2014). By contrast, M. spicatum (Eurasian watermilfoil) is widely distributed throughout its native range in Eurasia and northern Africa and its non‐native range of North America (Couch & Nelson, 1985), although it is unable to survive in very cold regions and has only been recorded south of 60°N in Europe (Faegri, 1982) and south of 55°N in Asia (Ceska & Ceska, 1986). However, there are overlapping areas, and the ecogeographic isolation between the two sister species is incomplete. Frequent hybridization between the species has been reported in North America (Moody & Les, 2002, 2007; Sturtevant et al., 2009), although similar investigations have not been conducted in their native ranges in Eurasia.
In China, M. sibiricum has a disjunctive distribution between north‐east China (NEC) and the Qinghai‐Tibetan Plateau (QTP), whereas M. spicatum is distributed throughout the country (Yu et al., 2002). Hybridization between these two species likely occurs in their overlapping areas in NEC and the QTP. However, these two areas have distinct climates due to differences in latitude and altitude. Therefore, under the influence of distinct climates, the climate niches of the two species may vary between NEC and the QTP, which in turn could influence the frequency and patterns of hybridization. To test these hypotheses, we used molecular markers and bioclimatic variables (i) to explore the occurrence and pattern of hybridization between the two watermilfoil species in NEC and the QTP and (ii) to compare the climate niches of the two species in each region to test whether climate niche differences influence hybridization patterns.
Materials and methods
Plant materials
A total of 474 individuals (ramets) from 50 sites were collected, of which 330 were from 24 sites in the QTP and 144 were from 26 sites in NEC (Fig. 1; see Appendix S1). The populations were morphologically identified in the field as parental species or hybrids by counting leaf segment numbers in the majority of collected individuals, because there are different ranges of leaf segment numbers in M. sibiricum (6–18), in M. spicatum (24–36) and in hybrids (16–28) (Moody & Les, 2002, 2007). One site, MDa, which contained both species, was coded as MDab and MDap for M. sibiricum and M. spicatum, respectively. Plants were collected arbitrarily at intervals of at least 10 m to avoid sampling different ramets from the same genet. The fresh leaves were dried in silica gel in the field and frozen at −20 °C after being transported to the laboratory. Voucher specimens were deposited in the herbarium of Wuhan University (WH).
Geographic distribution of the 51 Myriophyllum populations from the (a) QTP and (b) NEC regions. Circles are coloured based on morphology: filled black circles, M. sibiricum; open circles, M. spicatum; filled grey circles, hybrids; and half‐open and half‐filled black circle, two species co‐occur. The bar plot on the left side depicts the structure admixture coefficients when K = 2. A single horizontal bar displays the membership coefficient of each individual, with sample site labels shown on the left. The M. spicatum cluster is shown in white, and the M. sibiricum cluster is shown in grey.
Genetic analyses
Total genomic DNA was extracted using the DNA secure Plant Kit (Tiangen Biotech, Beijing, China). Half‐samples were genotyped at twenty microsatellite loci isolated for M. spicatum (Wu et al., 2013) in preliminary experiments, and then eight loci with species‐specific alleles were chosen (Myrsp2, Myrsp4–7, Myrsp12 and Myrsp15–16) for subsequent analysis. The reaction conditions and the genotyping processes were according to Wu et al. (2013). Both M. sibiricum and M. spicatum are hexaploid (Löve, 1961; Aiken et al., 1979) and their allele copies are ambiguous; therefore, identifying their exact genotypes based on allele sizes is difficult. Several methods have been used to address this problem (e.g. Bruvo et al., 2004; Esselink et al., 2004). In this study, microsatellites were treated as dominant markers, and the microsatellite data were converted into a binary matrix based on the presence/absence (1/0) of observed alleles (Lynch, 1990; Samadi et al., 1999), according to the methods of Andreakis et al. (2009). The transformation was performed using the r package polysat (Clark & Jasieniuk, 2011), and the exported binary matrix was used for the following analyses. Clone identification was performed using genotype and genodive (Meirmans & Van Tienderen, 2004) with the criterion of treating individuals with the same multilocus genotype as a clone (genet).
A Bayesian clustering method implemented in the program structure version 2.3.4 (Pritchard et al., 2000) was used to identify the lineage proportion of each individual and to detect interspecific hybridization or introgression between M. sibiricum and M. spicatum. The structure analyses were run on the QTP and NEC separately. The current version of structure was adapted to use dominant markers (Falush et al., 2007). The focus was on hybridization between the two watermilfoil species, and therefore, the number of clusters was set to K = 2. The format of the input file was sorted according to the document. A burn‐in period of 20 000 iterations and 100 000 Markov Chain Monte Carlo (MCMC) iterations under the admixture model was set without prior population information (e.g. characteristics of sampled individuals or geographical sampling locations), and 10 replicate runs were performed at K = 2. The parameter Q was used to denote the probability of the admixture proportions for each individual, assuming purebred individuals should have Q‐values near 0 or 1, whereas hybrids should have intermediate values. Another approach, principal coordinate analysis (PCoA), was also used to examine and verify the genetic clusters of Myriophyllum samples, which was conducted using the GenAlEx 6.5 software (Peakall & Smouse, 2012). We used polysat to calculate the pairwise genetic distances among individuals, which was measured by the ‘Lynch distance’ (Lynch, 1990). Only individuals with different genotypes were considered in this analysis.
Two chloroplast fragments were used to identify maternal kinship. The primers ‘c’ and ‘f’ reported by Taberlet et al. (1991) and ‘rpL32‐F’ and ‘trnL(UAG)’ of Shaw et al. (2007) were used to amplify and sequence the trnL‐F and the rpl32‐trnL regions, respectively. PCR was performed using 10–30 ng of genomic DNA, 0.1 μm of each primer, 0.2 mm of each dNTP, 2 mm of MgCl2 and 0.6 U of ExTaq DNA polymerase (TaKaRa, Otsu, Japan) in a volume of 25 μL under the following conditions: 3 min at 95 °C, followed by 35 cycles of 30 s at 95 °C, 30 s at 55 °C and 90 s at 72 °C, and then a final 5 min extension at 72 °C. Amplifications were conducted in a Veriti 96‐well thermal cycler (Applied Biosystems, Foster City, CA, USA). The PCR products were purified and sequenced in both directions by the Beijing Genomic Institute in Wuhan, China. Only individuals with different genotypes were amplified and sequenced. For populations that contained no hybrids, a maximum of six samples per population were sequenced. For populations that contained hybrids, all samples with distinct genotypes were sequenced. The two chloroplast fragments were combined for subsequent analyses, and the sequences were aligned using the program mafft 6.7 (Katoh & Toh, 2008). Identical sequences were collapsed into a single haplotype using dnasp 5.10 (Librado & Rozas, 2009), and all sequences of different haplotypes were deposited in GenBank (Accession Nos. KP031615–KP031642). To interpret the genealogical relationships among the haplotypes, a median‐joining network (Bandelt et al., 1999) was generated from the haplotype sequence data using network 4.5.1.6 (http://www.fluxus-engineering.com). In the network analysis, gaps with two or more base pairs were coded as single mutation events (Simmons & Ochoterena, 2000).
Environmental data and niche analysis
We obtained information on environmental variables across the entire study areas (the QTP: 22–43°N and 70–105°E; NEC: 40–65°N and 112–135°E) with a resolution of 10′ from 19 Bioclim layers (Bio1‐Bio19; http://worldclim.org/bioclim; Busby, 1991; Hijmans et al., 2005). A correlation test was conducted using the pasw statistics 18 software (SPSS Inc., Chicago, IL, USA) to remove highly correlated variables, which may have biased subsequent analyses (Graham et al., 2004; Walker et al., 2009). If two variables had a correlation coefficient >0.75, we considered them highly correlated and kept the variable that resulted in a greater amount of explained inertia. The remaining environmental variables included BIO2 = mean diurnal range; BIO3 = isothermality; BIO4 = temperature seasonality; BIO11 = mean temperature of coldest quarter; BIO12 = annual precipitation; BIO14 = precipitation of driest month; and BIO15 = precipitation seasonality.
The environmental space, representing the realized niche, was quantified by the first two components of a principal component analysis (PCA) based on the available environmental data (seven environmental variables) of the entire study areas (i.e. background), referred to as PCA‐env in Broennimann et al. (2012). Among the different ordination or environmental niche modelling techniques, the PCA‐env performed best with respect to assessing niche overlap. Population occurrences were converted into smooth densities using a kernel function and plotted in the gridded environmental space. Niche overlap was quantified with the D metric (Schoener, 1970; Warren et al., 2008; Broennimann et al., 2012), which varied between 0 (no overlap) and 1 (total overlap). Only genetically pure populations (Q‐value >0.9 or <0.1) were included in the niche overlap analysis, which were run on the two regions separately. Statistical tests of niche similarity between the two Myriophyllum species, which were adapted to environmental space according to Broennimann et al. (2012), were run on the two regions separately. The observed niche overlap was compared against a null distribution of 200 simulated overlap values. The null hypothesis was rejected (P < 0.05) when the value of observed overlap was greater than the simulated one, which indicated that the niches of the species were more similar than would be expected at random (Broennimann et al., 2012). Both the two niche analyses can identify the extent of niche similarity between parental species and help us understand the influence of niche similarity on hybridization patterns. We also performed niche similarity test between the parental species and their hybrids to explore whether the fitness of hybrids was associated with climate. The analysis was only conducted in the QTP, as hybrids were only found at two sites (see results). All of the above analyses were conducted in R version 3.1.1 (R Development Core Team, 2014) with the ‘ecospat’ package (Broenniman et al., 2014). The framework was built by Broennimann et al. (2012) and used in Petitpierre et al. (2012).
Results
Hybridization in the QTP and NEC
Of the 26 sites in NEC, nine populations were defined as M. spicatum, 15 as M. sibiricum, and two as hybrid based on leaf segment numbers. Of the 24 sites in the QTP, eight populations were defined as M. spicatum, nine as M. sibiricum, and eight as hybrid, as both parental species were found at the MDa site (Fig. 1).
A total of 148 amplification variants (alleles) were revealed at eight microsatellite loci, and the number of genotypes for each locus ranged from 12 to 158 (see Appendix S2). A total of 264 multilocus genotypes were found among the 474 individuals, and the genotype number of each population ranged from one to 17 (see Appendix S1). Under the assumption of K = 2, structure analysis revealed two genetic clusters corresponding to M. spicatum and M. sibiricum, which were approximately consistent with the groupings based on morphological traits with a few exceptions. Of the 17 M. spicatum populations, all individuals were assigned to the M. spicatum genetic cluster with high posterior probability (Fig. 1). Of the 24 M. sibiricum populations, most individuals were assigned to the M. sibiricum genetic cluster with high posterior probability; however, two individuals of population ARb were assigned to the M. spicatum genetic cluster with a probability higher than 0.95, and one of ARb and two of MDab were admixed with a probability range from 0.569 to 0.694 (Fig. 1 and Table 1). Of the 10 hybrid populations, all individuals of populations DJa, DJb and DRa, and two of FJ were assigned to the M. spicatum genetic cluster with probabilities higher than 0.95, one of MK was assigned to the M. sibiricum genetic cluster with a probability higher than 0.95, and the rest of the individuals were admixed with a probability range from 0.520 to 0.863 (Fig. 1 and Table 1).
Number of hybrids, Q‐value of structure results revealed by microsatellite markers, cpDNA haplotypes and number of leaf segments for the 12 populations involved in hybridization events
| Populations | Number of hybrid ramets/genets (sample size) | Q‐value | Haplotype | Leaf segments |
| SG | 15/3 (15) | 0.310–0.360 | A3 | 14–28 |
| ARb | 1/1 (20) | 0.431 | B4 | 14–16 |
| DRa | 19/7 (19) | 0.991–0.993 | B4 B5 | 22–34 |
| DJa | 15/8 (15) | 0.959–0.996 | B5 | 18–24 |
| DJb | 12/1 (12) | 0.993–0.994 | B5 | 18–22 |
| DX | 20/9 (20) | 0.539–0.802 | B1 | 18–32 |
| MK | 7/6 (7) | 0.158–0.388 | B8 B9 | 16–26 |
| MDa | 2/2 (30) | 0.306–0.346 | A1 | 18–26 |
| DL | 13/3 (13) | 0.303–0.333 | B1 B6 | 16–28 |
| KD | 13/8 (13) | 0.137–0.406 | B6 | 12–26 |
| LD | 3/2 (3) | 0.282–0.306 | A3 | 20–26 |
| FJ | 5/4 (7) | 0.373–0.480 | B6 | 20–26 |
| Populations | Number of hybrid ramets/genets (sample size) | Q‐value | Haplotype | Leaf segments |
| SG | 15/3 (15) | 0.310–0.360 | A3 | 14–28 |
| ARb | 1/1 (20) | 0.431 | B4 | 14–16 |
| DRa | 19/7 (19) | 0.991–0.993 | B4 B5 | 22–34 |
| DJa | 15/8 (15) | 0.959–0.996 | B5 | 18–24 |
| DJb | 12/1 (12) | 0.993–0.994 | B5 | 18–22 |
| DX | 20/9 (20) | 0.539–0.802 | B1 | 18–32 |
| MK | 7/6 (7) | 0.158–0.388 | B8 B9 | 16–26 |
| MDa | 2/2 (30) | 0.306–0.346 | A1 | 18–26 |
| DL | 13/3 (13) | 0.303–0.333 | B1 B6 | 16–28 |
| KD | 13/8 (13) | 0.137–0.406 | B6 | 12–26 |
| LD | 3/2 (3) | 0.282–0.306 | A3 | 20–26 |
| FJ | 5/4 (7) | 0.373–0.480 | B6 | 20–26 |
Q‐values close to 0 or 1 represent purebred M. sibiricum or M. spicatum, respectively.
Number of hybrids, Q‐value of structure results revealed by microsatellite markers, cpDNA haplotypes and number of leaf segments for the 12 populations involved in hybridization events
| Populations | Number of hybrid ramets/genets (sample size) | Q‐value | Haplotype | Leaf segments |
| SG | 15/3 (15) | 0.310–0.360 | A3 | 14–28 |
| ARb | 1/1 (20) | 0.431 | B4 | 14–16 |
| DRa | 19/7 (19) | 0.991–0.993 | B4 B5 | 22–34 |
| DJa | 15/8 (15) | 0.959–0.996 | B5 | 18–24 |
| DJb | 12/1 (12) | 0.993–0.994 | B5 | 18–22 |
| DX | 20/9 (20) | 0.539–0.802 | B1 | 18–32 |
| MK | 7/6 (7) | 0.158–0.388 | B8 B9 | 16–26 |
| MDa | 2/2 (30) | 0.306–0.346 | A1 | 18–26 |
| DL | 13/3 (13) | 0.303–0.333 | B1 B6 | 16–28 |
| KD | 13/8 (13) | 0.137–0.406 | B6 | 12–26 |
| LD | 3/2 (3) | 0.282–0.306 | A3 | 20–26 |
| FJ | 5/4 (7) | 0.373–0.480 | B6 | 20–26 |
| Populations | Number of hybrid ramets/genets (sample size) | Q‐value | Haplotype | Leaf segments |
| SG | 15/3 (15) | 0.310–0.360 | A3 | 14–28 |
| ARb | 1/1 (20) | 0.431 | B4 | 14–16 |
| DRa | 19/7 (19) | 0.991–0.993 | B4 B5 | 22–34 |
| DJa | 15/8 (15) | 0.959–0.996 | B5 | 18–24 |
| DJb | 12/1 (12) | 0.993–0.994 | B5 | 18–22 |
| DX | 20/9 (20) | 0.539–0.802 | B1 | 18–32 |
| MK | 7/6 (7) | 0.158–0.388 | B8 B9 | 16–26 |
| MDa | 2/2 (30) | 0.306–0.346 | A1 | 18–26 |
| DL | 13/3 (13) | 0.303–0.333 | B1 B6 | 16–28 |
| KD | 13/8 (13) | 0.137–0.406 | B6 | 12–26 |
| LD | 3/2 (3) | 0.282–0.306 | A3 | 20–26 |
| FJ | 5/4 (7) | 0.373–0.480 | B6 | 20–26 |
Q‐values close to 0 or 1 represent purebred M. sibiricum or M. spicatum, respectively.
Similar results were obtained with the PCoA analysis. The genets assigned to the M. spicatum cluster based on the structure analysis were gathered into one group, including all genotypes of the DJa, DJb and DRa populations, two of the FJ and two of the ARb populations. The genets assigned to the M. sibiricum cluster were gathered into another group that included one genotype of population MK. These two groups were apparently separated, and the admixed individuals identified in the structure analysis were located between them, including one individual of ARb and two of MDab (Fig. 2).
PCoA performed with pairwise genetic distances (Lynch distance) among all genotypes (genets) from all populations calculated by polysat. For abbreviations of hybrid populations, see Fig. 1 and Appendix S1.
A total of 225 individuals were sequenced in the trnL‐F and rpl32‐trnL regions. The length of trnL‐F and rpl32‐trnL ranged from 903 to 963 bp and from 1274 to 1307 bp, respectively. The aligned length of the combined sequences was 2199 bp with 27 polymorphic sites, including 15 substitutions and 12 indels (see Appendix S3). The sequences were collapsed into 14 haplotypes (A1–A5 and B1–B9). The haplotype network consisted of two distinct clades. Clade A included 5 haplotypes (A1–A5) corresponding to M. spicatum, and clade B included nine haplotypes (B1–B9) corresponding to M. sibiricum (Fig. 3). Among the individuals assigned to M. spicatum based on microsatellite data, most individuals had haplotypes from clade A, whereas all individuals of populations DJa, DJb and DRa had haplotypes from clade B, which indicated that chloroplast capture due to repeated backcrossing to M. spicatum had occurred (Fig. 3 and Table 1; see Appendix S1). All individuals assigned to M. sibiricum had haplotypes from clade B, including one individual of population MK (Fig. 3 and Table 1; see Appendix S1). Among the admixed individuals identified by microsatellite loci, individuals of populations SG and LD and two individuals of population MDab had haplotypes from clade A, whereas individuals of populations DL, DX, FJ, KD and MK, and one individual of ARb had haplotypes of clade B (Fig. 3 and Table 1).
The network of 14 haplotypes based on two chloroplast fragments. Ellipses represent different haplotypes (A1–A5 and B1–B9) with population names indicated inside. Black dots represent inferred interior nodes that were absent in the samples. For abbreviations of populations, see Fig. 1 and Appendix S1.
Overall, combining the microsatellite and cpDNA data revealed the distributions of M. sibiricum, M. spicatum and hybrids. Of the 26 sites in NEC, 15 sites only contained M. sibiricum, nine sites only contained M. spicatum, one site (LD) only contained hybrids, and one site (FJ) contained M. spicatum and hybrids. Of the 24 sites in the QTP, eight sites only contained M. spicatum, six sites only contained M. sibiricum, seven sites only contained hybrids, two sites (ARb and MDa) contained both parents and hybrids, and one site (MK) contained M. sibiricum and hybrids. Among all genets involved in hybridization events, approximately 90% (47 of 54 genets) contained M. sibiricum lineage haplotypes (Table 1).
Niche overlap and similarity
The percentage of inertia explained by the two PCA axes was 81% (Table 2). The first axis was primarily associated with the temperature variables, whereas the second axis was primarily associated with precipitation (Table 2). Visualization of the environmental spaces of M. sibiricum and M. spicatum is represented in Fig. 4. When quantified, the values of niche overlap (D) between M. spicatum and M. sibiricum were 0.387 and 0.103 in the QTP and NEC, respectively (Fig. 5a,b). Similarly, the niche similarity hypothesis was rejected in the QTP (P = 0.030) but not in NEC (P = 0.323) (Fig. 5a,b), indicating that the climate niches of the two species were more similar than would be expected at random in the QTP, but not in NEC. For the hybrids in the QTP, tests of niche similarity to both parental species were rejected (P = 0.045 compared to M. sibiricum and P = 0.005 compared to M. spicatum, Fig. 5c,d), indicating no distinct climate niche in the hybrids.
Contributions of the first two axes to environmental space
| Bioclimatic variable | PC1 (48.38%) | PC2 (32.59%) |
| BIO 2 (Mean diurnal range) | 0.1193946 | −0.57321328 |
| BIO 3 (Isothermality) | 0.5136370 | 0.03935316 |
| BIO 4 (Temperature seasonality) | −0.5028157 | −0.14567173 |
| BIO 11 (Mean temperature of coldest quarter) | 0.4932938 | 0.15532537 |
| BIO 12 (Annual precipitation) | 0.2358088 | 0.49516044 |
| BIO 14 (Precipitation of driest month) | −0.2052431 | 0.53638682 |
| BIO 15 (Precipitation seasonality) | 0.3578115 | −0.30271597 |
| Bioclimatic variable | PC1 (48.38%) | PC2 (32.59%) |
| BIO 2 (Mean diurnal range) | 0.1193946 | −0.57321328 |
| BIO 3 (Isothermality) | 0.5136370 | 0.03935316 |
| BIO 4 (Temperature seasonality) | −0.5028157 | −0.14567173 |
| BIO 11 (Mean temperature of coldest quarter) | 0.4932938 | 0.15532537 |
| BIO 12 (Annual precipitation) | 0.2358088 | 0.49516044 |
| BIO 14 (Precipitation of driest month) | −0.2052431 | 0.53638682 |
| BIO 15 (Precipitation seasonality) | 0.3578115 | −0.30271597 |
Contributions of the first two axes to environmental space
| Bioclimatic variable | PC1 (48.38%) | PC2 (32.59%) |
| BIO 2 (Mean diurnal range) | 0.1193946 | −0.57321328 |
| BIO 3 (Isothermality) | 0.5136370 | 0.03935316 |
| BIO 4 (Temperature seasonality) | −0.5028157 | −0.14567173 |
| BIO 11 (Mean temperature of coldest quarter) | 0.4932938 | 0.15532537 |
| BIO 12 (Annual precipitation) | 0.2358088 | 0.49516044 |
| BIO 14 (Precipitation of driest month) | −0.2052431 | 0.53638682 |
| BIO 15 (Precipitation seasonality) | 0.3578115 | −0.30271597 |
| Bioclimatic variable | PC1 (48.38%) | PC2 (32.59%) |
| BIO 2 (Mean diurnal range) | 0.1193946 | −0.57321328 |
| BIO 3 (Isothermality) | 0.5136370 | 0.03935316 |
| BIO 4 (Temperature seasonality) | −0.5028157 | −0.14567173 |
| BIO 11 (Mean temperature of coldest quarter) | 0.4932938 | 0.15532537 |
| BIO 12 (Annual precipitation) | 0.2358088 | 0.49516044 |
| BIO 14 (Precipitation of driest month) | −0.2052431 | 0.53638682 |
| BIO 15 (Precipitation seasonality) | 0.3578115 | −0.30271597 |
Climatic niche dynamics between M. sibiricum and M. spicatum in the (a) QTP and (b) NEC regions. The solid and dashed contour lines indicate 100% and 50%, respectively, of the available environment as backgrounds. Green, red and blue areas represent the unique niche of M. sibiricum, the unique niche of M. spicatum and the shared niche of the two species, respectively. The red solid arrows represent the change in the centre of the niche between the two species.
Results of niche similarity tests for pairwise comparisons between M. sibiricum and M. spicatum in NEC and among the two species and their hybrids in the QTP. Bars with a diamond represent the observed niche overlap, and histograms represent simulated niche overlaps (grey bars). The tests were repeated for 200 iterations, and the null hypothesis of niche similarity was rejected if the observed D fell outside the 95% probability threshold of the stimulated values. The D and P values of the tests are shown.
Discussion
Our study revealed reciprocal hybridization between M. spicatum and M. sibiricum in their indigenous range in China as was reported in the invasive range of M. spicatum in North America (Moody & Les, 2002). However, the hybridization was asymmetrical because the situation with M. spicatum as the female parent was only found in three sites of the 12 sites with hybrids, which was further supported by the chloroplast capture that resulted from repeated backcrossing to M. spicatum as the male parent (e.g. Jackson et al., 1999; Tsitrone et al., 2003; Fehrer et al., 2007; Acosta & Premoli, 2010), which was not mentioned in previous studies (e.g. Moody & Les, 2002, 2007; LaRue et al., 2013a). The asymmetry was also observed in the frequency of hybridization between the QTP and NEC regions. Among the 24 sites of the QTP, hybrids were found at 10 sites and repeated backcrossing occurred, which indicated a hybrid swarm in the QTP region. By contrast, among the 26 sites of NEC, hybrids were only found at two sites.
In Eurasia, with the exception of the QTP, M. sibiricum is distributed more northerly than M. spicatum with only a narrow range of overlap (Faegri, 1982; Ceska & Ceska, 1986), and this geographical pattern is also present in the NEC region (Fig. 1b). Accordingly, a small amount of overlap and distinct climate niches were observed between the two species in this region (Figs 4b and 5b). Only two hybrid populations (LD and FJ) were found in the narrow geographical or ecological transition region, and this pattern is commonly observed when competition pressures from parental species and environmental selection are strong (Barton & Hewitt, 1985; Arnold, 1997; Barton, 2001). In the QTP, hybrids were found 10 of 24 studied sites, indicating a high level of hybridization. The mosaic distribution of the two parental species (Fig. 1a) and the obvious overlap in climate niches (e.g. the niche of M. sibiricum was almost entirely embedded into that of M. spicatum) (Figs 4a and 5a) suggests that there is more opportunity for hybridization in the QTP. The hybrid populations were scattered throughout different regions of the QTP (Fig. 1a), and populations from different regions had different cpDNA haplotypes (Table 2), suggesting that multiple contact zones were present and that hybridization occurred independently in these zones. The extensive hybridization was closely associated with the mosaic distribution and niche overlap of the two parental species in the QTP. As only genetically pure populations were included in the niche overlap analysis, hybridization had little effect on genetic changes related to the climate niche changes in the two Myriophyllum species, which differs from some examples of niche shift resulting from genetic changes following hybridization during biological invasion (e.g. Mukherjee et al., 2012; Thornton & Murray, 2014). Therefore, for the two Myriophyllum species, climate niche overlap was a primary reason for the extensive hybridization observed in the QTP.
Among the 10 sites that contained hybrids in the QTP, eight were dominated by hybrids, indicating that the hybrids were fitter than the parental species at these sites (Arnold & Hodges, 1995). Niche separation between parental and hybrid lineages is often necessary for hybrids to establish and persist (Lewontin & Birch, 1966; Templeton, 1981); for example, extensive hybridization between Engelmann and scrub oaks was associated with the distinct climate niche of their hybrids (Ortego et al., 2014). However, for Myriophyllum in the QTP, the climate niche of the hybrids was more similar to that of either parental species than would be expected by random chance (Fig. 5c,d), suggesting that environmental factors other than climate contributed to the high fitness of the hybrids. Hybrid fitness is often related to transgressive characters (Schwarzbach et al., 2001). Two transgressive traits, faster growth and higher resistance to herbicides, were observed in hybrids between M. sibiricum and M. spicatum in North America (LaRue et al., 2013b). However, these same traits may not be present in hybrids in the QTP due to environmental differences between the two continents. More detailed investigations will be needed to explore the transgressive phenotypes and hybrid fitness in the QTP.
In the invasive range of M. spicatum in North America, extensive hybridization between M. spicatum and M. sibiricum was revealed by recent molecular studies (Moody & Les, 2007; LaRue et al., 2013a). The phenomenon was in accordance with the large overlapping areas of the distributions of these two species in North America (Brown et al., 2014). In China, with the exception of the QTP, the environmental space of M. sibiricum for January mean temperature was below approximately −19 °C (estimated by the Bioclim data of our sampling sites), whereas in North America, M. sibiricum is abundant north of the mean January isotherm of 0 °C (Aiken & Walz, 1979; Aiken, 1981). The M. sibiricum dispersal route was inferred to travel from Eurasia to North America, and the splitting time was estimated to be ca. 0.59 MYA (Chen et al., 2014). Although the climate niche was not completely quantified, the M. sibiricum populations on the two continents potentially occupy different climatic niches, perhaps due to long‐term ecological and evolutionary changes. Therefore, during the subsequent spread following its introduction into North America in the 1940s (Couch & Nelson, 1985), M. spicatum hybridized with M. sibiricum frequently, probably due to their partly overlapping climatic niches.
Taken together, our data reveal that the hybridization pattern of two closely related species, M. spicatum and M. sibiricum, in two co‐occurring regions was associated with overlap and similarity between their climate niches. Our study provides empirical evidence for climate niche overlap promoting hybridization in plants. The similarity between the climate niche of the hybrids and those of the two parental species suggest that environmental factors rather than climate are critical for the persistence of the hybrids. More studies are needed to elucidate the contributions of environmental factors to fitness and adaptation in these two watermilfoil species and their hybrids.
Acknowledgments
We thank Nikos Andreakis, Antoine Guisan and Sarah Goslee for their recommendations and guidance in the data analyses and members of Dan Yu's group for field assistance. This study was supported by grants from the National Natural Science Foundation of China to X. X. (31070190 and 31270265) and to D. Y. (30930011).
References
Author notes
Data deposited at Dryad: doi: 10.5061/dryad.c418b
