Original Research ARTICLE
Genetic Differentiation of Eastern Honey Bee (Apis cerana) Populations Across Qinghai-Tibet Plateau-Valley Landforms
- 1College of Life Sciences, Fujian Agriculture and Forestry University, Fuzhou, China
- 2College of Bee Science, Fujian Agriculture and Forestry University, Fuzhou, China
- 3Tibet Academy of Agricultural and Animal Husbandry Sciences, Lhasa, China
- 4Ganzi Tibetan Autonomous Prefecture Apiculture Management Station, Ganzi, China
- 5Sichuan Province Apiculture Management Station, Chengdu, China
Many species of high-altitude plateaus tend to be narrowly distributed along river valleys at lower elevations due to a limitation of suitable habitats. The eastern honeybee (Apis cerana) is such a species and this study explored the effects of long and narrow geographic distributions on honeybee populations. Genetic differentiation and diversity were assessed across populations of the southeastern Qinghai-Tibet Plateau. A total of 492 honeybee samples from eight sampling sites in four valleys were analyzed for the genetic differentiation and diversity of 31 microsatellite loci and mitochondrial tRNAleu-COII fragments. The following results were obtained: (1) Microsatellite genetic differentiation coefficients (FST) ranged from 0.06 to 0.16, and mitochondrial FST estimates ranged from 0.18 to 0.70 for different sampling sites in the same valley, indicating genetic differentiation. (2) Honeybees in adjacent valleys were also genetically differentiated. The FST of microsatellites and mitochondria were 0.04–0.29 and 0.06–0.76, respectively. (3) Likely a result of small population sizes, the observed genetic diversity was low. The observed impedance of honeybee gene flow among valleys increased both genetic differentiation and population numbers in the Qinghai-Tibet Plateau. This study contributes significantly to the current understanding of the mechanism underlying population genetic differentiation and highlights the potential effects of utilizing genetic resources that are subject to the ecological conditions of the long and narrow geographic distributions of plateau-valley landforms.
Genetic differentiation of populations is the theoretical basis of the study of species formation, biological evolution, and genetic resources. In response to the blocked gene flow, evolutionary forces, such as natural selection and genetic drift, promote genetic differentiation among populations. Adaptation to their respective environments leads to the development of unique genetic characteristics and traits in differentiated populations (Mayr, 1947; Bohonak, 1999; Friesen et al., 2007). As the basis of new species formation, the genetic differentiation of populations plays an important role in biological evolution, as well as for the study of the conservation and application of genetic resources. Ecological environments are important causes of gene flow restriction, which results in genetic differentiation (Seddon et al., 2002; Jacquet et al., 2014; Sobierajska et al., 2016). It has been reported that the orchid bee, the Amazonian frog, the tufted deer, the Alpine silver ant, and various amphibians exhibit genetic differentiation in different conditions, including marine (Boff et al., 2014; Warren et al., 2015), riverine (Fouquet et al., 2012; Sun et al., 2016), and mountainous habitats (Purcell et al., 2015). Gene flow is related to both the distribution and range of a given species. The gene flow within a species that is confined to a small geographic range and characterized by little movement is easily blocked, which increases its potential for genetic differentiation (Sun et al., 2016; Sánchez-Montes et al., 2017).
The study of genetic differentiation in the honeybee has shown that the particular characteristics of the ecological environment, including nearby water bodies, distance, elevation, and presence or absence of deserts and mountains, could lead to genetic differentiation. On islands or between islands and continents, honeybee have undergone genetic differentiation. Hainan and Damen Island in China (Zhu et al., 2009; Xu et al., 2013a,b; Zhou et al., 2018), Koh Samui Island in Thailand (Sihanuntavong et al., 1999; Sittipraneed et al., 2001; Songram et al., 2006), as well as the archipelagos of the Philippines, Indonesia, Malaysia, Cyprus, Malta, Balearic, Madeira, and the Azores are prime examples (Smith and Hagen, 1996; Sheppard et al., 1997; Smith et al., 2000; Rúa et al., 2001, 2006; De la Rúa et al., 2003; Kandemir et al., 2006; Radloff et al., 2010). Many cases of population genetic differentiation have been reported. For example, Apis cerana in the Loess Plateau of northern Shaanxi, China, exhibited population genetic differentiation due to the discontinuous distribution of the populations (Xu et al., 2013c). Similarly, A. cerana populations distributed throughout the low- and high-altitude plateaus of western Sichuan also differentiated (Zhu et al., 2017). In the Sahara Desert, genetic differentiation of the western honey bee (Apis mellifera) was reported (Franck et al., 2001; Shaibi and Moritz, 2010). In the Pyrenees, the European dark bee (Apis mellifera mellifera), and the Spanish bee (Apis mellifera iberiensis), inhabiting the northeast and southwest sides of the mountains, respectively, have evolved into distinct subspecies (Ruttner, 1989; Miguel et al., 2007). In the Central Asia, A. mellifera show genetic differentiation among subspecies as a result of long distances (Sheppard and Meixner, 2003).
A. cerana is naturally distributed throughout the eastern and southern regions of the Qinghai-Tibet Plateau, which is the highest plateau in the world. Here, elevation is an important environmental factor that limits the distribution of species. Limited by the availability of appropriate nesting conditions, low temperatures, and insufficient nectar sources, A. cerana are not found across the plateau at elevations exceeding 3,500 m and are only capable to survive in the valleys of the eastern Qinghai-Tibet Plateau, where the altitude is relatively low (Yang, 2001; China National Commission of Animal Genetic Resources, 2011). For example, a number of A. cerana populations are distributed throughout the valleys of the rivers Palongzangbu, the Jinsha, the Yalong, and the Dadu in the southeastern Qinghai-Tibet Plateau. The plateaus and mountains between these valleys have elevations exceeding 4,000 m (Shu et al., 1994; Li et al., 2011), thus restricting the mixing of A. cerana populations. The long and narrow distribution of A. cerana in this unique plateau landform results in characteristic patterns of gene flow, which is of particular value for the study of genetic differentiation, origin, and evolution of A. cerana populations. Here, using microsatellite and mitochondrial genetic markers, the genetic differentiation of A. cerana populations was examined within and between the valleys of the Qinghai-Tibet Plateau. The observations reported here provide a foundation for the exploration of the origin and evolution of A. cerana and for the search for genetic resources in this species.
Materials and Methods
A. cerana is primarily distributed throughout the valleys at lower elevations of the eastern and southern Qinghai-Tibet Plateau (Table S1). Eight sites were sampled in this study that were distributed among the valleys of the rivers Palongzangbu, Jinsha, Yalong, and Dadu. Bomi is located in the Palongzangbu River Valley (Palongzangbu R. V.) at the southern Tanggula Mountains of the Qinghai-Tibet Plateau. The sampling sites of Batang, Derong, and Diqing are located in the Jinsha River Valley (Jinsha R. V.) between the Shaluli and Mangkang Mountains of the Hengduan Mountain range. The Yajiang, Jiulong, and Muli sites are located in the Yalong River Valley (Yalong R. V.), with tributaries extending toward the Zheduo and Shaluli Mountains of the Hengduan Mountain ranges. The Xiaojin sampling site is situated in the Dadu River Valley (Dadu R. V.) between the Qionglai and Daxue Mountains of the Hengduan Mountain range (Figure 1, Table S1).
Figure 1. Locations of Apis cerana sampling sites on the Qinghai-Tibet Plateau. The numbers 1–8 on the map represent the sites of Bomi, Batang, Derong, Diqing, Yajiang, Jiulong, Muli, and Xiaojin, respectively.
A total of 492 A. cerana colonies were sampled from eight sampling sites: 38 samples from Bomi (XZBM), 78 samples from Batang (SCBT), 67 samples from Derong (SCDR), 50 samples from Diqing (YNDQ), 65 samples from Yajiang (SCYJ), 61 samples from Jiulong (SCJL), 67 samples from Muli (SCML), and 66 samples from Xiaojin (SCXJ). To eliminate the influence of captive bee colonies, all A. cerana samples were collected from wild colonies without captively bred queens (Figure 2).
Figure 2. Apis cerana samples collected from traditional colonies. (A) Hollow Tilia with traditional honeycomb. (B) Traditional honeycomb in a hollow Tilia. (C) Crosscut honeycomb in Tilia growing on the Qinghai-Tibet Plateau. (D) Traditional bee colony honeycomb in a crosscut Tilia.
DNA Extraction and Amplification
A single bee was collected from each bee colony to analyze microsatellite and mitochondrial diversities. Genomic DNA was extracted from the thorax of worker bees using the Ezup column genomic DNA extraction kit (Sangon Biotech Co., Ltd., Shanghai, China). Sample tissues were digested using protease K, extracted with phenol, then cleaned up using a column, and finally, DNA was dissolved in 60 μl CE buffer (pH 9.0) from the kit.
A total of 31 microsatellite loci were amplified using multiplex PCR kits (Solignac et al., 2003; Takahashi et al., 2009) (Table S2). The Tiangen multiplex PCR kit was purchased from the Tiangen Biotech Co., Ltd. (Beijing, China), and the Toptaq multiplex PCR kit was purchased from TransGen Biotech Co., Ltd. (Beijing, China). The following thermocycling protocol was used for amplification with the Tiangen multiplex PCR amplification kit: initial denaturation at 95°C for 15 min, 30 cycles of denaturation at 95°C for 30 s, annealing at 58°C for 40 s, and extension at 72°C for 40 s, followed by a final extension at 72°C for 10 min. The following amplification thermocycling profile was used with the Toptaq multiplex PCR kit: initial denaturation at 94°C for 5 min, 30 cycles of denaturation at 94°C for 30 s, annealing at 58°C for 30 s, and extension at 72°C for 30 s, followed by a final extension at 72°C for 10 min. The amplified products were analyzed by capillary electrophoresis with an ABI 3730xl automatic sequencer (Applied Biosystems Inc., Foster City, CA, USA) by the Suzhou Genewiz Biotech Co., Ltd. (Suzhou, China). As internal standard the GeneScan™ LIZ® was used, and electrophoresis data were analyzed using GeneMapper 4.0.
Fragments of mitochondrial sequences (tRNAleu-COII; 1–352 bp) were amplified using the rtaq amplification kit (Takara Biomedical Technology Co., Ltd., Dalian, China). Primer sequences were optimized versions of similar primers that are used to amplify the same sequences in A. mellifera (Garnery et al., 1992). The upstream primer (H2.W2R) sequence was 5′-TCAGGTATTCAGGATCA-3′, and the downstream primer (E.W1F) sequence was 5′-TTTAATATGCAGAATAGTG-3′. The following amplification thermocycling profile was used: initial denaturation at 95°C for 5 min, followed by 30 cycles of denaturation at 95°C for 30 s, annealing at 50°C for 30 s, and extension at 72°C for 30 s, followed by a final extension at 72°C for 8 min. The amplicons were sequenced by BioSune Biotech Co., Ltd. (Fuzhou, China) using an ABI 3730xl sequencer.
A principal coordinate analysis (PCoA) of microsatellite data was conducted by GenAlEx 6.5 (Peakall and Smouse, 2012). The spatial distribution of the principal coordinates was used to evaluate the genetic differentiation of samples, and information on genetic variations was extracted using the distance-standardized method where the first three principal coordinates were selected. To detect the genetic differences between different valleys and within the same valley, AMOVA was conducted on the microsatellite data using GenAlEx 6.5. Discriminant analysis of principal components (DAPC) was conducted on microsatellite data using R version 3.3.2 to detect genetic differentiation among samples (R Core Team, 2016). The samples were grouped a priori, and the principal components, which explained 85% of the total variation, were retained for discriminant analysis. To analyze the FST value (genetic differentiation coefficient) of microsatellite markers among sampled sites, the degree of genetic differentiation among samples was assessed, and statistical significance was determined using GenAlEx 6.5 (Nei, 1977; Peakall and Smouse, 2012). For equivalent analyses of mitochondrial genetic markers, Arlequin 3.5 was used (Excoffier and Lischer, 2010). Sample clustering was inferred using Structure 2.3.4 (Pritchard et al., 2000). The K value in Structure was set to 2–8 clusters, with 10 repetitions, and the admixture model was used. Markov chain Monte Carlo iterations were set to 100,000, and the burn-in was set to 10,000. After the operation, ΔK was assessed using the Structure Harvester program to obtain the best K value (Earl and Vonholdt, 2012), and the results were compiled with the CLUMPAK online software tool (Kopelman et al., 2015).
For microsatellite allele data, expected heterozygosity (He) (Nei and Li, 1979; Nei, 1987), observed heterozygosity (Ho), polymorphic information content (PIC), and allele number (Na) were calculated with Mstools (Park, 2001). The effective allele number (Ne) (Hartl and Clark, 1989) and Shannon index (I) (Verdu, 1998) were calculated using Popgene1.31 (Yeh and Boyle, 1997). Fragments of mitochondrial sequences were cut by Clustal X (Larkin et al., 2007), and the obtained sequence information was uploaded to NCBI to confirm haplotypes. To count the various observed haplotypes and to calculate haplotype diversity (Hd), the average number of nucleotide differences (k), and nucleotide diversity (π), DNAsp 5.0 was used (Librado and Rozas, 2009). The median-joining network of haplotypes was constructed using Network 5 (Fluxus Technology, http://www.fluxus-engineering.com/sharenet.htm). Ne estimator 2.01 was utilized to investigate the effective population size (LDNe), using the linkage disequilibrium method with a random mating model, where the critical value was set to 0.05 (Do et al., 2014).
Genetic Differentiation of A. cerana Within a Single Valley
All A. cerana samples of the same valley exhibited genetic differentiation. In the structure analysis of microsatellite markers, when K = 4 (the best K value), the Yajiang sample was genetically differentiated from Jiulong and Muli samples in the Yalong R. V. Samples from Diqing also showed genetic differentiation from Batang and Derong samples in Jinsha R. V. When K = 8, genetic differentiation was found between Batang and Derong and also between Jiulong and Muli (Figure 3).
Figure 3. Structure analysis of Apis cerana throughout the eight sampling sites. (A) According to the Structure Harvester program results, the best K value for clusters is 4. (B) Results for K = 4 and K = 8 clusters, with colors representing proportions of samples in each of the K inferred clusters. XZBM, Bomi; SCBT, Batang; SCDR, Derong; YNDQ, Diqing; SCYJ, Yajiang; SCJL, Jiulong; SCML, Muli; SCXJ, Xiaojin.
The genetic structure of mitochondrial haplotypes was analyzed for samples from Yajiang, Jiulong, and Muli in the Yalong R. V. The genetic structure of Yajiang samples was dominated by the unique Acmt01022 haplotype (65%). Both Acmt01152 and Acmt01001 haplotypes were found to be dominant in Jiulong and Muli. The Acmt01152 and Acmt01001 haplotypes detected in Jiulong comprised 77 and 10% of the population, respectively, and those detected in Muli comprised 24 and 37%, respectively. Analysis of the Batang, Derong, and Diqing (Jinsha R. V.) samples showed that Batang was composed of three haplotypes, of which Acmt01001 accounted for 94%. Derong included six haplotypes, 61% of which accounted for Acmt01001. Diqing was dominated by the two haplotypes Acmt01152 (65%) and Acmt01215 (20%) (Table 1). The mitochondrial network structure indicated that all Bomi haplotypes were unique, while haplotypes from the other sample sites were different with the exception of Acmt01001, Acmt01003, Acmt01152, and Acmt01215 (Figure 4, Table S3).
Table 1. Genetic diversity of Apis cerana throughout the investigated river valleys of the Qinghai-Tibet Plateau.
Figure 4. Median-joining network of mitochondrial tRNAleu-COII haplotypes in Apis cerana populations. Each circle represents one particular haplotype. Circle size reflects the number of individuals with that haplotype (not to scale); 388609 represents haplotype DQ388609, and the remaining three digits indicate the abbreviations of the name for a particular haplotype. For example, “265” indicates the Acmt01265 haplotype. Populations are distinguished by different colors. XZBM, Bomi; SCBT, Batang; SCDR, Derong; YNDQ, Diqing; SCYJ, Yajiang; SCJL, Jiulong; SCML, Muli; SCXJ, Xiaojin.
In the analysis of microsatellite FST, Yajiang, Jiulong, and Muli (Yalong R. V.) FST values were 0.06–0.16, while FST values of Batang, Derong, and Diqing in the Jinsha R. V. were 0.07–0.10. Analysis of mitochondrial FST of Yajiang, Jiulong, and Muli in the Yalong R. V. found FST values of 0.23–0.54, and FST values of Batang, Derong, and Diqing in the Jinsha R. V. were 0.18–0.70. Samples collected from within a single valley exhibited significant genetic differentiation, indicating that genetic differentiation could occur within the same valley (Table 2).
Table 2. Microsatellite FST and mitochondrial FST values of Apis cerana between different sampling sites.
Microsatellite AMOVA was conducted to verify genetic differentiation within a particular valley. The analysis indicated that the genetic variation among samples from the same valley represented 9% of the total variation. The differences were significant (P < 0.01), indicating the existence of high levels of genetic differentiation among bees of the same valley (Table S4).
Genetic Differentiation of A. cerana Between Valleys
Genetic differentiation was found between sampling sites of different valleys. PCoA of the microsatellite showed that the first and second principal coordinates accounted for 77.18% of the variation (51.26 and 25.92%, respectively). For principal coordinate 1, samples collected from Xiaojin (Dadu R. V.) were clearly different from other samples. The principal coordinates 1 and 2 of samples collected from Bomi (Palongzangbu R. V.) were different from other samples, and samples from Batang, Derong, and Diqing (Jinsha R. V.) differed from those of Yajiang, Jiulong, and Muli (Yalong R. V.) with respect to principal coordinate 2 (Figure 5). In the microsatellite DAPC, the discriminant functions (1 and 2) for Bomi (Palongzangbu R. V.) were clearly different from other samples, and those of Xiaojin (Dadu R. V.) were distinct as well. Furthermore, the discriminant function 1 of samples collected from Batang, Derong, and Diqing (Jinsha R. V.) differed from those of Yajiang, Jiulong, and Muli (Yalong R. V.) (Figure 6).
Figure 5. Scatterplot of principal coordinates 1 and 2 for Apis cerana among different sampling sites. Populations from different river valleys are indicated with different colored symbols. XZBM, Bomi; SCBT, Batang; SCDR, Derong; YNDQ, Diqing; SCYJ, Yajiang; SCJL, Jiulong; SCML, Muli; SCXJ, Xiaojin.
Figure 6. Discriminant analysis of principal components (DAPC) for Apis cerana among different sampling sites. Populations from different river valleys are indicated with different color symbols. XZBM, Bomi; SCBT, Batang; SCDR, Derong; YNDQ, Diqing; SCYJ, Yajiang; SCJL, Jiulong; SCML, Muli; SCXJ, Xiaojin.
Although introgression was found in A. cerana colonies (No. 1–16) sampled at the Xiaojin site when K = 4 in the Structure analysis of microsatellite markers, genetic differentiation was observed between Xiaojin (Dadu R. V.), Yajiang (Yalong R. V.), Bomi (Palongzangbu R. V.), and other sampling sites. Similarly, the sampling sites of Batang and Derong (Jinsha R. V.) were also genetically different compared to other samples (Figure 3).
The genetic structures of mitochondrial haplotypes from Bomi (Palongzangbu R. V.) and Xiaojin (Dadu R. V.) were dominated by the unique haplotypes Acmt01308 (74%) and Acmt01290 (56%), respectively. In populations from Yajiang, Jiulong, and Muli (Yalong R. V.), Acmt01001 (27%), Acmt01022 (22%), and Acmt01152 (32%) were predominant haplotypes. Furthermore, Acmt01001 (59%) and Acmt01152 (16%) were dominant haplotypes in samples collected from Batang, Derong, and Diqing (Jinsha R. V.) (Table 1).
Microsatellite and mitochondrial FST values were calculated to determine the degree of genetic differentiation of samples. Samples from Xiaojin (Dadu R. V.) and other valleys had microsatellite FST values of 0.18–0.29. The microsatellite FST values between Bomi (Palongzangbu R. V.) and other populations were 0.10–0.24, and the microsatellite FST values of samples collected from Batang, Derong, and Diqing (Jinsha R. V.) and those from Yajiang, Jiulong, and Muli (Yalong R. V.) were 0.04–0.20. The mitochondrial FST values between Xiaojin (Dadu R. V.) and the other valleys were 0.26–0.59. Comparison between samples collected from Bomi (Palongzangbu R. V.) and others had FST values of 0.38–0.76, and the mitochondrial FST values between Batang, Derong, and Diqing (Jinsha R. V.) and Yajiang, Jiulong, Muli (Yalong R. V.) were 0.06–0.72 (Table 2). All samples of the four valleys showed very significant genetic differentiation, which is indicative of genetic differentiation between different valleys.
Microsatellite AMOVA was conducted to confirm the results of genetic differentiation among different valleys. The genetic variation of the honeybees in the four different valleys accounted for 7% of the total variation. The observed differentiation was significant (P < 0.01), indicating strong genetic differentiation among honeybees in the different valleys (Table S4).
With regard to microsatellite markers, the expected heterozygosity of A. cerana in the valleys of the Qinghai-Tibet Plateau was 0.28–0.40. The number of effective alleles was 1.78–2.60, and the Shannon index was 0.52–0.84. For the tested mitochondrial markers, haplotypes ranged from 3 to 11, with observed haplotype diversities of 0.123–0.761.
Among the eight sampled populations, the effective population sizes of Diqing and Muli were infinite, indicating that these particular populations may be sufficiently large that the effect of genetic drift is negligible. However, the exact effective population size remains unknown. The expected heterozygosity of the effective populations in Diqing and Muli was 0.38–0.39, the observed heterozygosity was 0.35–0.37, the PIC value was 0.25–0.37, the allele number was 4.61–5.77, the effective allele number was 2.59–2.60, the Shannon index was 0.79–0.84, and the haplotype diversity was 0.539–0.761. The effective population sizes of Xiaojin, Yajiang, Bomi, Batang, Derong, and Jiulong were all <500, with expected heterozygosities of 0.28–0.40. The actual observed heterozygosities were 0.25–0.36, PIC values were 0.25–0.37, allele numbers were 3.03–4.65, effective allele numbers were 1.78–2.57, the Shannon index was 0.52–0.79, and haplotype diversities were 0.123–0.660 (Table 1).
In this study, the genetic differentiation of A. cerana inhabiting the long and narrow alpine valleys of the Qinghai-Tibet Plateau was analyzed. This will help to understand the population differentiation rules of species under this specific geographical condition, and will furthermore help to discover and protect the genetic resources of honeybees. Due to the characteristic narrowness of plateau canyons, the degree of continuous distribution of organisms throughout canyons determines the level of population differentiation, and also provides the necessary natural conditions for the formation of diverse germplasm resources (Gaudeul and Till-Bottraud, 2008; Raffl et al., 2008). The results of a study of population differentiation of grass snakes (Natrix natrix) in valleys of the Alps showed that the snake population in canyons was large and distributed continuously; therefore, gene flow could be maintained and no genetic differentiation was found (Meister et al., 2012). However, organisms such as ants and marmots show genetic differentiation among populations as a result of long-distance and potential discontinuous distribution (Goossens et al., 2001; Purcell et al., 2015). The results presented here indicate that A. cerana, with its long and narrow distributions along plateau-valley landforms, is prone to population genetic differentiation as a result from impeded gene flow caused by their discontinuous distribution. Therefore, this region was more likely to contain population diversity and diverse genetic resources. However, the obtained results show that the population sizes of these honeybees were small and had low genetic diversity, leading to a high risk of extinction. Therefore, particular attention should be focused on the conservation of honeybees in alpine valleys.
Genetic Differentiation of A. cerana Within the Same Valley
Within the valleys of the Yalong and Jinsha Rivers, three A. cerana sampling sites were sampled per valley, which were separated by distances of approximately 114–176 km and 85–110 km, respectively. The FST values of microsatellites and mitochondria in the samples from the Yalong R. V. were 0.06–0.16 and 0.23–0.54, respectively, and 0.07–0.10 and 0.18–0.70, respectively, for samples collected from the Jinsha R. V. According to Wright (Wright, 1978), FST values of 0–0.05 indicate weak differentiation, FST values of 0.05–0.15 indicate moderate differentiation, FST values of 0.15–0.25 indicate great differentiation, and an FST value exceeding 0.25 indicates extremely great differentiation. Therefore, FST patterns within the same valleys are indicative of at least moderate differentiation. Based on other A. cerana research, population genetic differentiation was found when FST was 0.00427–0.51852 based on mitochondrial markers and 0.06–0.12 based on microsatellite markers (Xu et al., 2013a; Yin and Ji, 2013). This suggests that the A. cerana populations in the same valley have genetically differentiated. The restricting environment of plateau-valley landforms prevents the spread of A. cerana to both sides of a particular valley. A. cerana populations are therefore only able to reside in narrow valleys with their relatively limited space. This results in the emergence of small, discrete populations in the valleys. Furthermore, the elevation of the valley floor is approaching the upper limit of the environmental range of A. cerana. With increasing elevation, the numbers of tall trees suitable for nesting as well as both nectar and target pollen-producing plants (i.e., food sources) decrease. These ecological factors limit population sizes, especially upstream of valleys in the areas of Batang, Derong, Yajiang, and Jiulong. In these areas, the observed populations consisted of <500 colonies. Due to these factors, A. cerana is vulnerable to catastrophic natural events and environmental fluctuations, resulting in discontinuous distribution. The particular distribution of A. cerana makes them prone to genetic differentiation. These characteristics are of great value for an improved understanding of the process and mechanisms that contribute to population differentiation of other species in similar environments.
Genetic Differentiation of A. cerana in Different Valleys
When the movement of individuals is blocked by plateaus that are detrimental to A. cerana survival, populations in adjacent valleys are more likely to undergo genetic differentiation as a result of the blockage of gene flow. This study indeed showed genetic differentiation between A. cerana populations in different valleys of the Qinghai-Tibet Plateau. Based on the conducted PCoA and DAPC, A. cerana in different valleys showed obvious differentiation. Furthermore, the FST of microsatellites and mitochondria ranged from 0.04–0.29 to 0.06–0.76, respectively. According to Wright (Wright, 1978) and compared to other A. cerana research (Xu et al., 2013a; Yin and Ji, 2013), these results indicate population genetic differentiation. The highest reported elevation for a A. cerana population in the literature is 3250 m (Hepburn et al., 2001; Yang, 2001; Radloff et al., 2005). The observations of A. cerana at 3,040 m reported here are consistent with previously reported observations (Zhu et al., 2017). Based on this information, it could be inferred that the upper limit of the hospitable zone for A. cerana does not exceed 3,500 m. The valleys that were chosen in this study are separated by mountains with elevations > 4,000 m, such as the Boshula and Taniantaweng Mountains (Liu et al., 2016; Yang et al., 2016). Two main environmental characteristics are specific for this area. The first is that the area lacks vegetation, and is mostly covered by bare rock. The second is the presence of plateau meadows. In the eastern valleys of the Qinghai-Tibet Plateau, the lack of tree holes for nesting prevents A. cerana survival, thus blocking gene flow throughout the valleys. The resulting genetic differentiation between valleys suggests that nesting conditions are important ecological factors for A. cerana. Between valleys, nectar and pollen plants suitable as food sources for A. cerana have been found, along with bumblebees nesting in the ground. However, the apparent lack of suitable nesting places, such as caves or holes in tall trees, prevents the survival of A. cerana populations in such environments.
Genetic Diversity and Resource Conservation
Comparison of A. cerana in this study with other A. cerana shows their genetic differentiation, which reflects the special genetic structure and potential as germplasm resource of A. cerana in the alpine valley. The obtained samples show genetic differentiation with A. cerana from the Loess Plateau, the Qinling-Daba Mountains, and the Hainan Island as indicated by the FST value with an average is 0.14 in both utilized loci (Table S5) (Xu et al., 2013a,c; Guo et al., 2016). Similarly, the FST values between the samples of the current study and A. cerana from Changbai Mountains and Fujian Province ranged from 0.31 to 0.72 (with an average of 0.45) (Zhu et al., 2011; Yu et al., 2013). Genetic differentiation between the investigated samples and A. cerana in Guizhou is corroborated by FST values (with an average of 0.08) (Yu et al., 2017). Moderate or strong genetic differentiation was found in loci Ap085, AP313, Ac-2, Ac-5, Ac-26, Ac-1, Ac-35, UN117, SV039, BI314, K0715, AP243, AP066, AC011, AP189, BI225, UN244T, and AT004, which indicates that the investigated sample has a distinct genetic structure in these loci. These analyses indicate the special genetic structure of honeybees in the alpine valleys of the Qinghai-Tibet Plateau, which is a consequence of selection and genetic drift influenced by long-time isolation. Therefore, A. cerana in the valleys of the Qinghai-Tibet Plateau is a unique and precious genetic resource.
In the valleys of the Qinghai-Tibet Plateau, the gene flow of A. cerana is easily blocked, resulting in genetic divergence among populations. The diversity of these populations is relatively high; however, the genetic divergence between populations is low. Comparison with similar A. cerana research indicates that He ranges from 0.2066 to 0.8305 (Chen et al., 2011; Ji et al., 2011), PIC ranges from 0.28 to 0.81 (Cao et al., 2013; Xu et al., 2013a), Na ranges from 1.81 to 10.90 (Ji et al., 2011; Xu et al., 2013c), Hd ranges from 0.171 to 0.905 (Zhou et al., 2012; Ren et al., 2018), and π ranges from 0.00049 to 0.03034 (Zhou et al., 2012; Li et al., 2018). The genetic diversity of each population examined in this study is relatively low. This is mainly a result of the effect of small population sizes (Xu et al., 2013b; Zhao et al., 2017). The environment of the Qinghai-Tibet Plateau determines the natural distribution of A. cerana. Here, ecological factors such as elevation and nesting environments have approached the limits of suitability for this species. Under specific environmental constraints, the effective population size of A. cerana is generally lower than 500 colonies. At a given mutation rate, this results in low genetic diversity due to the small numbers of individuals in any given population (Vrijenhoek, 1997; Amos and Harwood, 1998; Frankham et al., 2002; Ellis et al., 2006). The ecological environments of the valleys restrict the expansion of small A. cerana populations into larger populations, with results in low genetic diversity.
The populations of A. cerana in the Qinghai-Tibet Plateau are thus prone to genetic differentiation and exhibit a high level of population diversity. Therefore, the various populations need to be protected and may be utilized as genetic resources. However, both the effective population size and genetic diversity of A. cerana in the valleys are small, making the populations, and the species within this range in general, vulnerable to effects of genetic drift and selection. Potentially, this can result in local extinctions of populations. Consequently, further studies are required to determine how best to protect the genetic resources of these particular A. cerana populations.
This research is not related to biomedical research. In addition, our research object was honeybee, an invertebrate. Commonly, such a study is exempt from ethics approval.
BZ, YY, SZ, XZ, XX, and LH participated in the whole research project, including the experimental design and data analysis. WW, LZ, PW, JW, KL, and SW participated in part of the experimental design and sample collection.
This research was supported by an earmarked fund from the China Agriculture Research System (No. CARS-45-KXJ11) and funding from The Education Department of Fujian Province (No. JT180119). This research was also supported by a Special Fund for Scientific and Technological Innovation in FAFU (No. CXZX2017338, No. CXZX2017339, No. CXZX2017061) and Research Development Fund in FAFU (No. KF2015025).
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
We would like to thank all of the beekeepers who allowed us to collect samples and the many people who helped us to collect samples, specifically Nianzhou Li, Kun Wang, Qing Wang, Xijian Xu, and Zhenbang Hao. We are also grateful to Huiping Guo, Heng Wang, Daoyin Chen, and Kai-chieh Yang, who helped with our experiments and data analysis. Elevation data was obtained from the Data Center of Resources and Environmental Science, Chinese Academy of Sciences.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2019.00483/full#supplementary-material
Boff, S., Soro, A., Paxton, R. J., and Alves-dos-Santos, I. (2014). Island isolation reduces genetic diversity and connectivity but does not significantly elevate diploid male production in a neotropical orchid bee. Conserv. Genet. 15, 1123–1135. doi: 10.1007/s10592-014-0605-0
Cao, L. F., Su, X. L., Zhao, D. X., Hua, Q. Y., and Hu, F. L. (2013). Genetic diversity of microsatellite DNA for Apis cerana cerana in Zhejiang. Apicult. China 64, 10–11. doi: 10.3969/j.issn.0412-4367.2013.02.002
Chen, Y., Zong, C., Yu, L. S., and Ji, T. (2011). Study on microsatellite DNA genetic diversity of Apis cerana cerana in wannan mountain area and Wangxi Da Bie Mountain Area. Apicult. China 62, 8–11. Available online at: http://zgyf.cbpt.cnki.net/WKB3/WebPublication/paperDigest.aspx?paperID=1343d780-e136-44dd-a48a-a2561939eb7a#
De la Rúa, P., Galián, J., Serrano, J., and Moritz, R. F. (2003). Genetic structure of Balearic honeybee populations based on microsatellite polymorphism. Genet. Select. Evol. 35, 1–12. doi: 10.1186/1297-9686-35-3-339
Do, C., Waples, R. S., Peel, D., Macbeth, G. M., Tillett, B. J., and Ovenden, J. R. (2014). NEESTIMATOR v2: re-implementation of software for the estimation of contemporary effective population size (Ne) from genetic data. Mol. Ecol. Resour. 14, 209–214. doi: 10.1111/1755-0998.12157
Earl, D. A., and Vonholdt, B. M. (2012). STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv. Genet. Res. 4, 359–361. doi: 10.1007/s12686-011-9548-7
Ellis, J. S., Knight, M. E., Darvill, B., and Goulson, D. (2006). Extremely low effective population sizes, genetic structuring and reduced genetic diversity in a threatened bumblebee species, Bombus sylvarum (Hymenoptera: Apidae). Mol. Ecol. 15, 4375–4386. doi: 10.1111/j.1365-294X.2006.03121.x
Excoffier, L., and Lischer, H. E. (2010). Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol. Ecol. Resour. 10, 564–567. doi: 10.1111/j.1755-0998.2010.02847.x
Fouquet, A., Ledoux, J.-B., Dubut, V., Noonan, B. P., and Scotti, I. (2012). The interplay of dispersal limitation, rivers, and historical events shapes the genetic structure of an Amazonian frog. Biol. J. Linn Soc. 106, 356–373. doi: 10.1111/j.1095-8312.2012.01871.x
Franck, P., Garnery, L., Loiseau, A., Oldroyd, B. P., Hepburn, H. R., Solignac, M., et al. (2001). Genetic diversity of the honeybee in Africa: microsatellite and mitochondrial data. Heredity 86, 420–430. doi: 10.1046/j.1365-2540.2001.00842.x
Garnery, L., Cornuet, J. M., and Solignac, M. (1992). Evolutionary history of the honey bee Apis mellifera inferred from mitochondrial DNA analysis. Mol. Ecol. 1, 145–154. doi: 10.1111/j.1365-294X.1992.tb00170.x
Gaudeul, M., and Till-Bottraud, I. (2008). Genetic structure of the endangered perennial plant Eryngium alpinum (Apiaceae) in an alpine valley. Biol. J. Linn Soc. 93, 667–677. doi: 10.1111/j.1095-8312.2008.00958.x
Goossens, B., Chikhi, L., Taberlet, P., Waits, L. P., and Allainé, D. (2001). Microsatellite analysis of genetic variation among and within Alpine marmot populations in the French Alps. Mol. Ecol. 10, 41–52. doi: 10.1046/j.1365-294X.2001.01192.x
Guo, H. P., Zhou, S. J., Zhu, X. J., Xu, X. J., Yu, Y. L., Yang, K. J., et al. (2016). Population genetic analysis of Apis cerana cerana from the Qinling-Daba Mountain Areas based on microsatellite DNA. Acta Entomol. Sin. 59, 337–345. doi: 10.16380/j.kcxb.2016.03.011
Hepburn, H. R., Radloff, S. E., Verma, S., and Verma, L. R. (2001). Morphometric analysis of Apis cerana populations in the southern Himalayan region. Apidologie 8, 435–447. doi: 10.1051/apido:2001142
Jacquet, F., Nicolas, V., Colyn, M., Kadjo, B., Hutterer, R., Decher, J., et al. (2014). Forest refugia and riverine barriers promote diversification in the West African pygmy shrew (Crocidura obscurior complex, Soricomorpha). Zool. Scr. 43, 131–148. doi: 10.1111/zsc.12039
Ji, T., Yin, L., and Chen, G. H. (2011). Genetic diversity and population structure of Chinese honeybees (Apis cerana) under microsatellite markers. Afr. J. Biotechnol. 10, 1712–1720. doi: 10.5897/AJB10.753
Kandemir, I., Meixner, M. D., Ozkan, A., and Sheppard, W. S. (2006). Genetic characterization of honey bee (Apis mellifera cypria) populations in northern Cyprus. Apidologie 37, 547–555. doi: 10.1051/apido:2006029
Kopelman, N. M., Mayzel, J., Jakobsson, M., Rosenberg, N. A., and Mayrose, I. (2015). Clumpak: a program for identifying clustering modes and packaging population structure inferences across K. Mol. Ecol. Resour. 15, 1179–1191. doi: 10.1111/1755-0998.12387
Larkin, M. A., Blackshields, G., Brown, N. P., Chenna, R., McGettigan, P. A., McWilliam, H., et al. (2007). Clustal W and clustal X version 2.0. Bioinformatics 23, 2947–2948. doi: 10.1093/bioinformatics/btm404
Li, B., Lai, K., Xian, F. H., Xiang, G. D., Zhou, M., Ouyang, B., et al. (2018). Genetic diversity of Apis cerana cerana based on mitochondrial DNA in tangjiahe national nature reserve, Sichuan, China. J. Sichuan Agric. Univ. 36, 386–391. doi: 10.16036/j.issn.1000-2650.2018.03.017
Li, Z. X., He, Y. Q., Wang, C. F., Wang, X. F., Xin, H. J., Wei, Z., et al. (2011). Spatial and temporal trends of temperature and precipitation during 1960–2008 at the Hengduan Mountains, China. Q. Int. 236, 127–142. doi: 10.1016/j.quaint.2010.05.017
Liu, W., Zhao, Y., You, J., Qi, D., Zhou, Y., Chen, J., et al. (2016). Morphological and genetic variation along a north-to-south transect in Stipa purpurea, a dominant grass on the Qinghai-Tibetan Plateau: implications for response to climate change. PLoS ONE 11, 1–22. doi: 10.1371/journal.pone.0161972
Miguel, I., Iriondo, M., Garnery, L., Sheppard, W. S., and Estonba, A. (2007). Gene flow within the M evolutionary lineage of Apis mellifera: role of the Pyrenees, isolation by distance and post-glacial re-colonization routes in the western Europe. Apidologie 38, 141–155. doi: 10.1051/apido:2007007
Peakall, R., and Smouse, P. E. (2012). GenAlEx 6.5: genetic analysis in Excel. Population genetic software for teaching and research-an update. Bioinformatics 28, 2537–2539. doi: 10.1093/bioinformatics/bts460
R Core Team (2016). R: A Language and Environment for Statistical Computing (Version 3.3.2). Vienna: R Foundation for Statistical Computing. Available online at: https://www.R-project.org/
Radloff, S. E., Hepburn, C., Hepburn, H. R., Fuchs, S., Hadisoesilo, S., Tan, K., et al. (2010). Population structure and classification of Apis cerana. Apidologie 41, 589–601. doi: 10.1051/apido/2010008
Radloff, S. E., Hepburn, H. R., Fuchs, S., Otis, G. W., Hadisoesilo, S., Hepburn, C., et al. (2005). Multivariate morphometric analysis of the Apis cerana populations of oceanic Asia. Apidologie 36, 475–492. doi: 10.1051/apido:2005034
Raffl, C., Holderegger, R., Parson, W., and Erschbamer, B. (2008). Patterns in genetic diversity of Trifolium pallescens populations do not reflect chronosequence on alpine glacier forelands. Heredity 100, 526–532. doi: 10.1038/hdy.2008.8
Ren, Q., Cao, L. F., Zhao, H. X., Wang, R. S., Cheng, S., Luo, W. H., et al. (2018). Analysis of genetic diversity of main Apis cerana populations in China. J. Henan Agric. Univ. 52, 91–95. doi: 10.16445/j.cnki.1000-2340.2018.01.015
Rúa, P. D. L., Galián, J., Bo, V. P., and Serrano, J. (2006). Molecular characterization and population structure of Apis mellifera from Madeira and the Azores. Apidologie 37, 699–708. doi: 10.1051/apido:2006044
Rúa, P. D. L., Galián, J., Serrano, J., and Moritz, R. F. A. (2001). Molecular characterization and population structure of the honeybees from the balearic islands (Spain). Apidologie 32, 417–427. doi: 10.1051/apido:2001141
Sánchez-Montes, G., Wang, J., Ariño, A. H., and Martínez-Solano, Í. (2017). Mountains as barriers to gene flow in amphibians: quantifying the differential effect of a major mountain ridge on the genetic structure of four sympatric species with different life history traits. J. Biogeogr. 45, 318–331. doi: 10.1111/jbi.13132
Seddon, J. M., Santucci, F., Reeve, N., and Hewitt, G. M. (2002). Caucasus Mountains divide postulated postglacial colonization routes in the white-breasted hedgehog, Erinaceus concolor. J. Evol. Biol. 15, 463–467. doi: 10.1046/j.1420-9101.2002.00408.x
Shu, F., Zhao, Y. H., Lei, H., Jin, Y., and Ma, C. Q. (1994). Boundaries and characteristics of arid regions in mountain valleys in Southwestern China. Mount. Res. Dev. 38, 73–84. doi: 10.1659/MRD-JOURNAL-D-17-00010.1
Sihanuntavong, D., Sittipraneed, S., and Klinbunga, S. (1999). Mitochondrial DNA diversity and population structure of the honey bee, Apis cerana, in Thailand. J. Apicult. Res. 38, 211–219. doi: 10.1080/00218839.1999.11101012
Sittipraneed, S., Laoaroon, S., Klinbunga, S., and Wongsiri, S. (2001). Genetic differentiation of the honey bee (Apis cerana) in Thailand: evidence from microsatellite polymorphism. J. Apicult. Res. 40, 9–16. doi: 10.1080/00218839.2001.11101043
Smith, D. R., Villafuerte, L., Otis, G., and Palmer, M. R. (2000). Biogeography of Apis cerana F. and A. nigrocincta Smith: insights from mtDNA studies. Apidologie 31, 265–279. doi: 10.1051/apido:2000121
Sobierajska, K., Boratynska, K., Jasinska, A., Dering, M., Ok, T., Douaihy, B., et al. (2016). Effect of the aegean sea barrier between europe and asia on differentiation in Juniperus drupacea (Cupressaceae). Bot. J. Linn. Soc. 180, 365–385. doi: 10.1111/boj.12377
Solignac, M., Vautrin, D., Loiseau, A., Mougel, F., Baudry, E., Estoup, A., et al. (2003). Five hundred and fifty microsatellite markers for the study of the honeybee (Apis mellifera L.) genome. Mol. Ecol. Notes 3, 307–311. doi: 10.1046/j.1471-8286.2003.00436.x
Songram, O., Sittipraneed, S., and Klinbunga, S. (2006). Mitochondrial DNA diversity and genetic differentiation of the honeybee (Apis cerana) in Thailand. Biochem. Genet. 44, 256–269. doi: 10.1007/s10528-006-9030-5
Sun, Z., Pan, T., Wang, H., Pang, M., and Zhang, B. (2016). Yangtze River, an insignificant genetic boundary in tufted deer (Elaphodus cephalophus): the evidence from a first population genetics study. PeerJ 4, 1–21. doi: 10.7717/peerj.2654
Takahashi, J., Shimizu, S., Koyama, S., Kimura, K., Shimizu, K., and Yoshida, T. (2009). Variable microsatellite loci isolated from the Asian honeybee, Apis cerana (Hymenoptera; Apidae). Mol. Ecol. Resour. 9, 819–821. doi: 10.1111/j.1755-0998.2009.02268.x
Warren, B. H., Simberloff, D., Ricklefs, R. E., Aguilée, R., Condamine, F. L., Gravel, D., et al. (2015). Islands as model systems in ecology and evolution: prospects fifty years after MacArthur-Wilson. Ecol. Lett. 18, 200–217. doi: 10.1111/ele.12398
Xu, X. J., Zhou, S. J., Zhu, X. J., and Zhou, B. F. (2013a). Microsatellite DNA analysis of genetic diversity of Apis cerana cerana in Hainan Island, southern China. Acta Entomol. Sin. 56, 554–560. doi: 10.16380/j.kcxb.2013.05.012
Xu, X. J., Zhou, S. J., Zhu, X. J., and Zhou, B. F. (2013c). Microsatellite DNA genetic diversity of Apis cerana cerana from the Loess Plateau, Northweest China. J. Fujian Agric. For. Univ. 42, 638–642. doi: 10.13323/j.cnki.j.fafu(nat.sci.).2013.06.020
Xu, X. J., Zhu, X. J., Zhou, S. J., Wu, X. D., and Zhou, B. F. (2013b). Genetic differentiation between Apis cerana cerana populations from Damen Island and adjacent mainland in China. Acta Ecol. Sin. 33, 122–126. doi: 10.1016/j.chnaes.2013.02.001
Yang, H., Lin, C. P., and Liang, A. P. (2016). Phylogeography of the Rice Spittle Bug (Callitettix versicolor) implies two long-term mountain barriers in South China. Zoolog. Sci. 33, 592–602. doi: 10.2108/zs160042
Yu, Y. L., Zhou, S. J., Xu, X. J., Zhu, X. J., Yang, K. C., Chen, D. Y., et al. (2017). Genetic diversity and genetic differentiation of Apis cerana in Guizhou Province of southwest China. J. Fujian Agric. For. Univ. 46, 323–328. doi: 10.13323/j.cnki.j.fafu(nat.sci.).2017.03.015
Yu, Y. L., Zhou, S. J., Xu, X. J., Zhu, X. J., and Zhou, B. F. (2013). Analysis on genetic diversity of Apis cerana cerana in Changbai Mountains. J. Fujian Agric. For. Univ. 42, 643–647. doi: 10.13323/j.cnki.j.fafu(nat.sci.).2013.06.021
Zhao, W. Z., Wang, M., Liu, Y. Q., Gong, X. Y., Dong, K., Zhou, D. Y., et al. (2017). Phylogeography of Apis cerana populations on Hainan island and southern mainland China revealed by microsatellite polymorphism and mitochondrial DNA. Apidologie 48, 63–74. doi: 10.1007/s13592-016-0450-x
Zhou, S. J., Xu, X. J., Zhu, X. J., Gao, J. L., and Zhou, B. F. (2012). Genetic diversity of Apis cerana cerana in Hainan based on mitochondrial DNA. J. Fujian Agric. For. Univ. 41, 170–175. doi: 10.13323/j.cnki.j.fafu(nat.sci.).2012.02.022
Zhou, S. J., Zhu, X. J., Xu, X. J., Gao, J. L., and Zhou, B. F. (2018). Multivariate morphometric analysis of local and introduced populations of Apis cerana (Hymenoptera: Apidae) on Hainan Island, China. J. Apicult. Res. 57, 374–381. doi: 10.1080/00218839.2018.1455439
Zhu, X. J., Zhou, B. F., Wu, X., Xu, X. J., and Wang, Q. (2009). Morphometric genetic analysis: population differentiation of Apis cerana cerana in Damen Island. Apicult. China 60, 5–7. doi: 10.16380/j.kcxb.2011.05.002
Zhu, X. J., Zhou, B. F., Wu, X. D., Xu, X. J., and Xia, X. C. (2011). Genetic diversity of Apis cerana cerana in Fujian based on microsatellite markers. J. Fujian Agric. For. Univ. 40, 407–411. doi: 10.3969/j.issn.0412-4367.2009.01.001
Zhu, X. J., Zhou, S. J., Xu, X. J., Wang, J. W., Yu, Y. L., Yang, K.-C., et al. (2017). Morphological differentiation in Asian honey bee (Apis cerana) populations in the basin and highlands of southwestern China. J. Apicult. Res. 56, 1–7. doi: 10.1080/00218839.2017.1306374
Keywords: Apis cerana, population genetics, genetic differentiation, Qinghai-Tibet Plateau, population diversity
Citation: Yu Y, Zhou S, Zhu X, Xu X, Wang W, Zha L, Wang P, Wang J, Lai K, Wang S, Hao L and Zhou B (2019) Genetic Differentiation of Eastern Honey Bee (Apis cerana) Populations Across Qinghai-Tibet Plateau-Valley Landforms. Front. Genet. 10:483. doi: 10.3389/fgene.2019.00483
Received: 27 September 2018; Accepted: 06 May 2019;
Published: 22 May 2019.
Edited by:Jianke Li, Institute of Apiculture Research (CAAS), China
Reviewed by:Konstantinos M. Kasiotis, Benaki Phytopathological Institute, Greece
Xiangqian Guo, Henan University, China
Fuliang Hu, College of Animal Science, Zhejiang University, China
Copyright © 2019 Yu, Zhou, Zhu, Xu, Wang, Zha, Wang, Wang, Lai, Wang, Hao and Zhou. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Bingfeng Zhou, firstname.lastname@example.org