<figure><img src="image_1.png" /></figure>

<figure><img src="image_2.png" /></figure>

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

M. Guichard*† ID, B. Dainat* ID, S. Eynard‡§ ID, A. Vignal‡*,§ ID, B. Servin‡*,§ ID,
the Bestrong Consortium¶ and M. Neuditschko† ID

*Agroscope, Swiss Bee Research Centre, Schwarzenburgstrasse 161, Bern 3003, Switzerland. †Agroscope, Animal GenoPhenomics, Rte de la Tioleyre 4, Posieux 1725, Switzerland. ‡GenPhySE, INRA, INPT, INPENVT, Université de Toulouse, Castanet-Tolosan 31320, France. §UMT PrADE, Protection des Abeilles Dans l'Environnement, Avignon 84914, France. ¶Labogena, Domaine de Vilvert Bat 224 CS80009, Jouy-en-Josas CEDEX 78353, France.

## Summary

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.

**Keywords** ataxin-10, *Dscam*, genome-wide association study, honey bee, pool sequences, recapping, Wnt7

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. carica* (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

## Address for correspondence

M. Guichard, Agroscope, Swiss Bee Research Centre,
Schwarzenburgstrasse 161, 3003 Bern, Switzerland.
E-mail: matthieu.guichard@agroscope.admin.ch

Accepted for publication 10 October 2021
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).

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

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,

<figure><img src="image_3.png" /><figcaption>**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.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**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.</figcaption></figure>

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.

## Acknowledgements

Financial support for this study was provided by Bundesamt für Landwirtschaft BLW (Swiss Federal Office for Agriculture; grant no. 627000708) and Labogena and FranceAgriMer (Programme d'Investissements d'Avenir). We thank Lucie Genestout (Labogena) for coordinating the Beestrong project, Christina Sann and Yannick Poquet (Labogena) for organising the sampling in Savoie and Alice Revel (Agroscope) for her help in organising the sampling in Switzerland. We thank Rachid Mahla (Labogena) for DNA extraction and Olivier Bouchez (GeT-PlaGe) for sequencing. Sequencing was performed in collaboration with the GeT platform, Toulouse (France), a partner of the National Infrastructure France Genomique, with the support of the Commissariat aux Grands Investissements (ANR-10-INBS-0009). Sequence alignments and allele frequencies estimations were performed on the computer cluster of the Genotoul bioinformatics platform Toulouse Midi-Pyrénées (www.bioinfo.genotoul.fr). We are grateful to all the experts for colony evaluation and to all the beekeepers for supporting the sample collection. This study was performed in collaboration with the Beestrong Consortium. Beestrong Consortium: Sonia Eynard, Alain Vignal, Bertrand Servin (INRAE, GenPhyse, 31326 Castanet-Tolosan, France); Lucie Genestout, Rachid Mahla, Christina Sann, Yannick Poquet, Yves Agez (LABOGENA DNA, Laboratoire d'Analyses Génétiques pour les Espèces Animales, Domaine de Vilvert, 78353 Jouy-en-Josas, France); François Guillaume (Evolution XY, rue Eric Tabarly, 35538 Noyal-Sur-Vilaine, France); Benjamin Basso, Yves Le Conte, Fanny Mondet (INRAE, Abeilles et Environnement, 84914 Avignon, France & UMT PrADE, 84914 Avignon, France); Benjamin Basso, Axel Decourtye (ITSAP, 84914 Avignon, France); Florence Phocas, Jean-Pierre Bidanel (Université Paris-Saclay, INRAE, AgroParisTech, GABI, 78350, Jouy-en-Josas, France); Olivier Bouchez (GeT-PlaGe, Genotoul, INRAE Auzeville, 31320, Castanet-Tolosan, France). Open Access Funding provided by Agroscope.

## Conflict of interest

The authors declare no conflicts of interest.

## Data availability statement

The data that support the findings of this study remain the property of Agroscope (Swiss samples) and the Beestrong Consortium (French samples). However, the data are available from the authors upon reasonable request.

## References

Arechavaleta-Velasco M.E., Alcala-Escamilla K., Robles-Rios C., Tsuruda J.M. & Hunt G.J. (2012) Fine-scale linkage mapping reveals a small set of candidate genes influencing honey bee grooming behavior in response to *Varroa mites*. *PLoS One* 7, 1–7.

Büchler R., Costa C., Mondet F., Kezic N. & Kovacic M. (2017) Screening for low Varroa mite reproduction (SMR) and recapping in European honey bees. Available at: https://www.beebreeding.net/index.php/2017/09/01/new-smr-protocol/. Accessed 27 February 2019.

Dearden P.K., Wilson M.J., Sablan L. et al. (2006) Patterns of conservation and change in honey bee developmental genes. *Genome Research* 16, 1376–84.

Diao Q., Sun L., Zheng H. et al. (2018) Genomic and transcriptomic analysis of the Asian honeybee *Apis cerana* provides novel insights into honeybee biology. *Scientific Reports* 8, 822.

Guichard M., Dainat B., Eynard S., Vignal A., Servin B. The Beestrong Consortium & Neuditschko, M. (2021) Identification of quantitative trait loci associated with calmness and gentleness in honey bees using whole-genome sequences. *Animal Genetics* 52, 472–81.

Guichard M., Dietemann V., Neuditschko M. et al. (2020) Advances and perspectives in selecting resistance traits against the parasitic mite *Varroa destructor* in honey bees. *Genetics Selection Evolution* 52, 71.

Khongphinitbunjong K., de Guzman L.I., Tarver M.R., Rinderer T.E., Chen Y. & Chantawannakul P. (2015) Differential viral levels and immune gene expression in three stocks of *Apis mellifera* induced by different numbers of *Varroa destructor*. *Journal of Insect Physiology* 72, 28–34.

Le Conte Y., Alaux C., Martin J.F., Harbo J.R., Harris J.W., Dantec C., Severac D., Cros-Arteil S. & Navajas M. (2011) Social immunity in honeybees (*Apis mellifera*): transcriptome analysis of varroa-hygienic behaviour. *Insect Molecular Biology* 20, 399–408.

Martin S.J., Hawkins G., Brettell L., Reece N., Correia-Oliveira M.E. & Allsopp M.H. (2019) *Varroa destructor* reproduction and cell recapping in mite-resistant *Apis mellifera* populations. *Apidologie* 51, 369–81.

März P., Probst A., Lang S., Schwager M., Rose S., Otten J.U. & Özbek S. (2004) Ataxin-10, the spinocerebellar ataxia type 10 neurodegenerative disorder protein, is essential for survival of cerebellar neurons. *Journal of Biological Chemistry* 279, 35542–50.

Mondet F., Beaurepaire A., McAfee A., Locke B., Alaux C., Blanchard S., Danka R. & Le Conte Y. (2020) Honey bee survival
mechanisms against the parasite *Varroa destructor*: a systematic
review of phenotypic and genomic research efforts. *International
Journal for Parasitology* **50**, 433–47.

Navajas M., Migeon A., Alaux C., Martin-Magniette M.L., Robinson G.E., Evans J.D., Cros-Arteil S., Crauser D. & Le Conte Y. (2008) Differential gene expression of the honey bee *Apis mellifera* associated with *Varroa destructor* infection. *BMC Genomics* **9**, 1–11.

Oddie M.A.Y., Büchler R., Dahle B., Kovacic M., Le Conte Y., Locke B., De Miranda J., Mondet F. & Neumann P. (2018) Rapid parallel evolution overcomes global honey bee parasite. *Scientific Reports* **8**, 7704.

Purcell, S., Neale, B. & Todd-Brown, K. (2007) PLINK: a toolset for whole-genome association and population-based linkage analysis. *American Journal of Human Genetics*, **81**, 559–75.

Spötter A., Gupta P., Mayer M., Reinsch N. & Bienefeld K. (2016) Genome-wide association study of a *Varroa*-specific defense behavior in honeybees (*Apis mellifera*). *Journal of Heredity* **107**, 220–7.

Traynor K.S., Mondet F., de Miranda J.R., Techer M., Kowallik V., Oddie M.A.Y., Chantawannakul P. & McAfee A. (2020) *Varroa destructor*: a complex parasite, crippling honey bees worldwide. *Trends in Parasitology* **36**, 592–606.

Turner, S. D. (2014) qqman: an R package for visualizing GWAS results using Q-Q and manhattan plots. *BioRxiv*.

Wallberg, A., Bunikis, I., Pettersson, O. V., Mosbech, M. B., Childers, A. K., Evans, J. D., Mikheyev, A. S., Robertson, H. M., Robinson, G. E. & Webster, M. T. (2019) A hybrid de novo genome assembly of the honeybee, Apis mellifera, with chromosome-length scaffolds. *BMC Genomics*, **20**, 275.