Two quantitative trait loci are associated with recapping of Varroa destructor-infested brood cells in Apis mellifera mellifera.

Recapping of Varroa destructor-infested brood cells is a trait that has recently attracted interest in honey bee breeding to select mite-resistant Apis mellifera colonies. To investigate the genetic architecture of this trait, we evaluated a sample of A. mellifera mellifera colonies (N = 155) from Switzerland and France and performed a genome-wide association study, using a pool of 500 workers per colony for next-generation sequencing. The results revealed that two QTL were significantly (P < 0.05) associated with recapping of V. destructor-infested brood cells. The best-associated QTL is located on chromosome 5 in a region previously found to be associated with grooming behaviour, a resistance trait against V. destructor, in A. mellifera and Apis cerana. The second best-associated QTL is located on chromosome 4 in an intron of the Dscam gene, which is involved in neuronal wiring. Previous research demonstrated that genes involved in neuronal wiring are associated with recapping and varroa sensitive hygiene. Therefore, our study confirms the role of a gene region on chromosome 5 in social immunity and simultaneously provides novel insights into genetic interactions between common mite resistance traits in honey bees.

The invasive parasitic mite Varroa destructor remains a major threat to the global survival of the honey bee Apis mellifera (Traynor et al. 2020). Thus far, various resistance mechanisms have been investigated with the aim of selecting V. destructor‐resistant A. mellifera colonies (Guichard et al. 2020; Mondet et al. 2020). Currently, recapping, a trait observed in several natural A. mellifera populations (Oddie et al. 2018; Martin et al. 2019), is increasingly gaining the attention of scientists and beekeepers, with research suggesting this trait could provide resistance against V. destructor. Worker bees expressing recapping behaviour open and then re‐seal brood cells, which probably disturbs the reproduction cycle of V. destructor mites (Oddie et al. 2018). Compared to varroa sensitive hygiene (VSH), where workers remove the infested brood, recapping does not cause brood destruction (Oddie et al. 2018), which could favour colony survival.

The genetic background of recapping is yet not well understood. In this study, to investigate the genetic architecture of this trait, we derived pooled sequence information of 155 A. m. mellifera colonies, originating from a Swiss selection programme (referred to as SL_CH) and two conservation areas in Switzerland and France (CS_CH and CS_FR respectively), as well as 28 A. m. carnica (CAR) colonies, from a recently described dataset (Guichard et al. 2021). We applied the same quality control criteria as used in the study by (Guichard et al. 2021), which resulted in 1 355 136 genome‐wide SNPs for subsequent analyses.

From each colony, a single worker brood sample was collected during the summer. Only cells containing pupae at least at the purple eye stage (7 days post‐capping) were included. The status of the cell cap, either untouched or recapped, was assessed based on a standard protocol (Büchler et al. 2017). Following the removal of the pupa, the presence or absence of at least a single founder mite in the cell was investigated. Phenotype evaluation of a given colony was terminated after 35 single infested cells were identified in the corresponding brood sample. The status of the cell cap and presence of mites were combined to calculate the percentage of infested and recapped cells in each colony.

To identify QTL involved in recapping, we performed a genome‐wide association study (GWAS) on sequence‐derived SNP dosage data on 155 A. m. mellifera colonies using a linear regression model implemented in plink (Purcell et al. 2007). We adjusted the model for covariates capturing population stratification and significant effects on the trait, following the method described by Guichard et al. (2021). Significant associated SNPs were determined based on a 5% genome‐wide Bonferroni‐corrected threshold. The results of the GWAS were visualised using a Manhattan plot and quantile–quantile plot with the r package qqman (Turner 2014). We also explored the effect and allele frequency of the best associated QTL within each sub‐population, including the CAR colonies.

After verification of data normality, significant sub‐population effects on allele frequencies were identified by an analysis of variance, followed by a Tukey multiple comparison of means, with a 95% confidence interval. Genes within the identified QTL regions were determined using the NCBI Genome Data Viewer (https://www.ncbi.nlm.nih.gov/genome/gdv/browser/genome/?id=GCF_003254395.2) and the reference genome assembly Amel_HAv3.1 (Wallberg et al. 2019).

Figure 1 summarises the observed recapping values in the different sampled sub‐populations. It shows a significantly higher recapping rate in the CS_FR sub‐population as compared with that in the other two A. m. mellifera subpopulations (CS_CH and SL_CH). In contrast, the difference in recapping values between the CS_FR and CAR colonies was not significant (P < 0.05, Tukey’s multiple comparison of means).

Figure 1

Jitterplot of uncorrected observations for the percentage of infested cells that are recapped, according to subgroup. Black rhombus corresponds to mean value, bars correspond to standard deviation. Different letters indicate significant (P < 0.05) differences between groups following a Tukey multiple comparison of means with a 95% confidence interval.

The GWAS on recapping was adjusted for population stratification using two principal components (PCs), which accounted for 99% of the total variance (PC1 = 98%, PC2 = 1%), and two covariates showing a significant effect on the trait (year and apiary). After adjustment of the data, recapping was significantly associated with two QTL on chromosomes 4 and 5 (Fig. 2a).

Figure 2

Genome‐wide association study. (a) Manhattan plot and quantile–quantile plots for percentage of infested cells recapped (N = 155 MEL colonies, outliers removed). The red line is the threshold for SNPs having a significant (P < 0.05) effect on phenotype. Two SNPs have a highly significant effect. The best SNP located at 805,163 base pairs on chromosome 5 does not correspond to a gene. The other significant SNP located at 11,852,817 base pairs on chromosome 4, is situated in the Dscam gene. (b) Percentage of infested cells recapped (uncorrected phenotype) according to percentage of A allele of the best SNP (chromosome 5) in the three Apis mellifera mellifera subgroups (selected SL_CH, conserved CS_CH and CS_FR). (c) Mean percentage of A allele for best SNP and associated standard deviation in each subgroup. Different letters indicate significant (P < 0.05) differences between groups following a Tukey multiple comparison of means with a 95% confidence interval.

The best‐associated QTL on chromosome 5 (805 163 bp) is not embedded in a gene region. The two nearest genes, LOC726806 and LOC411919, are located 5 kb downstream and 10 kb upstream of the QTL respectively. Dividing the observed recapping rate of the colonies into two groups according to the allelic frequency of the best‐associated QTL showed that colonies segregating the A allele at high frequency (>50%) expressed a high recapping rate, whereas in colonies carrying the A allele at low frequency (<50%), the recapping rate was relatively low (Fig. 2b). The associated A allele was highly segregated within the CAR and CS_FR colonies. In contrast, in the SL_CH and CS_CH colonies, the frequency of the allele was below 50% and 25% respectively (Fig. 2c).

LOC726806 is a gene coding for the protein ataxin‐10, which is involved in the functioning of the nervous system (März et al. 2004). Genes coding for ataxin‐10 were reported to be associated with grooming behaviour in both the western honey bee, A. mellifera (Arechavaleta‐Velasco et al. 2012), and the eastern honey bee, Apis  cerana (Diao et al. 2018). The second nearest gene, LOC411919, codes for the Wnt7 protein, which was shown to be involved in cell signalling pathways in A. mellifera (Dearden et al. 2006). Therefore, our results confirm that the previously identified gene region on chromosome 5 is associated with social immunity in A. mellifera.

The second best‐associated QTL identified in the present study (chromosome 4, 11 852 817 bp) is located in an intron of the Dscam gene. Previous studies demonstrated that this gene is downregulated in naturally surviving A. mellifera colonies and in A. mellifera colonies selected for VSH (Navajas et al. 2008; Le Conte et al. 2011). Dscam is involved in neuronal development and causes a different neuronal wiring in the brain of VSH bees (Le Conte et al. 2011). A recent GWAS reported that the cdk5alpha gene, located on chromosome 3, and also involved in neuronal wiring, is associated with the detection and uncapping of V. destructor‐infested cells (Spötter et al. 2016), suggesting that workers with specific neuronal abilities could better detect mites present in the brood, and target them by recapping or VSH. Interestingly, Dscam regulation in A. mellifera pupae is affected by the presence or absence of a parasitising mite. A study that compared different stocks in North America found out that Dscam expression was downregulated in mite‐infested pupae from an Italian line, as well as in a line selected for VSH, although at a lower level (Khongphinitbunjong et al. 2015). In contrast, in a resistant Russian A. mellifera population, no association was found between infestation and downregulation of Dscam (Khongphinitbunjong et al. 2015). The potential effects of Dscam regulation and infestation status at the pupal stage on the expression of the Dscam gene in adult bees remain unknown.

In this study, we identified two QTL associated with recapping of infested brood using whole‐genome sequences of 155 A. m. mellifera colonies. It should be noted that for two additional investigated mite‐related traits, including the infestation level of the worker brood and the infestation level of adult workers, no QTL were detected based on our dataset (data not shown). The two identified candidate genes for recapping are involved in the nervous system of A. mellifera and associations with grooming and VSH, respectively, have already been established. Our study provides additional evidence for the presence on chromosome 5 of a major QTL involved in social immunity. However, further data and research are needed to better understand the interrelationship and genetic architecture of VSH, recapping and grooming. In the meantime, we recommend using the best‐associated QTL identified in the present study for marker‐assisted selection to improve the selection of V. destructor‐resistant A. mellifera colonies.

Conflict of interest

The authors declare no conflicts of interest.