Association between the Microsatellite Ap243, AC117 and SV185 Polymorphisms and Nosema Disease in the Dark Forest Bee Apis mellifera mellifera

The microsporidian Nosema parasites, primarily Nosema ceranae, remain critical threats to the health of the honey bee Apis mellifera. One promising intervention approach is the breeding of Nosema-resistant honey bee colonies using molecular technologies, for example marker-assisted selection (MAS). For this, specific genetic markers used in bee selection should be developed. The objective of the paper is to search for associations between some microsatellite markers and Nosema disease in a dark forest bee Apis mellifera mellifera. For the dark forest bee, the most promising molecular genetic markers for determining resistance to nosemosis are microsatellite loci AC117, Ap243 and SV185, the alleles of which (“177”, “263” and “269”, respectively) were associated with a low level of Nosema infection. This article is the first associative study aimed at finding DNA loci of resistance to nosemosis in the dark forest bee. Nevertheless, microsatellite markers identified can be used to predict the risk of developing the Nosema disease.


Introduction
In the last few decades, negative processes such as massive losses of bee colonies and hybridization have been observed in honey bee populations worldwide. The honey bee colony losses, called colony collapse disorder (CCD), as a result of reduced adaptation of honey bees to environmental factors, pose a threat to beekeeping worldwide [1,2]. It has been suggested that CCD can be caused by many causes, including various diseases such as nosemosis, environmental pollution, exposure to pesticides, weather and agricultural and beekeeping practices [3][4][5]. On the one hand, in order to avoid a catastrophic population decline from pests and diseases, it is necessary to maintain a high level of genetic diversity in honey bee populations [6,7]. On the other hand, molecular technologies, for example marker-assisted selection (MAS) can be used to identify bee colonies carrying specific traits of interest (e.g., resistance to pathogens and parasites, gentleness and high honey productivity) or the lack of undesirable traits (e.g., aggression and swarming) [8][9][10][11]. Marker-assisted selection is a new technology in beekeeping, and no specific genetic markers that could be used in bee breeding have been proposed [12,13].
The search for informative DNA loci/genes associated with economically useful and other traits (associative mapping) is highly relevant. The preferred strategy is to genotype a high number of genetic markers in linkage or association studies in order to identify genomic regions and discover the causative genes [14]. The identification of genetic markers associated with the phenotype can also immediately be used to selectively breed colonies that are more resistant [15].
Like varroosis, nosemosis is also one of the most dangerous diseases [32][33][34], but the study of associations between molecular genetic markers and nosemosis is rare [35,36].
Nosemosis is a serious disease in adult honey bees caused by microsporidia, which are obligate intracellular eukaryotic parasites [37]. Nosema parasites multiply and develop within the host-cell cytoplasm causing extensive and even total destruction of the midgut epithelial layer [38,39].
Using microsatellite markers, four QTLs associated with low spore load were revealed in a Danish selected Nosema-resistant honey bees line [35]. These Buckfast honey bee colonies have been selectively bred for the absence of Nosema over decades, resulting in a breeding line that is tolerant toward Nosema [59,60]. Unlike pure race breeding, the Buckfast breeding system mixes different stocks to establish a hybrid bee with desired characteristics. The Buckfast contains heritage from mainly A. m. ligustica and A. m. mellifera and from other subspecies. Since the Buckfast bee is a hybrid bee, the expression of its notable characteristics can vary greatly within the stock [61].
It is assumed that the honey bee subspecies (lines and colonies) differ in their resistance to disease, which may be determined by social immunity including hygienic and other types of behavior [13,15,30,[62][63][64][65][66][67][68][69][70]. For example, in A. mellifera, hygienic and grooming behaviors are expressed more highly in Africanized honey bees than in European ones. Perhaps this explains the higher resistance of Africanized bees to V. destructor compared to European bees [71].
Although no significant effect of N. ceranae infections on hygienic behavior was detected [72], it is clear that natural resistance of honey bees to Nosema depends on many factors and the genetic variants of honey bees can play a relevant role. The purpose of this study was to identify associations between genetic variants of some microsatellite loci and Nosema infection/disease resistance in the dark forest bee Apis mellifera mellifera.

Bee Samples
In the present study, samples of a dark forest bee Apis mellifera mellifera from Siberian populations (longitude 81 • 29 −92 • 08 E and latitude 50 • 44 −65 • 47 N) were examined. The dark forest bee is a native bee that was introduced to Siberia about 230 years ago and has adapted well to the local climate and plant communities. In Siberia, the bee population is an artificial population; wintering of bees is controlled by people [73]. Honey bees were collected from 12 apiaries between the end of May 2016 and August of the same year. A total of 226 workers from twenty-eight bee colonies (from 8 to 10 bees from each bee colony) were examined.
For the diagnosis and detection of Nosema infection, the oldest honey bees (forager bees) were collected outside the entrance of hive, because they have the greatest infection and the highest proportion of infected bees [39]. Bee samples were stored in a freezer at −20 • C until further processing.

Study Design
The present research has conducted in several stages. At the first stage of the study, the presence of Nosema spp. in honey bees was investigated using both light microscopy and polymerase chain reaction (PCR).
At the second stage of the study, the genetic diversity of honey bees with different degrees of Nosema infestation was examined using polymorphic microsatellite loci. Earlier, the genetic diversity of local dark forest bees (Siberian population) was studied using a complex of microsatellite loci and identified polymorphic microsatellite loci [74]. In this study, 23 polymorphic microsatellite loci were used to search for genetic markers of Nosema disease resistance in honey bees.
Finally, the associations of polymorphic variants of microsatellite loci studied with nosemosis were analyzed using the odds ratio method (OR).

Experimental Procedures
To carry out individual analysis of the bees, for each sample, the bee's midgut was isolated and divided in two. One part of the midgut was used for light microscopy. For this, the midgut was ground in 0.5 mL of sterile, distilled water and the number of Nosema spores was counted using a Zeiss Axio Lab.1 light microscope.
DNA was extracted from another part of the midgut using a DNA purification kit, PureLink™ Mini (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's protocol. PCR was performed using a thermal MyCycler T100 (BioRad, Foster City, CA, USA).
For the diagnosis of nosemosis, duplex-PCR was used [42]. The primer sequences used to amplify the 321 bp fragment corresponding to the 16S ribosomal gene of N. apis were 321APIS-FOR 5 -GGGGGCATGTCTTTGACGTACTATGTA-3 and 321APIS-REV 5 -GGGGGGCGTTTAAAATGTGAAACAACTATG-3 . The primer sequences utilized to amplify the 218 bp fragment corresponding to the 16S ribosomal gene of N. ceranae were 218MITOC-FOR 5 -CGGCGACGATGTGATATGAAAATATTAA-3 and 218MITOC-REV 5 -CCCGGTCATTCTCAAACAAAAAACCG-3 [42]. PCR was performed in a reaction volume of 20 µL containing 5-10 ng of template DNA, 1 × PCR buffer, 200 µM of each dNTP, 1.5 mM MgCl 2 , 0.2 µM of each primer and 1U Taq polymerase (Fermentas, Thermo Fisher Scientific, Chelmsford, MA, USA). The routine consisted of an initial denaturation step at 94 • C for 2 min, followed by 35 cycles of 94 • C for 30 s, 58 • C for 30 s and 72 • C for 1 min, and a final extension step at 72 • C for 5 min. PCR products were analyzed on 1.5% agarose gels. Gels were stained with ethidium bromide and visualized using UV illumination (Gel Doc XR+, BioRad, Foster City, CA, USA). For each PCR, positive control (reference N. apis and N. ceranae DNA extracts as template) was used. Negative control (ddH2O) was also included in each run of PCR amplification to detect possible contamination.

Estimation of the Nosemosis Level
According to PCR, most of the bees examined were coinfected with two Nosema species, N. apis and N. ceranae. In this regard, we considered the general infestation of honey bees with microsporidia, without division into Nosema species. A light microscope using 400× magnification was used for counting Nosema spores in macerated bee preparations.
Honey bees were divided into the following groups, uninfected (Nosema-negative) and Nosema-infected (Nosema-positive). Since our study did not require a high degree of precision, we used an arbitrary infection scale and divided the Nosema-positive bee group into two variants, Nosema-positive low and Nosema-positive high. Nosema-positive low bees were with a small amount of microsporidial spores on microscopic analysis (less than 100 spores in the field of view of the microscope). Nosema-positive high bees were with a significant number of microsporidial spores detected by microscopic analysis (more than 500 spores in the field of view of the microscope). An intermediate variant of the bee infection (100-500 Nosema spores in the field of view of the microscope) was not identified. In total, three groups of bees were formed: Nosema-negative, Nosema-positive low and Nosema-positive high.

Statistical Analysis
The genotypes obtained for each of the honey bee were used to estimate population parameters. The allelic and genotypic observed frequencies by Hardy-Weinberg equilibrium (HWE), the number of alleles, observed (Ho) and expected (He) heterozygosity were estimated employing the GENEPOP v.4.1 package [76]. To assess genetic diversity, for each infectious category, the observed and expected heterozygosity for each locus were compared using a Student's test. Comparison of allele and genotype frequencies between bee samples that differed in Nosema infestation was performed using the chi-square test. In the case of a small number of one of the comparison classes, the chi-square test with Yates' correction was used.
To assess the association of polymorphic variants of the microsatellite loci studied with Nosema infection in bees, the odds ratio (OR) and 95% confidence interval (95% CI) corresponding to the p-value was calculated [77]. The association of the genetic marker with the tested trait was determined by the OR value when the differences in the allele frequency between the compared groups reached the level of statistical significance, p < 0.05. If the OR > 1, the assumption about the association of the analyzed genetic variant (allele/genotype) with the studied pathology is considered (increased chance of developing disease). When the OR < 1, a protective role of the corresponding genetic variant (allele/genotype) is assumed.

Genetic Diversity of the A. m. mellifera Honey Bees in Siberia on the Microsatellite Loci
Honey bees with different degrees of Nosema infestation were genotyped using 23 microsatellite markers. When comparing the distribution of allele frequencies between three groups of the A. m. mellifera bees (Nosema-negative, Nosema-positive low and Nosemapositive high), several loci (AC117, A113, Ap243, A024, A007, Ap049 and SV185) were identified that are promising for further analysis. The parameters of genetic diversity of these loci, such as the frequencies of alleles and genotypes, expected heterozygosity and observed heterozygosity, are presented in Table 1. All microsatellite loci studied were polymorphic: the minimum number of alleles was found for locus A007 (3 alleles), and the maximum number of alleles was for loci A113 and Ap243 (8 alleles); the average number of alleles per locus was 5 (Table 1).
An assessment of the heterozygosity of most of the studied loci (except for loci A007 and SV185 in Nosema-positive high bees) revealed lower values of the observed heterozygosity (Ho) compared with the expected heterozygosity (He) ( Table 1). A statistically significant level of differences between the values of the observed and expected heterozygosity is shown for the following loci: A113 (t = 3.82, p < 0.001 and t = 3.16, p < 0.05), Ap243 (t = 4.35, p < 0.001 and t = 2.09, p < 0.05) in Nosema-positive low and high bees, respectively; SV185 (t = 2.77, p < 0.05) in Nosema-positive low bees; AC117 (t > 3.69, p < 0.001) in all infection categories of bees.
Thus, the analysis of the variability of 23 microsatellite loci in the A. m. mellifera honey bees allowed us to identify the most promising loci for searching for associations of DNA markers with the Nosema disease.

Comparative Characteristics of the Genetic Diversity of A. m. mellifera Bees from Different Infectious Categories
The heterogeneity of bee groups differing in the degree of Nosema infection was evaluated. To do this, statistically significant differences between the bee groups in allele frequencies of microsatellite loci were determined.
Despite the fact that no significant differences were found in the distribution of all alleles (in total) of loci A007, A024 and SV185, there existed differences for some alleles of these loci. For example, allele "269" of the SV185 locus determined the differences between groups of Nosema-uninfected and Nosema-positive high bees (χ 2 = 5.65, df = 1, p < 0.05).

Assessment of Associations of Genetic Markers with Nosema Infection/Resistance in the Dark Forest Bee A. m. mellifera
In order to search for alleles, possibly associated with Nosema disease in the dark forest bees (A. m. mellifera), the odds ratio, OR, was calculated ( Table 2).
For loci AC117, A113, Ap243, A024, A007, Ap049 and SV185, statistically significant differences in allele and/or genotype frequencies between the compared groups (infection categories) were shown. However, according to OR calculations, only some genetic variants of these loci showed associations with Nosema infestation. Based on the calculated OR, it can be concluded that alleles "177" of the AC117 locus, "263" of the Ap243 locus and "269" of the SV185 locus have protective properties (that is, they reduce the risk of Nosema infection). Table 2. Comparative analysis of the frequency of potentially significant genetic variants in the formation of nosemosis resistance in the dark forest bee A. m. mellifera. It is worth pointing out that the frequency of homo-and heterozygous genotypes with an allele "177" of the AC117 locus decreased in the sequence: uninfected bees (25.9% of Nosema-negative individuals) were weakly infected (14.1% of Nosema-positive low bees) and significantly infected (4.9% of Nosema-positive high individuals) ( Table 1). Between the extreme compared groups of bees (Nosema-negative and Nosema-positive high) differences in the frequency of individuals with genotypes having this allele reach the level of statistical significance (χ 2 = 16.61, df = 2, p < 0.01). For the Ap243 locus, the frequency of genotypes with the "263" allele in Nosema-positive high bees (11.1%) differed significantly from both uninfected bees (43.3%; χ 2 = 13.45, df = 6, p < 0.05) and Nosema-positive low individuals (43.5%; χ 2 = 18.86, df = 7, p < 0.01) ( Table 1). For the SV185 locus, there were no differences between the compared groups at the genotype level (p > 0.05). In addition, the results obtained suggest that the "218" allele of the A113 locus is associated with an increased risk of nosemosis in A. m. mellifera bees (no differences were recorded at the genotype level).

Discussion
The genetic diversity of the A. m. mellifera bees, differing in the degree of Nosema infection, was investigated, and a search for associations between genetic variants of microsatellite loci and Nosema disease was carried out. The present results show that some alleles of microsatellite loci are associated with nosemosis in the A. m. mellifera bees living in the Siberian region. For example, alleles "177" of the AC117 locus, "263" of the Ap243 locus and "269" of the SV185 locus reduce the risk of Nosema infection.
This study found associations of genetic markers with Nosema infection/disease resistance in honey bees and was an exploratory study. Therefore, it is necessary to understand whether the results obtained were random or reflected some general biological regularities. Possible solutions to this issue are expanding the bee samples and/or analyzing different bee subspecies and studying them in terms of the importance of microsatellite loci in determining resistance to nosemosis.
Thus, among the investigated microsatellite loci, the Ap243 locus located on chromosome 1 (group 1.1) was shared for two bee subspecies (A. m. mellifera and A. m. carpathica) and is probably associated with the incidence/resistance to nosemosis. However, in these subspecies, different alleles associated with Nosema disease were identified. For A. m. mellifera bees, allele "263" probably reduced the risk of infection with nosemosis (protective role), and statistically significant differences were also found at the genotype level in bees of different groups (Nosema-positive high bees differed from uninfected and Nosema-positive low bees). In A. m. carpathica bees, statistically significant differences were found for two alleles ("256" and "253") between uninfected bees and Nosema-infected individuals, and the allele "256" is likely to have a protective meaning, while the allele "253", on the contrary, is associated with a disease.
It should be noted that for two bee subspecies, various microsatellite loci and alleles associated with the Nosema incidence/resistance have been identified, which may be determined by the different resistance of bee subspecies to nosemosis and/or by different habitats of bees (geographic, natural, climatic and nutritional conditions) [58,[79][80][81][82]. Perhaps, in addition to the subspecies-specific features to Nosema resistance, the revealed differences can be determined by the structure of the chromosomal region where the QTL is located. For example, not this locus itself, but another, which is in linkage disequilibrium with it, may be involved in the determination of resistance to nosemosis, i.e., different alleles are in the same linkage group with a favorable variant in different bee subspecies.
The few hygienic behavior-related studies focused on finding quantitative trait loci controlling Varroa resistance hygienic behavior of honey bees have also identified different chromosomal regions related to hygienic behavior and reduced mite reproduction [14,[20][21][22]31]. QTL studies of Varroa resistance behavior have identified over 20 suggestive chromosome regions associated with linkage groups 2, 3, 4, 5, 6, 7, 9, 10, 13, 15, 16 and 22. For example, using RAPD markers, Lapidge et al. (2002) found seven suggestive QTLs controlling hygienic behavior of honey bees [22]. Using microsatellite loci, Oxley et al. (2010) identified three significant and three suggestive QTLs that influence a honey bee worker's propensity to engage in hygienic behavior [31] whereas Behrens et al. (2011) identified three QTLs controlling reduced mite reproduction in the Swiss Varroa mite tolerant honey bee lineage [20]. In the QTL study presented by Spötter et al., six SNPs showed significant genome-wide associations with hygienic behavior against Varroa at the genotype level [14].
It is worth pointing out that the presented studies did not identify the same (common) QTLs associated with hygienic behavior of honey bees.  [20,31]. Additionally, Tsuruda et al. (2012) also reported one major QTL on chromosome 9, but in a different region, for hygienic behavior against Varroa using a small-scale SNP-Chip [21]. Thus, the identified genomic regions related to hygienic behavior are different, which can be explained by different bee materials (freeze-killed brood, brood and worker bees), different research methods used and DNA markers analyzed (RAPD markers, microsatellite loci or SNP), different genetic maps (map based on RAPD or microsatellite markers). Finally, the use of different bee subspecies can significantly affect the research outcome. For example, in QTL studies presented by Oxley et al. (2010) and Spötter et al. (2016), QTLs were identified in the same chromosomes (chromosomes 2 and 5), but at different ends of the chromosome [14,31]. It is assumed that different bee material used in these studies contributed to the differences in the genomic regions identified [14]. On the one hand, a freeze-killed brood instead of brood that was artificially infested with Varroa was used in Oxley et al.'s research [31]. On the other hand, two different bee subspecies (Apis mellifera ligustica and Apis mellifera carnica) have been investigated [14].
A comparison between the QTLs involved in N. ceranae infection tolerance according to Huang et al.'s study [35] and our own trait-associated regions found no agreement. So, four QTLs located on chromosomes 3, 10, 6 and 14 were significantly associated with low Nosema spore load, explaining 20.4% of total spore load variance in the selected Nosema-resistant Danish honey bee strain. The significant QTL on chromosome 14 explains 7.7% of the total variance and may be responsible for the resistance to nosemosis in the selected Danish honey bee. A candidate gene Aubergine (Aub) within this QTL region was significantly more overexpressed in drones with a low spore load than in those with a high spore load [35].
From the data presented, in the dark forest bee, microsatellite loci Ap243 (chromosome 1), SV185 (chromosome 5) and AC117 (chromosome 12) are associated with resistance to nosemosis. It is interesting to note that in the chromosome region 1.1, the microsatellite locus Ap243 and the honey bee microRNA (ame-miR-2b) are located. As shown, host microRNAs respond to infection by the parasite N. ceranae [36]. In honey bees, 17 miRNAs were differentially expressed during N. ceranae infection, which may target over 400 predicted genes for ion binding, signaling, the nucleus, transmembrane transport and DNA binding. MicroRNA ame-miR-2b is particularly interesting because 11 out of 27 enzymes were significantly correlated with its expression level [36]. In addition, in the same region on chromosome 1, a suggestive QTL associated with the performance of Varroa sensitive hygiene was identified [21]. This locus contains more than candidate 30 genes [21], including puromycin-sensitive aminopeptidase involved in proteolytic events essential for cell growth and viability [83], selenoprotein F located in the endoplasmic reticulum and regulated by cell stress conditions [84], transcription and splicing factors and other genes. Since Microsporidian Nosema species are intracellular parasite [37], of particular interest is the cell wall integrity and stress response component 1. In Schizosaccharomyces pombe, the homologous gene wsc1 is responsible for such biological processes as a cell surface receptor signaling pathway, intracellular signal transduction and regulation of cell wall organization or biogenesis [85].
Despite the fact that numerous QTLs associated with the disease resistance of honey bee have been identified, it is assumed that the variations in this trait is controlled by a small number of loci. For example, a strong genetic component is involved in the control of hygienic behavior, although it may also be influenced by environmental factors to some extent [14]. Therefore, the identification of the gene variants that are responsible for disease resistance in honey bees, and then breeding resistant bee colonies is one promising approach in the fight against bee disease.

Conclusions
The present study examined the associations between several variants of microsatellite loci and Nosema disease. In the dark forest bee, genetic markers promising for the assessment of nosemosis resistance, such as the allele "177" of the locus AC117, the allele "263" of the locus Ap243, the allele "269" of the locus SV185 have been identified. At the same time, some issues, for example, differences in the spectrum of loci and/or alleles that determine resistance to nosemosis in different bee subspecies/breeds, remained unresolved. In this regard, additional research both for the same bee subspecies/breeds bred in other regions, and for other bee subspecies/breeds is needed. However, already at this stage, these markers can be used to predict the risk of developing nosemosis, but with the obligatory consideration of the specificity of diagnostic markers for different bee subspecies/breeds.