Diversity of the Bosmina (Cladocera: Bosminidae) in China, revealed by analysis of two genetic markers (mtDNA 16S and a nuclear ITS)

Background China is an important biogeographical zone in which the genetic legacies of the Tertiary and Quaternary periods are abundant, and the contemporary geography environment plays an important role in species distribution. Therefore, many biogeographical studies have focused on the organisms of the region, especially zooplankton, which is essential in the formation of biogeographical principles. Moreover, the generality of endemism also reinforces the need for detailed regional studies of zooplankton. Bosmina, a group of cosmopolitan zooplankton, is difficult to identify by morphology, and no genetic data are available to date to assess this species complex in China. In this study, 48 waterbodies were sampled covering a large geographical and ecological range in China, the goal of this research is to explore the species distribution of Bosmina across China and to reveal the genetic information of this species complex, based on two genetic markers (a mtDNA 16S and a nuclear ITS). The diversity of taxa in the Bosmina across China was investigated using molecular tools for the first time. Results Two main species were detected in 35 waterbodies: an endemic east Asia B. fatalis, and the B. longirostris that has a Holarctic distribution. B. fatalis had lower genetic polymorphism and population differentiation than B. longirostris. B. fatalis was preponderant in central and eastern China, whereas B. longirostris was dominated in western China. The third lineage (B. hagmanni) was only detected in a reservoir (CJR) of eastern China (Guangdong province). Bosmina had limited distribution on the Tibetan plateau. Conclusions This study revealed that the biogeography of Bosmina appear to be affected by historical events (Pleistocene glaciations) and contemporary environment (such as altitude, eutrophication and isolated habitat). Electronic supplementary material The online version of this article (10.1186/s12862-019-1474-4) contains supplementary material, which is available to authorized users.


Background
China is an important biogeographical zone, where the genetic legacies of the Tertiary and Quaternary periods are abundant. On the one hand, the region (~25-50°N) in the warm temperate-subtropical zone retained moderately high temperatures and provided a good shelter environment for endemics during Tertiary climatic changes [1,2]. On the other hand, the existence of different geographic refugia of Pleistocene (an era in the Quaternary period) did shelter varied lineages with diversity that are also affected by vicariance processes and allopatric speciation [3,4]. At contemporary time, China, a vast territory with different North-South monsoons and East-West altitudes, has a complex geographic environment, the effect of which on the species distribution is also important. Therefore, many biogeographical studies have focused on the organisms of the region [5][6][7][8]. In particular, zooplankton, a group of freshwater invertebrate with a key role in the formation of biogeographical principles, attracted researchers' attention [9][10][11]. Previous genetic analyses have supported the generality of endemism and provincialism, and reinforced the need for detailed regional studies of zooplankton [12,13].
Species of the genus Bosmina Baird, 1845 (Anomopoda: Bosminidae) is small cladocerans (approximation 0. 2540 .319 mm) that exist in lakes and reservoirs in every biogeographical region [14,15]. The genus not only purifies the water by the highly efficient filtration of bacteria, organic detritus and phytoplankton (top-down control) in the waterbodies, but also is rich in protein and fat, which can be used as bait for fish (bottom-up control). Therefore, Bosmina plays an important role in energy flow, the carbon cycle and the nutrient salt mineralization of waters [16]. Moreover, Bosmina may achieve high biomass under a variety of conditions that are unfavorable for other cladocerans. And its dominant position has become increasingly more apparent in recent decades [17]. Bosmina is also among the best preserved zooplankton in limnological sediments due to their robust carapace, making them ideal paleolimnological indicators of ecosystem change [18,19].
Three species of this genus were reported in China based on a morphological study, including B. longirostris O. F. Müller, 1785, B. coregoni Baird, 1857 and B. fatalis Burckhardt, 1924 [20]. Due to the phenotypic plasticity induced by biotic/abiotic signals and the intermediate hybrid morphology that is affected by introgression, assignment of taxa within the Bosmina based on the general phenotype is unreliable [21][22][23]. Fortunately, DNA sequencing of a standard gene region or "DNA barcoding" might speed a solution [24]. Recent studies based on DNA markers have contributed significantly in the resolution of taxonomic units within Bosmina and have allowed questions regarding the diversity and distribution of species to be readdressed in other Holarctic regions (mainly in Europe and in North America) [25,26]. However, no genetic data are available to assess the biogeography of Bosmina in China, which is rich in endemics that are awaiting description [2].
In the present study, zooplankton samples were analyzed from 48 Chinese waterbodies (from 2 to 5135 m above sea level; from oligotrophic to eutrophic conditions), in which 35 populations of Bosmina was detected. Because the mtDNA gene (partial 16S rDNA, further abbreviated as 16S) is highly informative in biogeographic studies, and the nuclear gene region (internal transcribed spacer) is most informative for Bosmina that involve comparisons above the species level [19]. Both fragments were chosen as molecular markers to detect distribution of different lineages for Bosmina and to reveal the DNA molecular information of those lineages found in China.

Sampling methodology
Previous study showed that density of Bosmina is relatively high in summer and positively correlated with water temperature [27], so all samples were collected in summer (June, July and August) of 2016 across China. We sampled waterbodies to cover as much as possible in five major regions, which could fully represent different geographic and limnological types of waterbodies in China. Eastern Plain Region (EPR) is located in the middle and lower reaches of the Yangtze River, where most waterbodies have a high degree of eutrophication. Mongolia-Xinjiang Plateau (MXP) is located in arid and semiarid regions with a sandy/dusty weather and less precipitation. Northeast China region (NCR) is located in a humid and sub-humid area with a continental monsoon climate. Waterbodies in this region are rich in organic materials and humus in sediments. Yunnan-Guizhou Plateau (YGP) is located in the southwest of China, in a large karst area, and there are > 30 waterbodies with a low degree of eutrophication. Tibetan Plateau (QTP) is located at the highest altitude and is the most recently formed region, which is the least polluted area [28]. The sampling situation is as follows: nineteen of these waterbodies were located on the EPR, seven waterbodies were situated on the MXP, four waterbodies lied on the NCR, eight waterbodies were seated on the YGP, and ten remaining waterbodies were distributed on QTP (Table 1).

Environmental variables
Several variables were measured in situ: the sampling position (Latitude, Longitude and Altitude) was recorded with GPS-S6 (Newsmy), Secchi disk visibility (SD) was determined with a Secchi disk M369488 (Beijing Zhongxiyuanda Technology Co., LTD), and water temperature (WT) was measured using FG4-FK (Mettler Toledo Co., Greifensee, Switzerland). Water samples were collected with a Van Dorn bottle for subsequent laboratory analyses of total phosphorus (TP) concentrations, and chlorophyll-a (Chl. a) according to the methods described in detail by Huang [29]. Using the values of SD, TP, and Chl a, a classic trophic state index for waterbodies is calculated to characterize the trophic level of the sampled sites [30].

Zooplankton sampling
Zooplankton samples were collected (Additional file 1: Tables S1, S2 and Fig. 1) with a 125-mm plankton net hauled through the whole water column at three different sites per waterbody within the deep basin or from the shore; these samples were later pooled and preserved in 95% ethanol for further analyses. By using a stereomicroscope, Bosmina was found in 35 waterbodies and ten adult females were selected per waterbody for genetic analyses.

Molecular analysis
DNA of the Bosmina specimens was extracted by the proteinase K method [31] and stored at − 4°C. Each 50-μl polymerase chain reaction (PCR) consisted of 40 μl irradiated H 2 O, 5 μl 10× PCR buffer, 1.5 mM MgCl 2 , 2 mM of each dNTP, 1 μM of each primer, 0.5-1 units Taq DNA polymerase, and 1/10th DNA extract. The primers and PCR temperature profiles of the 16S and ITS followed those of a previous study [19]. Next, the PCR products were purified and sequenced in the forward direction on an ABI PRISM 3730 DNA capillary sequencer by Invitrogen Trading Co., Ltd. (China). The chromatograms were carefully checked for scoring errors.

Sequence alignment
Sequences of the target individuals were aligned with reference species (Additional file 1: Tables S1 and S3) using the Clustal W algorithm [32] at locus 16S and ITS in MEGA.7 [33]. The genera of Bosminopsis and Ilyocryptus are distantly related to Bosmina; thus, the level of divergence between them could lead to an alignment error, long-branch attraction bias and incorrect model estimation. Therefore, in this study, model parameters and phylogenetic trees with in-group taxa only were estimated as suggested by Taylor et al. (2002) [19].

Phylogeny
The phylogenetic relationships among lineages were visualized for both the 16S and ITS datasets. In cases where several individuals carried an identical haplotype, only one was included subsequently. First, the best-fit for both datasets (16S: GTR + G + I and ITS: GTR + G) were identified by jModelTest version 2.1.3 [34], and we then applied these models to run the maximum likelihood (ML) analysis in RAxML [35]. The assignment of species identity (or genetic lineage) of studied individuals was based on clustering with reference sequences in the phylogenetic trees, which was supported by 1000 bootstrapping. All subsequent calculations were based on the taxonomic units that were thereby defined.

Genetic polymorphism and differentiation
The followed analyses were based on data (Additional file 1: Table S1) to learn the genetic polymorphism of two main lineages (i.e. B. fatalis and B. longirostris) found in this study (Fig. 1). In DnaSP v5, haplotype diversity (Hd) and nucleotide diversity (Pi) were calculated for the 16S (249 individuals) and the ITS (305 individuals) genes of main lineages, to determine the level of sequence diversity within the Bosmina [36]. Because 16S is more highly informative for biogeographic studies within species than ITS [19], we calculated p-distances (in MEGA.7) of 16S among all the geographical populations for B. fatalis and B. longirostris respectively, plus standard error calculated on 1000 bootstrap replicates [33]. In addition, the pairwise genetic differentiation of 16S among geographical populations was estimated (in Arlequin v. 3.11) using the fixation index Fst, which included information on mitochondrial haplotype frequency and genetic distances. The  Table S1). Green symbols: B. fatalis; purple symbols: B. hagmanni; red symbols: B. longirostris. Taxon identification is based on 16S and ITS sequence variation significance was tested by 1000 permutations for each pairwise comparison [37].

Haplotype network
To describe the intraspecific variation and relationships among populations from different regions, two haplotype networks (16S and ITS) were constructed based on previous and our studies (Additional file 1: Tables S1 and S3) by HapView with the input of maximum-likelihood trees [38].

Environmental data
Among the 35 sampling sites where Bosmina found, most central and eastern waterbodies (15 of EPR, 2 of MXP, 3 of NCR and 6 of YGP) usually have very high trophic levels, while western waterbodies (2 of EPR, 5 of MXP and 2 of YGP) are mesotrophic or oligotrophic (Table 1).

Phylogeny and distribution
Altogether, 249 individuals were successfully sequenced at the 16S locus, and 305 were sequenced at the ITS region (because of a failure in either amplification or sequencing; Table 1). Among these individuals, 17 unique 16S haplotypes and 12 unique ITS haplotypes were found (Additional file 1: Table S1). Phylogenetic analysis revealed that 291 out of 305 Bosmina individuals sampled in this study definitely belong to two species clades: B. fatalis and B. longirostris. Moreover, one subclade (with 5 individuals) was detected from B. fatalis and another (9 individuals) from the species of B. hagmanni (Figs. 2 and 3, Additional file 1: Table S1). The Bosmina detected in China had different geographic distributions: B. fatalis was present in eutrophic waterbodies from central and eastern China, whereas B. longirostris was dominant in oligotrophic/mesotrophic waterbodies of the western region (Table 1 and Fig. 1). The lineage close to B. fatalis was distributed in four waterbodies, including CAH, CJR, DOH and PYH. The lineage close to B. hagmanni was rarely present in China from this study and was only detected in a reservoir (CJR) of the Guangdong province ( Fig. 1 and Additional file 1: Table S1). No Bosmina was detected in 10 waterbodies on the Tibetan Plateau, indicating its limited distribution in the region ( Fig. 1 and Additional file 1: Table S2).

Genetic polymorphism and differentiation
The level of nucleotide polymorphism within B. fatalis was lower than B. longirostris at the 16S. However, this was not the case for the ITS, for which the level of nucleotide polymorphism between the two species was comparable ( Table 2). The pairwise p-distance showed that most population of B. fatalis had lower differentiation than B. longirostris. In addition, most of the Fst population pairwise comparisons of B. fatalis were low and not significant. However, values of Fst in B. longirostris were generally high and most significant Fst values (p < 0.05) were included. The different Fst also displayed the different differentiation within B. fatalis and B. longirostris (Tables 3 and 4).

Haplotype network
The haplotypes of B. fatalis were restricted to eastern Asia (central and eastern China; Japan). In contrast, those of B. longirostris were widely distributed in other Holarctic regions (China, Denmark, Germany, Russia, Japan and the USA) (Figs. 4 and 5).

Different spatial pattern of B. fatalis and B. longirostris
Our phylogenetic analysis and network showed that the Chinese B. fatalis were closely related to those in Japan, confirming B. fatalis (also known as Sinobosmina fatalis) is endemic in Eastern Asia [39]. It has been early recognized that endemic taxa are not randomly distributed. Range restriction as a consequence of survival in refugia has been regarded as a driving force generating distributional patterns of endemics [12]. East Asia has been recognized as a major glacial Pleistocene refugium [40], where several new endemic cladocerans have found in Japan and neighboring areas [23,41]. It seems that the endemism of B. fatalis is also a result of survival in refugia. The waves of Pleistocene glaciation, resulted in mass extinction of biotas in many region [1], while East Asia region remained comparatively unchanged [40]. Therefore, B. fatalis may survive in this region until now. In addition, the vicariance effects of the Pleistocene glaciations have played a major role in shaping the intraspecific divergence of freshwater zooplankton [12]. A lineage close to B. fatalis found in only four waterbodies (CAH, CJR, DOH and PYH) may be a geographically isolated species differentiated from B. fatalis.
Some cladocerans, such as B. fatalis, may be relicts with restricted-range (endemic), whereas others, commonly and widely distributed within their extensive primary ranges, possiblely are not [2]. We found close relatives among B. longirostris lineages from China, Denmark, Germany, Russia, Japan and the USA, which suggested that B. longirostris is widespread in the Holarctic region [42]. Therefore, B. longirostris should be an advanced lineage and have a non-relict status, which affected by current climate [43]. In all, the different spatial patterns of B. fatalis and B. longirostris can be explained by their different evolutionary status (relict and non-relict).

Genetic polymorphism of B. fatalis and B. longirostris
Due to small population size and past bottlenecks, many relicts tend to have low within-population genetic polymorphism [44,45]. Accordingly, the lower nucleotide diversity of B. fatalis than B. longirostris at 16S locus could be another evidence of their different evolutionary status (Table 2). Also, homogenization of refugial population after glaciation is unlikely in water fleas because they should possess pronounced refugial inertia. That is, successful colonization of an occupied water body is theoretically very difficult because of the strong priority effects from rapid lake-specific selection, massive egg banks and large existing populations [26]. The effect can essentially extend the time of genetic isolation of refugia even in the face of dispersal [46]. This may be why their difference of nucleotide diversity at 16S locus have been preserved to the present. Furthermore, this refugial inertia could provide more time for reproductive isolation [26]. Therefore, after thousands of years of unglaciated conditions, B. fatalis seems remain distinct and did not hybridize with B. longirostris, although they co-occur in some waterbodies (Fig. 1).
However, as a non-coding sequence which do not involve in the formation of the ribosome, ITS have small selective pressure and high divergent rate. For example, ITS evolution rate is 2.5 times of mitochondrial gene (COI) in Tetranychus [47] and 2 times of rDNA in Nasonia [48]. In this study, the nucleotide diversity for B. fatalis is similar to this index for B. longirostris at the   (Table 2). That may be an evolutionary result of ITS between B. fatalis and B. longirostris.
Genetic differentiation of B. fatalis and B. longirostris "Three Gradient Terrains" is a general description of Chinese terrain characters on relief amplitude. It portrays an outline of terrain changes like ladders along West-East direction. The first gradient terrain has the highest topography, with an altitude of over 2,600 m, mainly on the QTP. The second step is between 660 m and 2,600 m above sea level, with many large plateaus and basins. The third stair is below 660 m above sea level. In this study, B. fatalis was on the third stair (Table 1), where hills, low mountains and plains are interlaced [49]. The lower and similar altitude of waterbodies is good for gene communication between different population of zooplankton [9]. That may be the reason for the lower differentiation within B. fatalis than B. longirostris (Tables 3 and 4). On the other hand, it has reported that the geological isolation induced by different altitude can limit gene flow between subpopulation of zooplankton and lead to their differentiation [50]. B. longirostris was distributed on both the first and second steps ( Table 1). The topographic drop between these  steps, about 1000 m, reflects a strong gradient nature and may limit gene flow [49]. That may be responsible for the higher differentiation within B. longirostris (Tables 3 and 4). Due to differing mutation rates and increased rate of transitional saturation rate, 16S will be more highly informative in biogeographic differentiation than ITS. Therefore, the higher differentiation may be shown in 16S, but not in ITS [19]. That should be explained for the different patterns of the haplotype networks obtained for B. longirostris at 16S and ITS genes (Figs. 4 and 5).

Different ecological preponderance of B. fatalis and B. longirostris
A number of physical, chemical and biological factors in the habitat can influence the community structure of zooplankton, and the factors recognized as the most important by majority of investigators were temperature, the quality and availability of food, competition and predation [50][51][52]. Nevertheless, there were no marked differences between the growths of B. fatalis and B. longirostris in response to temperatures ranging from 10°to 30°C in a laboratory experiment [53]. Our temperature data of sampled waterbodies did not exceed this range (13.4°C to 29.8°C) ( Table 1), which indicated that food, competition or predation seemed to more directly affect the population densities of B. fatalis and B. longirostris. The central and eastern waterbodies in our dataset (those on the NCR, YGP and the eastern of MXP and EPR) usually had very high trophic levels, while western waterbodies (those on the western of MXP and EPR) were mostly mesotrophic or oligotrophic and were less affected by human activities (Table 1). Eutrophication often causes an increased density of phytoplankton, especially the bloom of Microcystis [28]. In such condition, B. fatalis can overcome B. longirostris because of more efficient capture and utilization [54].
On the other hand, eutrophication favours cyclopoids, which can prey on Bosmina [55,56]. Previous study has shown that B. fatalis is a superior competitor against B. longirostris and is more resistant to predation [57]. Therefore, competitive edge induced by food and predation in eutrophic waterbodies can be responsible for the preponderance of B. fatalis in the central and eastern Fig. 4 Haplotype networks of the 16S reconstructed in HapView. Circles represent haplotypes, with the sector size being proportional to the number of individuals; the lengths of the connecting lines reflect the number of mutations between them. The multicolored circles with letters a-g stand for shared haplotypes CNm1-7; the monochromatic circles stand for private haplotypes, whose subordinate waterbodies are shown in the legend (Table 1, Additional file 1: Tables S1 and S3) region. On the contrary, dominance of B. longirostris was observed in western waterbodies, which are mesotrophic/oligotrophic.

Appearance of B. hagmanni in China
We only found one lineage close to B. hagmanni in a single locality, the Changjiang Reservoir (CJR) in Guangdong province. B. hagmanni was thought to be confined to the American continent [58]. It seems that an American local lineage of B. hagmanni spread to the new continent, which is not surprising, since previous studies have shown that Bosmina can invade new continents by ship [59]. Guangdong province, which is in frequent maritime trade with America, may provide convenient conditions for alien Bosmina invasion [60]. Nevertheless, the invasion range in China should be quite limited (only present in one reservoir of our samples). A previous study suggested that the interaction between the trophic state and species identity influenced the invasion success of Cladocera [61,62]. Although the waterbodies of Guangdong province are indeed eutrophic [63], we have little knowledge regarding the biological identity of B. hagmanni. Whether eutrophication also restricted the invasion range of B. hagmanni could be an interesting issue for future studies.

The limited distribution of Bosmina on the Tibetan plateau
Although previous reports showed that Bosmina is a cosmopolitan genus of Cladocera, we did not find any species in 10 waterbodies of Tibetan Plateau, where the waterbodies have relatively limited phytoplankton for planktivorous zooplankton [28]. Due to the phenotypic plasticity of the filter screens, Daphnia adapted to the low-food environment [64] and may generate the effective monopoly of food resources, thereby avoiding the spread of other zooplankton [46]. The limited distribution of Bosmina on the Tibetan Plateau may be relative to the extensive distribution of Daphnia in the same region [11]. Moreover, the Tibetan Plateau, as the roof of the world (an elevation of more than 4,000 m, average above sea level) surrounded by high mountains, is ecologically isolated from other regions of China [65]. The more isolated the habitat is, the smaller the probability is that species will colonize it, especially such species as zooplankton that disperse by passive means [66], which may be another reason why Bosmina spread into the Tibetan Plateau finitely.

Conclusion
Based on two genetic markers (mtDNA 16S and a nuclear ITS), we identified individuals from the Bosmina that were sampled from 35 Chinese waterbodies as descendants of two primary species (B. fatalis and B. longirostris) with different spatial pattern, which may be related to their different evolutionary status (relict and non-relict). Their different genetic polymorphism also seems to an inertia of their own evolutionary status. Nevertheless, their different genetic differentiation may Fig. 5 Haplotype networks of the ITS reconstructed in HapView using an ML tree. Circles represent haplotypes, with the sector size being proportional to the number of individuals; the lengths of the connecting lines reflect the number of mutations between them. The multicolored circles with letters h-l stand for shared haplotypes CNn1 to CNn6; the monochromatic circles stand for private haplotypes, whose subordinate waterbodies are shown in the legend (Table 1, Additional file 1: Tables S1 and S3) be related to the altitude of their habitat. Furthermore, B. fatalis was preponderant in central and eastern China, whereas B. longirostris was dominant in western China. Their different ecological dominance may result from competitive interaction induced by food and predation in different trophic waterbodies. Among the waterbodies sampled, the lineage close to alien B. hagmanni was found only at a single locality (CJR), which might be invaded. The limited distribution of Bosmina on the Tibetan Plateau is a feature, which could be explained by resource monopoly of Daphnia in isolated habitat.

Additional files
Additional file 1: Table S1. Numbers of observed haplotypes (of the 16S and ITS) in the investigated Chinese waterbodies;