Lineage diversity, morphological and genetic divergence in Daphnia magna (Crustacea) among Chinese lakes at different altitudes

The biogeography and genetic structure of aquatic zooplankton populations remains understudied in the Eastern Palearctic, especially the Qinghai-Tibetan Plateau. Here, we explored the population-genetic diversity and structure of the cladoceran waterflea Daphnia magna found in eight (out of 303 investigated) waterbodies across China. The three Tibetan D. magna populations were detected within a small geographical area, suggesting these populations have expanded from refugia. We detected two divergent mitochondrial lineages of D. magna in China: one was restricted to the Qinghai-Tibetan Plateau and the America Moderate to high levels of genetic differentiation within D. magna are often found between different continents and within continents (e.g., De Gelas & De Meester, Fields For example, moderate overall genetic divergence (based on a mitochondrial gene marker) was detected across the European range, and high genetic differentiation was even found on a local scale (De Gelas & De Meester, 2005). By using restriction-site associated DNA sequencing, a clear spatial genetic structure of D. magna was apparent across its Eurasian range, suggesting a geographical distance component in the genetic differentiation of D. magna on a continental scale 2015). Applying synonymous substitutions and microsatellites as genetic markers, higher levels of divergence other was present in lowland China. Several different haplotypes in the Qinghai-Tibetan Plateau were most similar to those from various parts of Siberia, suggesting that as a source region. We also found substantial genetic differentiation between D. magna populations from the Qinghai-Tibetan Plateau and those from lowland China. Moreover, significant morphological differences were identified: D. magna from the Qinghai-Tibetan Plateau had a larger head length, body length and body width than did those from lowland China. Geographical and environmental factors were correlated with the observed morphological variation and genetic divergence of D. magna in China. Our data offer an insight into the divergence of freshwater zooplankton due to the uplift of the Qinghai-Tibetan Plateau.


Introduction
Freshwater invertebrates have frequently been studied to investigate biogeographical principles. Of particular interest has been the genetic diversity across their geographical range, which is often very wide (Taylor et al., 1998) as a consequence of efficient dispersal mechanisms (e.g., Havel & Shurin, 2004;Mayr, 1963). Indeed, many freshwater invertebrate species were assumed in the past to be cosmopolitan because of morphological similarities of specimens inhabiting different continents (Mayr, 1963). However, in-depth morphological analyses, especially when integrated with molecular data, have identified the diagnostic characters of specific lineages and thus led to recognition and description of new species (e.g., Juracka et al., 2010;Kotov et al., 2006;Zuykova et al., 2018b). The application of molecular tools has also demonstrated substantial genetic divergence among populations of freshwater zooplankton taxa (e.g., Rotifera and Cladocera) not only at global (e.g., Petrusek et al., 2004;Xu et al., 2009), but also at regional scale (e.g., Gomez et al., 2000;Ni et al., 2019;Penton et al., 2004).
The many species in the cladoceran genus Daphnia Müller, 1785 (Anomopoda: Daphniidae) are among the most important com-ponents of freshwater ecosystems, being key grazers of phytoplankton as well as main prey for planktivorous fish (Lampert, 2011). Daphnia magna Straus, 1820 has been frequently subjected to ecotoxicological, ecological, evolutionary, biogeographical and physiological studies (e.g., De Gelas & De Meester, 2005;Koussoroplis et al., 2019;Lee et al., 2019;Seyoum & Pradhan, 2019). This species is widely distributed, being detected in Africa, Asia, Europe and North America (e.g., Bekker et al., 2018;Brooks, 1957;Fields et al., 2015;Xu et al., 2018). Moderate to high levels of genetic differentiation within D. magna are often found between different continents and within continents (e.g., De Gelas & De Meester, 2005;Fields et al., 2015). For example, moderate overall genetic divergence (based on a mitochondrial gene marker) was detected across the European range, and high genetic differentiation was even found on a local scale (De Gelas & De Meester, 2005). By using restriction-site associated DNA sequencing, a clear spatial genetic structure of D. magna was apparent across its Eurasian range, suggesting a geographical distance component in the genetic differentiation of D. magna on a continental scale (Fields et al., 2015). Applying synonymous substitutions and microsatellites as genetic markers, higher levels of divergence were found among D. magna populations from northern Eurasia relative to those from southern/central Eurasia (Walser & Haag, 2012). A more recent study explored a deep split between East Asian and Western Eurasian D. magna lineages by using the whole mitochondrial genomes (Fields et al., 2018).
In China, D. magna has been found in the Eastern Plain, the Inner Mongolia-Xinjiang Plateau, the Northeastern Plain and the Qinghai-Tibetan Plateau in the 1970s (Chiang & Du, 1979). This distribution suggests that D. magna can occupy a wide range of habitats, from relatively high-altitude locations in the west (the Qinghai-Tibetan Plateau) to relatively low-altitude regions in the east (lowland China). The Qinghai-Tibetan Plateau, in particular, has a special environment: extremely low temperatures and strong ultraviolet radiation (Clewing et al., 2016;Niu et al., 2019). This unique environment should lead to the adaptive divergence of local species (Favre et al., 2015;Hoorn et al., 2013). A recent study using DNA barcoding evaluated the species diversity of Daphnia from the Qinghai-Tibetan Plateau in China and identified six species (or species complexes), including D. tibetana Sars, 1903, D. longispina species complex, D. magna, D. pulex species complex, D. cf. himalaya Manca, 2006and D. similoides Hudec, 1991(Xu et al., 2018. Most recently, by applying a set of high-resolution microsatellite makers, we found that D. sinensis Gu, 2013 populations from Eastern China and the Qinghai-Tibetan Plateau of Western China were genetically separated from each other (Ma et al., 2019a). However, there have been no specific studies on genetic diversity or population structure of D. magna from the Qinghai-Tibetan Plateau, or even from China overall.
In this study, we analyzed D. magna populations found in eight out of 303 waterbodies sampled across China. These populations differ in their types of habitat, from high-altitude populations in the west (the Qinghai-Tibetan Plateau) to the relatively low-altitude populations in the east (here collectively termed lowland China). First, we investigated morphological variation of D. magna populations from the Qinghai-Tibetan Plateau relative to populations from lowland China. Then, we used the mitochondrial COI marker and a set of nuclear microsatellite loci to investigate the genetic diversity and structure of D. magna populations from China. We expected to detect substantial genetic divergence between D. magna populations from the Qinghai-Tibetan Plateau and those from lowland China, given the significant differences in geographical and environmental factors between these regions.

Sample collection
Daphnia magna samples were recovered from eight of 303 waterbodies across China from 2012 to 2018 (during the growing seasons of Daphnia), using a 125-μm plankton net hauled vertically at two or three sites per locality, from a boat or from the shore. Samples from the same waterbody were pooled and preserved in 95% ethanol. All D. magna specimens were identified morphologically (Benzie, 2005;Chiang & Du, 1979). The eight waterbodies containing this species were two natural lakes, one artificial reservoir and five ponds, representing five main geographical regions of China: Eastern Plain, Inner Mongolia-Xinjiang Plateau, Northeastern Plain, Yunnan-Guizhou Plateau and Qinghai-Tibetan Plateau (table 1 and fig. 1A). For each locality, the following information was collected: geographical position (latitude and longitude), altitude, habitat origin (natural or artificial), surface area, maximum depth, predators (presence or absence of fish) and whether the waterbody freezes in winter (this infor mation was obtained from local residents). Additionally, Downloaded from Brill.com07/11/2020 07:23:36PM via free access the trophic status (eutrophic, mesotrophic or oligotrophic) was assigned based on Chinese literature or information from local environmental protection bureaus (table 1).

Morphometric analyses
The morphological characters of adult female D. magna, including head length (defined as the distance between the upper edge of the head and the lower edge of the head, HL), body width (defined as the distance between the left edge of the body and the right edge of the body, BW) and body length (defined as the distance between the lower edge of the head and the base of the tailspine, BL; supplementary fig. S1), were recorded. Five to 9 adult female D. magna per population were randomly selected for the morphological analyses. The ratios of head length to body length (HL/BL) and body length to body width (BL/BW) were additionally calculated to compensate for size-dependent variations in head length, body width and body length. Comparisons of these measurements between populations were done using one-way ANOVA in Graph-Pad Prism 5 (GraphPad Software, San Diego, CA, USA). A t-test was performed to evaluate the differences in D. magna populations from the Qinghai-Tibetan Plateau versus those from lowland China (combined populations from Eastern Plain, Inner Mongolia-Xinjiang Plateau, Northeastern Plain and Yunnan-Guizhou Plateau), in GraphPad Prism 5. If data did not follow a normal distribution, a Rankit transformation was applied (Conover & Iman, 1981).

DNA extraction
Between 15 and 46 adult females of D. magna per population were randomly selected for the genetic analyses. DNA was extracted from each specimen using 60 μL H3 buffer (Schwenk et al., 1998) with proteinase K (10 mg/ml; MERCK, Germany), containing 10 mM Tris-HCl, 0.05 M KCl, 0.005% Tween 20 and 0.005% NP-40. Individuals were incubated at 55°C for 10-16 hours in a water-bath with mild shaking. The proteinase K was then denatured via a 12 min incubation at 95°C, the tube centrifuged briefly and stored at 4°C before use.

mtDNA sequencing
Six to ten individuals per locality were selected for mitochondrial DNA analysis (i.e., 75 out of 289 individuals; table 1). The mitochondrial cytochrome c oxidase subunit I (COI) gene was amplified using standard primer pairs LCO1490 and HCO2918 (Folmer et al., 1994). The PCR reaction was performed in a total volume of 30 μL, consisting of 3 μL of genomic DNA, 0.5 μM of each primer, 2.5 mM of each dNTP, 10 mM of buffer, and 2 units of Taq HS (TaKaRa, Japan). The PCR protocol was as follows: incubation at 94°C for 1 min, followed by 40 cycles of 1 min at 94°C, 1.5 min at 40°C and 1.5 min at 72°C. This was followed by a final incubation at 72°C for 6 min. The PCR products were then purified and sequenced using the forward primer on an ABI PRISM 3730 DNA capillary sequencer at BGI Tech. (China). Only chromatograms with high quality sequences (Phred quality scores > 40) were used for the genetic analysis. All newly obtained COI sequences from the present study have been submitted to GenBank under accession numbers MN310895-MN310900.

Phylogeny and divergence time estimation
All COI sequences were aligned together with all 655 published COI sequences of D. magna using Clustal W (Thompson et al., 1994) in MEGA 6.0 (Tamura et al., 2013). Published sequences, retrieved from GenBank, were from D. magna populations in Asia, Europe, North Africa and North America (supplementary table S1). Twenty published reference sequences were excluded because of their extremely  Downloaded from Brill.com07/11/2020 07:23:36PM via free access short length. A sequence from a member of the D. similis group was used as an outgroup (GenBank ID: LC389171). The alignment was visually inspected and all the sequences were trimmed to the same length (413 bp). Unique haplotypes were identified in DNASP 5.10 (Librado & Rozas, 2009), and one representative of each was used in phylogenetic analysis. The COI phylogenetic tree was constructed using the Bayesian method in BEAST 2 (Bouckaert et al., 2014), run for 107 generations and a tree recorded every 1,000 generations. Twenty-five percent of trees were discarded as burn-in, and the final 104 trees summarized with Tree-Annotator. Here, HKY+G+I was estimated to be the best-fit substitution model according to the corrected Akaike information criterion (AIC), using jModelTest v. 2.1.7 (Darriba et al., 2012). The tree priors were left at their default values and a strict molecular clock was applied. To ensure that enough generations were computed, Tracer v1.7 (Rambaut et al., 2018) was used. The fossil-based calibration point for the most-recent common ancestor was set at 100 Mya and a 10% standard deviation was assigned. This calibration point was estimated according to the fossil evidence for D. magna and the D. similis group (Cornetti et al., 2019).

Detection of genetic lineages and phylogeographic analyses
The general mixed Yule coalescent model (GMYC, Pons et al., 2006) was applied to the COI phylogenetic tree to identify different lineages within D. magna. The GMYC model, which is a likelihood-based method, uses an ultrametric tree to delimit different species/ lineages, and within-or between-species branching models were fitted to reconstruct gene trees.

Genetic diversity
Only individuals with a complete multilocus genotype (MLG) (i.e., across all 11 microsatellite loci) were used here, resulting in a final dataset of 236 individuals. To explore the level of genetic diversity, relative clonal richness R was calculated per population (R = (G -1) / (N -1), where G is the number of MLGs and N represents sample size; Dorken & Eckert, 2001). Only populations of sample size of 10 or more were included.

Population genetic structure
The identification of MLGs, based on 11 microsatellite loci, was performed in GENALEX 6.5 (Peakall & Smouse, 2006). To distinguish genetically differentiated population groups of D. magna, a factorial correspondence analysis, FCA (GENETIX 4.05, Belkhir et al., 1996Belkhir et al., -2004, was performed. To identify the most likely number of such genetic clusters, a Bayesian algorithm was applied in STRUCTURE V2.3.4 (Pritchard et al., 2000), assuming the existence of K groups. The value of K was set from 1 to 8, ten independent runs were performed for each value of K and, for each run, 105 iterations were carried out after a burn-in period of 105 iterations. The most likely K was determined by the distribution of ΔK, following the methods of Evanno et al. (2005). To display the genetic relationships among the populations, the unweighted pair-group method with arithmetic mean (UPGMA), based on pairwise Nei's genetic distances (Nei, 1978), was applied in MEGA 6.0. An AMOVA test comparable to that applied on COI sequence data (see above) was run on microsatellite data. Pairwise F st values were calculated among populations based on the microsatellite data, the approach being similar to that applied to the COI sequence data. Finally, the correlation between pairwise geographical distance and pairwise F st values (based on microsatellites) was calculated with and without populations from the Qinghai-Tibetan Plateau using a Mantel test.

Environmental preferences
To assess the contribution of geographical and environmental factors in explaining variation between regions (Qinghai-Tibetan Plateau and lowland China), a principal component analysis (PCA) was conducted in Past version 3.12 (Hammer et al., 2001). Next, to relate the presence or absence of D. magna lineages to combinations of these parameters, a PCA was calculated based on geographical position (i.e., longitude and latitude) and seven environmental factors, including altitude, maximum depth, origin (i.e., natural or artificial), predator (i.e., presence or absence of fish), surface area and trophic status (i.e., eutrophic, mesotrophic or oligotrophic) and whether or not the waterbody was frozen in winter. Finally, to determine whether D. magna lineages was non-randomly distributed along each component of the PCA, a one-way analysis of variance (ANOVA) was conducted. All statistical analyses were performed using SPSS 22.0.

Morphology
Significant differences in the head length, body width and body length were found between adult female D. magna from Qinghai-Tibetan Plateau and those from lowland China, see fig. 2. Specifically, adult female D. magna from the Qinghai-Tibetan Plateau were larger in all measured dimensions than those from other regions of China. However, there were no significant differences in the ratios of head length to body length and body length to body width of adult D. magna from the Qinghai-Tibetan Plateau and those from other regions ( fig. 2). Among the populations, DJT had the shortest head length, body width and body length, whereas DL2P had the greatest head length and QGZ had the greatest body width and body length (supplementary fig. S2).

Phylogeny and divergence time estimation
Seventy-five individuals were successfully sequenced at the COI locus (413 bp in the aligned dataset); among them, six unique haplotypes were detected (table 1). The Bayesian tree together with the GMYC lineage-delimitation showed that the six Chinese haplotypes from this study fell in two different mitochondrial lineages ("B2" and "B3"), but belonged to a single group "B". Two haplotypes (i.e., HSD and NLHA), which belonged to the lineage "B2", were restricted to Qinghai-Tibetan Plateau; while the other four haplotypes (i.e., DAHA, HSA, HSB and HSC), belonging to the lineage "B3", were present in lowland China (combined populations from Eastern Plain, Inner Mongolia-Xinjiang Plateau, Northeastern Plain and Yunnan-Guizhou Plateau; fig. 1B). The most recent common ancestor of all included D. magna lineages based on COI was estimated to be around 17.8 Mya (95% HPD: 8.9-27.1 Mya). Daphnia magna from Qinghai-Tibetan Plateau and those from lowland China diverged around 6.6 Mya (95% HPD: 2.5-9.9 Mya; supplementary fig. S3).

Phylogeography
Four Chinese haplotypes (two from the present study and two from GenBank; belonging to "B2") were detected in the Qinghai-Tibetan Plateau: three were restricted to this region, and another one was present in both the Qinghai-Tibetan Plateau of China and West Siberia of Russia. Another six Chinese haplotypes (four from the present study and two from GenBank; belonging to "B3") were detected in lowland China (i.e., mix of Eastern Plain, Inner Mongolia-Xinjiang Plateau, Northeastern Plain and Yunnan-Guizhou Plateau). One haplotype (HSB) was positioned in the center of the star-like haplotype network in the lineage "B3" and shared by four regions of China, including Eastern Plain, Inner Mongolia-Xinjiang Plateau, Northeastern Plain, and Yunnan-Guizhou Plateau ( fig. 3). No haplotype was shared between the Qinghai-Tibetan Plateau and the other regions of China. The AMOVA test based on COI sequences revealed that the genetic variation component at the among-region level was about 10 times higher than at the within-region level (i.e., among populations), or 14 times higher than within populations (i.e., among individuals; Table 2). Pairwise F st values, based on COI, ranged from 0 to 1 among D. magna populations. There was no significant correlation between pairwise F st and geographical distance for D. magna based on COI (R2 = 0.05077, P = 0.249; fig. 4A).

Genetic diversity
There was no evidence for large-allele dropout in any of 11 microsatellite loci in this study. Scoring error due to stuttering might have been present in 2 of the 88 analyzed population/locus cases, and null alleles may be present at loci A001, A002, B031, B075 and B107, as suggested by the general excess of homozygotes for most allele size classes. Overall, D. magna in China had a high clonal diversity    (table 1). Identical MLGs were detected within populations, but never between populations.

Population genetic structure
Assignment tests in STRUCTURE showed that the best estimated number of groups was K = 2, separating Qinghai-Tibetan Plateau and lowland China (combined populations from Eastern Plain, Inner Mongolia-Xinjiang Plateau, Northeastern Plain and Yunnan-Guizhou Plateau; fig. 1C). One has to be cautious in cases with a best K = 2 (Janes et al., 2017). However, specifying other values of K (i.e., K = 3 to K = 8) did not yield biologically reasonable results (i.e., to cluster the populations from lowland China; supplementary fig. S4), as evident from the FCA (as shown in fig. 5). Individuals that were genotyped at up to 11 microsatellite loci were obviously separated by the first FCA axis: one group restricted to the Qinghai-Tibetan Plateau and the other group detected in the remaining four regions (fig. 5). The FCA results were further supported by the UPGMA tree based on genetic distances, in which D. magna populations from the Qinghai-Tibetan Plateau were distinctly separated from those in other regions (supplementary fig. S5). The geographical separation of D. magna populations was further supported by AMOVA: although most of the variation component was within-populations (58.80%), values for among-region (23.59%) and within-regions (17.61%) components were also significant (table 2). Pairwise F st values ranged from 0.14 to 0.54 among D. magna populations (averaged over all loci) based on microsatellite data. There was a significant relationship between pairwise F st and geographical distance for D. magna based on microsatellite markers   (R2 = 0.1822, P = 0.0235; fig. 4B) when all populations were included. However, after excluding populations from the Qinghai-Tibetan Plateau, there was no significant correlation (R2 = 0.07527, P = 0.4430; fig. 4C).

Discussion
In the present study, two divergent mitochondrial lineages of D. magna were detected in China. Specifically, one was restricted to the Qinghai-Tibetan Plateau while the other was present in lowland China. Similarly, microsatellite data demonstrated substantial population-genetic differentiation between D. magna populations from the Qinghai-Tibetan Plateau and those from the lowlands of China. Additionally, a significant morphological difference was found in D. magna populations: individuals from the Qinghai-Tibetan Plateau had a larger body size than did those from lowland China. Such morphological and genetic differences between D. magna populations from the Qinghai-Tibetan Plateau and those from the lowlands were attributed to differences in geographical and environmental factors, e.g., longitude, altitude and trophic status of lake occupied. Phenotypic plasticity, such as predatorinduced phenotypic changes (e.g., Rabus & Laforsch, 2011;Tollrian, 1995), have been observed in numerous Daphnia species. For example, Daphnia reproduce at smaller or larger sizes in response to threats from visually hunting (e.g., fish) or gap-limited (e.g., small invertebrates) predators, respectively Riessen, 1999). Temperature can also significantly affect the body shape of D. magna, D. pulex andD. longispina Muller, 1776 (Ranta et al., 1993). Therefore, the morphological variation we observed between Daphnia populations from the Qinghai-Tibetan Plateau and those from lowland China might be due to phenotypic plasticity in response to the different environmental conditions, such as temperature and predator pressures. Alternatively, the different morphology of high altitude D. magna populations, when compared with those from lowland China, might be due to the unique geographical and environmental conditions there (e.g., Chen et al., 2013;Lin et al., 2008). Individual D. magna with longer heads, wider and longer bodies are more likely to survive and reproduce in the harsher environments (i.e., extremely low temperature and strong ultraviolet radiation) of the Qinghai-Tibetan Plateau. Indeed, the frog Rana kukunoris Xie, 2000 in the Qinghai-Tibetan Plateau produces larger eggs than do those from the lowlands, which allows for increased rates of embryonic development leading to earlier hatching of tadpoles (Chen et al., 2013). The larger females of R. kukunoris in high-altitude habitats produce more gelatinous matrix surrounding the eggs, which could protect embryos from temperature fluctuations and ultraviolet radiation (Chen et al., 2013).
In agreement with a recent study (Bekker et al., 2018), our data suggested that the included clades of D. magna diverged long ago. Both studies used a calibration point estimated based on fossil evidence for Daphnia (Kotov & Taylor, 2011). However, this divergence date differs significantly from findings by Fields et al. (2018). The latter authors applied a fully informative Gaussian prior to substitutions at the third codon position of four-fold degenerate codons. The substantial differences in the estimated coalescence time of D. magna clades is due to the different time calibrations used. Appropriate paleontological records will be required to resolve these differences.
Highly divergent lineages within species are often observed in Cladocera (e.g., Bekker et al., 2018;Ni et al., 2019;Woltereck, 1920). For example, using mitochondrial and nuclear genetic markers, multiple lineages were detected within Moina species across China (Ni et al., 2019). Another study suggested that dark-colored populations of D. pulicaria Forbes, 1893 from the Alps had at least two genetically divergent lineages: one genetically close to boreal haplotypes and another related to refugial populations that survived in southern Europe (Bellati et al., 2014). We found two divergent mitochondrial lineages of D. magna in China with a clear geographical separation: one was restricted to the Qinghai-Tibetan Plateau while the other was present in lowland China. Such geographically distinct phylogroups (based on part or all of the mitochondrial genome) within D. magna were also observed in Europe (De Gelas & De Meester, 2005), in Northern Eurasia ( Bekker et al., 2018) and even between East Asia and Western Eurasia (Fields et al., 2018). The observed phylogenetic structure of D. magna could be explained by strong founder effects that are persistent through time (De Meester et al., 2002), by genetic drift (De Gelas & De Meester, 2005), or local adaptation (Boersma et al., 1999). Indeed, a strong founder effect has been recently demonstrated in high-altitude D. longispina populations from the Altai and Sayan Mountains (Zuykova et al., 2019).
By applying three independent statistical approaches (i.e., FCA, UPGMA and STRUC-TURE) to microsatellite data, we discovered a clear geographical separation between D. magna populations from the Qinghai-Tibetan Plateau and those from lowland China. The geographical distance and differences in ecology and environment might explain the genetic differences. Specifically, lowland D. magna mostly occur at lower altitudes in eutrophic lakes, whereas D. magna populations in the Qinghai-Tibetan Plateau are found at higher altitudes in oligotrophic lakes. Many studies have shown substantial genetic divergence in plant and animal populations between lowland China and the Qinghai-Tibetan Plateau (reviewed in Qiu et al., 2011;Yang et al., 2009). The rapid uplift of the Qinghai-Tibetan Plateau must have had profound effects on the evolution of resident populations and created challenging environments for immigrants, especially for freshwater zooplankton. Altitudinal gradients are characterized by steep changes in physical factors, such as temperature, oxygen concentration and UV radiation, inducing strong diversifying selection. For example, the saker falcon Falco cherrug Gray, 1834 shows obvious genetic divergence between populations from the Qinghai-Tibetan Plateau and those inhabiting grasslands with a lower altitude (Pan et al., 2017). Hypoxia and low temperatures were believed to be the main forces driving such genetic divergence (Pan et al., 2017). Interestingly, we only detected eight D. magna populations from a large number of lakes (303) sampled across China, and the three Tibetan D. magna populations were detected within a small geographical area from 92 lakes we sampled from the Qinghai-Tibetan Plateau. It could be that these populations on the Qinghai-Tibetan Plateau have expanded from Pleistocene refugia and have successfully recolonized the area following the last glacial maximum. Moreover, several different haplotypes were found in the Qinghai-Tibetan Plateau, and they were most similar to those from various parts of Siberia, suggesting that as a source region. However, more intensive sampling in the Qinghai-Tibetan Plateau will be required to address this issue in the future.
Interestingly, there was a significant association between genetic distance and geographical distance for D. magna populations from China. A similar result was found by Fields et al. (2015), who applied RAD sequencing to a large sample of D. magna from Eurasia. But this was not the case in the central/ southern Eurasian D. magna populations (Walser & Haag, 2012). In addition to isolationby-distance, ecological differences among habitats, such as differences in trophic status, could also affect the genetic composition of the Daphnia assemblages (Keller et al., 2008). Therefore, a pattern of isolation-by-distance will not always be a rule, because of the selective constraints. Here, one haplotype (HSB) was shared by D. magna populations in all regions except the Qinghai-Tibetan Plateau. This haplotype was also present in Northeastern Siberia, Southeastern Siberia and Mongolia, being the center of the star-like network of lineage "B3". This star-like network pattern suggests that recent dispersal and colonization events of D. magna have occurred in Asia, and that haplotype HSB was ancestral in the lineage "B3". Similar evidence of rapid expansion has also been observed in other Daphnia species in Asia, such as D. mitsukuri Ishikawa, 1896 (Ma et al., 2019b) and D. galeata Sars, 1864 (Ma et al., 2015). In agreement with previous observations in other Daphnia species, for example, D. galeata (Thielsch et al., 2009;Yin et al., 2018) and D. mitsukuri (Ma et al., 2019b), we found that D. magna populations, regardless of their origin, had high relative clonal richness values (average R = 0.93), suggesting frequent sexual reproduction. However, high genetic diversity is not always a rule in Daphnia, as different species have different demographic histories (e.g., Zuykova et al., 2018a).
In conclusion, we detected substantial genetic divergence (at both mitochondrial DNA and microsatellite markers) and morphological variation between D. magna populations from the Qinghai-Tibetan Plateau and those from lowland China. However, this result needs to be interpreted with caution because the small number of populations and microsatellite loci in this study might limit the power of analysis. Geographical and environmental factors (e.g., altitude and trophic status) were likely to have contributed to this. Our study calls for future investigations on the local adaption of zooplankton to the special habitats in the Qinghai-Tibetan Plateau, especially at the genome level.