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).

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).

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®[500] was used, and electrophoresis data were analyzed using GeneMapper 4.0.

Fragments of mitochondrial sequences (tRNA^leu-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 F_ST 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).

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).

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 F_ST 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), F_ST values of 0–0.05 indicate weak differentiation, F_ST values of 0.05–0.15 indicate moderate differentiation, F_ST values of 0.15–0.25 indicate great differentiation, and an F_ST value exceeding 0.25 indicates extremely great differentiation. Therefore, F_ST patterns within the same valleys are indicative of at least moderate differentiation. Based on other A. cerana research, population genetic differentiation was found when F_ST 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.

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 F_ST 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.

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 F_ST value with an average is 0.14 in both utilized loci (Table S5) (Xu et al., 2013a,c; Guo et al., 2016). Similarly, the F_ST 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 F_ST 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.

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.