Heritability estimates of the novel trait ‘suppressed in ovo virus infection’ in honey bees (Apis mellifera)

Honey bees are under pressure due to abnormal high colony death rates, especially during the winter. The infestation by the Varroa destructor mite and the viruses that this ectoparasite transmits are generally considered as the bees’ most important biological threats. Almost all efforts to remedy this dual infection have so far focused on the control of the Varroa mite alone and not on the viruses it transmits. In the present study, the sanitary control of breeding queens was conducted on eggs taken from drone brood for 4 consecutive years (2015–2018). The screening was performed on the sideline of an ongoing breeding program, which allowed us to estimate the heritabilities of the virus status of the eggs. We used the term ‘suppressed in ovo virus infection’ (SOV) for this novel trait and found moderate heritabilities for the presence of several viruses simultaneously and for the presence of single viral species. Colonies that expressed the SOV trait seemed to be more resilient to virus infections as a whole with fewer and less severe Deformed wing virus infections in most developmental stages, especially in the male caste. The implementation of this novel trait into breeding programs is recommended.

The economic value of the honey bee can be almost entirely attributed to its pollination services for agricultural crops 1 . However, honey bees are under pressure due to abnormal high colony death rates, especially during the winter 2 . There are several possible explanations for the increased mortality, but the infestation by the Varroa destructor mite and the viruses that this ectoparasite transmits are generally considered the most important biotic threats 3,4 . The mite was found to be an efficient vector of viruses even to the extent that initially the clinical symptoms caused by the Deformed wing virus (DWV) were wrongly attributed to the mite 5 . The Varroa mite and DWV interact in many different ways, resulting in an increase in DWV virulence 6,7 and colony mortality [8][9][10][11][12] . These effects are partially due to the mutualistic symbiotic relationship between both, in which the mite provides transmission of the virus when it feeds on the bee, whereas the virus undermines the immunity of the honey bee by interfering with NF-κB signaling, possibly facilitating the mite's trophic activity 13 . Almost all efforts to remedy this dual infection have so far focused on the control of the Varroa mite alone. Indeed, several mitecontrol strategies have been set up in beekeeping practice that rely on acaricide medication or other treatments, biotechnical apicultural methods and selection for Varroa resistance 14 . The control of the viruses transmitted by the mite is therefore achieved almost exclusively indirectly, by limiting the severity of the mite infestation. One exception to this is the treatment with double-stranded RNA to provoke a targeted antiviral immune response based on RNA interference, but this is rarely used in beekeeping practice 15,16 .
The transmission of viruses between bees by the vectoring Varroa mite is described as horizontal transmission, as it is from one individual to another from the same generation by either direct or indirect contact. Virus transmission through bodily contact, trophallaxis, feeding and common flower visits are also categorized under horizontal transmission. On the other hand, vertical transmission refers to the passing of an infectious agent from parent to offspring via eggs and semen. It has been demonstrated that DWV vertical transmission occurs predominantly by virus particles adhering to the surface of the egg (transovum) rather than intracellularly 17 18 , which seems to be an important cause of colony collapse 23,24 . The queen's ovaries can become heavily infected 18,25 leading to typical pathological reactions with discoloration 26 . Because virus infection of a queen can influence her performances and represents an important source of infection for her progeny by vertical transmission, we started in 2012 a sanitary control of the breeding queens from Flemish beekeepers in the North of Belgium following a non-destructive approach. Freshly laid eggs taken from worker brood were examined for the commonly occurring bee viruses and led to the observation that vertical transmission occurs rather frequently 27 . Here we present the continuation of this study in which the sanitary control of the queens was conducted on eggs taken from drone brood for 4 consecutive years (2015-2018), before finally estimating the heritabilities of the virus status of the eggs. The term 'suppressed in ovo virus infection' (SOV) was given to this novel trait. We hypothesized that queens that were capable of clearing a visceral virus infection and consequently deposited virus-free eggs, would pass this predisposition of protective anti-viral immune responses on to their descendants. Drone brood eggs were our target samples as they are unfertilized and thus carry only the alleles of the breeding queen.

Results
Descriptive statistics. Fifty four Flemish queen breeders volunteered drone eggs to determine the virus status of the queens (from 2015 to 2018). The dataset contained 625 samples of 560 Carnica queens, 25 Buckfast queens and 23 queens of unknown subspecies (see Supplementary Table S1 online). Sixteen Carnica queens provided two samples and one queen three samples. The samples were analyzed for the presence of DWV, Acute bee paralysis virus (ABPV), Black queen cell virus (BQCV) and Sacbrood virus (SBV) and were scored 0 or 1 (0 = absence; 1 = presence) for each of the viruses and for total virus status (TVS). At the start (2015), the number of breeding queens involved in the sanitary screening was rather low (see Supplementary Fig. S1A online). This changed in 2016, when for the first time close to 200 queens were included. In the beginning, the vast majority of tested queens were bred from mother queens with an unknown virus status (UQ subgroup). But, the proportion of descendants of virus-negative queens (DV-Q subgroup) grew steadily and represented in 2018 almost 40% of the tested queens (see Supplementary Fig. S1A online).
The percentage of TVS-positive samples fluctuated slightly over the successive years (see Supplementary  Fig. S1B online). And when the entire period 2015-2018 was taken in consideration, the percentage of viruspositive samples (any virus) differed slightly between the different subgroups ( Fig. 1). However, for the virus species separately, some remarkable differences could be observed. The occurrence of DWV and SBV in eggs was considerably lower in the DV-Q subgroup as compared to the DV + Q (descendants of virus-positive queens) and UQ subgroups. That was not the case for ABPV and BQCV (Fig. 1). And we found the most ABPV, DWV and SBV, in the DV-Q, UQ and DV + Q subgroups, respectively.
When we look at the successive years, the DV-Q subgroup mostly scored better (lower) for DWV (except in 2016) and SBV and the occurrence of these viral species fell to its lowest level in 2018 with only 6.1 and 4.5% of the samples being positive, respectively (Fig. 2). On the contrary, no clear pattern could be found for ABPV and BQCV, though in 2018 the number of BQCV-positive samples was only half in DV-Q compared to DV + Q. Although these data have no burden of proof, they suggest that the DWV and SBV virus status of the daughter queens (as determined by the virus status of the drone eggs) can at least partially be explained by genetic factors. Heritabilities. The estimation of heritabilities was restricted to drone eggs of 560 Carnica queens. The majority of these queens with their colonies were part of an ongoing international breeding program (Beebreed) with some 10,000 colonies annually. Table 1 summarizes the heritabilities estimated for the status of individual viruses and for TVS. We used a threshold model, assuming that the phenotypes had a binomial distribution, while the underlying genetic liability followed a normal distribution. We used the threshold model because the phenotypes were scored 0 and 1 and the trait incidence was low. The heritability for SOV/TVS was estimated to be 0.25. The estimated heritabilities for three individual viruses were of the same order of size. The estimate for SOV/BQCV was larger but due to the limited size of the dataset the standard errors of all estimates were high. Nevertheless it may be concluded that the different forms of SOV were moderately heritable.
Environmental similarity of dams and daughters. The statistical model used to estimate the heritability also included the virus score of the dam. This is a fixed effect and represents the non-genetic (environmental) effect of the virus score of the dam on that of the daughter, as expressed in their respective drone eggs.  www.nature.com/scientificreports/ The significance levels (p-values) for the fixed effect for queen's virus loads were presented in Table 1. Various routes of this environmental similarity of dams and daughters might be envisaged like the direct transfer of virus from dams to daughters through eggs or various developmental stages. Also contamination in the internal environment as wax or royal jelly might be possible, or in the external environment as through the beekeeper or foraging. The statistical analysis did not confirm the relevance of those routes.
correlations. Ideally one would like to estimate genetic correlations between the different virus traits. In a dataset of this size that is not useful. Alternatively, phenotypic correlations were computed between the traits, albeit that trait variables were adjusted for the effects of the test location classes to which they belong. Table 2 summarizes these correlations between the traits. Generally speaking the correlations between individual viruses were low while SOV/BQCV and SOV/DWV correlated higher to SOV/TVS than the other two viruses.
Virus load of all developmental stages. We subsequently determined the impact of the SOV trait on the DWV and SBV infection levels of four developmental stages of drones and workers of 4 SOV colonies (SOV group) in comparison to 4 control colonies (control group). The colonies were kept at 4 apiaries, with 1 SOV colony and 1 control colony each. Due to the extremely low SBV infection levels, these results were not informative (see Supplementary Table S2 online). The DWV infection levels, on the other hand, allowed comparison between the two groups (see Supplementary Table S3  Together with the equal distribution of the SOV group and the control group between the apiaries we could assure that the trait rather than the environment or the circumstances at the apiary was here the determining factor.

Discussion
The present study suggests for the first time that the honey bee's potential to control virus infections is heritable.
There are some observations of isolated Varroa-resistant honey bee populations with lowered viral loads or diminished clinical signs 28 , but none of these studies demonstrated a genetic cause. Our findings are based on the sanitary control of breeding queens conducted on the sideline of an ongoing breeding program. The estimated heritabilities were moderate and were all about 0.25. For comparison, traditional honey bee traits like honey yield, gentleness and calmness, but also those related to Varroa resistance/tolerance, have heritabilities which vary from 0.17 to 0.91 for the effect of workers and from 0.10-0.70 for the effect of queens [29][30][31][32][33] . Indeed, several traditional traits proved to be affected by the genotypes of workers and queens as well, while the novel trait described in this study obviously solely is affected by the queen. The fact that the estimated heritabilities of SOV described in the present study are of the same size of order justifies its implementation in breeding programs. www.nature.com/scientificreports/ Given the fact that testing for the SOV trait demands laboratory analyses and the instability of the viral RNA genome under ambient temperature, the implementation of this trait requires to set in place a cold chain for the egg samples. This could be part of a concerted action of the responsible organisation or (beekeepers) association. We hypothesized that the virus status of the egg is a reflection of the health status-and by extension of the immune potential-of the queen, though we never verified this. As this study was interwoven in a running breeding program, in which breeders tested the offspring of their most precious queens, this was impossible to do. However, we did verify what the SOV trait means in terms of virus infections in the different developmental stages and castes of the colony. And the outcome was very promising: it seems that the honey bee egg holds the key of resilience to virus infection. The SOV trait provides no sterile immunity against viruses as they are still present, but some yet unknown mechanism seems to avoid that later exposures to-for instance-DWV will end with harmful infection levels.
The immunity of the queen against viruses certainly has also a social and behavioral component, but once a queen has acquired a virus infection, the propagation of the virus and the disease progression happens as in worker bees 34 . Bees like other insects have developed an antiviral immunity that relies on RNA interference (RNAi) among other pathways 35 . In an earlier study we found that the expression level of key components of the RNAi machinery are downregulated in heavily virus infected colonies 36 . We then suggested that viral RNAi www.nature.com/scientificreports/ suppressors mute the anti-viral immunity, a process that previously was described for the bees' Israeli acute paralysis virus 37 . The differences in expression profiles could, however, also be explained as a functional RNAi response that keeps viruses under control in the subgroup with the lowest virus load. No such responses have been studied in honey bee queens yet, but it seems plausible that RNAi plays a role in the expression of the SOV trait. In Caenorhabditis elegans it has been demonstrated that the antiviral RNAi response not only inhibits vertical transmission, but also promotes transgenerational inheritance of antiviral immunity. Such heritable antiviral immunity seems to be mediated by small RNAs 38 . If similar mechanisms would occur in honey bees, the predisposition of an effective antiviral immunity of the queens could be inherited to the worker bees and protect the colony from virus related pathologies, as suggested by our in depth study of the virus load of four SOV colonies and their controls. Given the fact that viruses play a crucial role in the collapse of honey bee colonies, the implications could be tremendous. Moreover, if worker bees are capable of clearing a virus infection the overall virus contamination level of the colony will drop gradually, a process from which even bees of neighboring colonies that lack the protective genetic profile will benefit. Our data are suggesting this is happening as also in the DV + Q subgroup both the DWV and SBV infections dropped considerably over the last 2 years.
In conclusion, we discovered a new trait that renders honey bee eggs free of virus infection. We used the term 'suppressed in ovo virus infection' and demonstrated that the trait is heritable through the genotype of the queen, and not that of the worker, and can be expressed against several viruses simultaneously or against each of the viruses individually. The estimated heritability seems to be moderate with a value of about 0.25. The trait has a beneficial effect on the virus load of the colony as a whole with fewer and less severe DWV infections and its implementation into breeding programs is recommended.

Sample collection. Flemish queen breeders provided between 1 and 45 samples each (see Supplementary
Table S1 online). From each selected colony, 10 eggs from drone cells were pooled in a pre-labeled 1.5 ml tube. After collection the eggs were immediately stored at − 20 °C and transported to the laboratory maintaining the cold chain where they were stored at − 80 °C until they were processed for virus status determination. All samples were collected in June (from 2015 to 2018).
At the end of the study (between June, 20 and July, 15, 2019) and with the queens' sanitary reports in hands, we selected 4 apiaries for an in depth study of the viral infection levels across different developmental stages and castes of the honey bee. The apiaries were managed by experienced beekeepers following standard beekeeping practices, including Varroa control. In each apiary one SOV colony and one control colony was selected based on the sanitary report of the queen in the previous year (2018) as determined by RT-PCR (see below) on a pooled egg sample. Of each colony 10 freshly laid eggs (pooled sample), 10 stretched larvae, 10 red-eyed pupae and 10 adults were collected for each caste (drones and workers). Where possible guard worker bees were selected in order to minimize the risk of sampling bees drifted from other colonies. Three drone larva samples could not be collected as one queen showed an irregular drone laying pattern and eventually ceased laying drone brood around the beginning of the sampling period. Placement of a drone brood frame in the hive to evoke the onset drone brood laying yielded insufficient success. Samples were kept in a cold chain as described here above.

Virus analysis by RT-PCR. The eggs were examined by uniplex RT-PCR assays for the presence of viruses
of the ABPV complex, BQCV, DWV and SBV. We used honey bee β-actin as a control gene. All PCR reaction mixtures contained: 2 μM of each primer (see Supplementary Table S4) . Each plate included a no target control (NTC). Viral loads of SBV and DWV were quantified using absolute quantification methods based on a standard curve obtained through a eightfold 5 × serial dilution of viral plasmid loads that ranged between 10 4 -10 10 copies/µl. All data were analyzed using CFX Manager 3.1 software (Bio-Rad). Baseline correction and threshold setting were performed using the automatic calculation offered by the same software. Maximum accepted quantification cycle (Ct) difference between replicates was set to one Ct. The successful amplification of the reference gene β-actin was used to confirm the integrity of samples throughout the entire procedure.
Quantitative genetic analysis. The purpose of the quantitative genetic analysis was to estimate heritabilities, to look at the environmental similarity of virus scores of queens and their dams, and to estimate correlations between the scores of different viruses. Because for this analysis the genetic relationships between queens are required the analysis was restricted to the records of 560 Carnica queens. This included 176 queens without a pedigree, but of the 384 other queens the dams were known and in 58% of the cases also the sires. These sires were drone-producing queens on mating islands (54% of the cases) or used to produce drones for instrumental insemination. The pedigrees were traced back to 2005 using historical information in the data base www.beebr eed.eu. Beebreed is an international data base that annually grows with some 10,000 queens with pedigrees and breeding values for traits like honey yield, gentleness and varroa-resistance. The extended pedigree was used to derive the additive genetic relationships matrix between individuals in the pedigree to allow estimation of heritabilities 39 . The linear animal model describing the observations was as in Eq. (1).

In this model
• y ik stands for the observation (SOV/ABPV, SOV/BQCV, SOV/DWV, SOV/SBV or SOV/TVS); • µ stands for the overall mean; • q i is a random effect representing the breeding value of the ith queen producing the drone eggs, giving an estimate for the additive genetic variance component σ 2 q ; • vd j is a fixed effect for the jth virus-score class of the dam. There are three of those classes: unknown, no virus present and virus present; • t k is a fixed effect for the kth test location in which an apiary in another year is considered to be a different test location. There were 101 test locations with 1 to 18 records each; • e ijk is a random residual term, giving a variance component for residual environmental effect σ 2 e .
The presence of both q i and vd j in the model enables to distinguish between genetic causes that the virus scores of queens resemble those of their dams on the one hand and non-genetic causes (environmental similarity) on the other hand. In theory, a model not including vd j might lead to overestimation of σ 2 q and consequently of the heritability.
Although in the model all effects were estimated simultaneously, one might say that all queens with observations contributed to estimate the effect of test location, that cases where queen and dam have an observation contributed to the estimation of environmental similarity and that the genetic relationships between queens allowed to estimate σ 2 q . The statistical analyses were carried out using ASReml 40 . We analysed the models assuming a binomial distribution of the phenotypes (Table 3) with an underlying normally-distributed genetic liability linked to the observed scale by a logit function. A binomial distribution seemed a logical choice as the scores were 0 or 1, and the low incidence makes a linear model not suited.
(1) y ijk = µ + q i + vd j + t k + e ijk where variance components σ 2 q were the result of the analyses and the implicit residual variance σ 2 e for the underlying logistic distribution was π 2 3 ≈ 3.3 for the logit link function 40 . correlations. The dataset was too small to estimate genetic correlations. Conceptually for each observation the deviation was computed from the estimate of the fixed effect of test location and pairwise Pearson correlations were computed between those deviations. The correlations and their standard errors were estimated using ASReml 40 with a model that only included the overall mean, the effect of test location and the random residual effect.

Statistics on qRT-PCR data.
To improve data compliance with statistical tests and data visualization, viral loads of all samples were log10 transformed. Data analysis and visualization were performed using 'R' version 3.6.1. The level of 10 9 DWV genomic copies per bee was used as a threshold to distinguish 'severely infected bees' . This threshold is reported in multiple studies as the transition between covert to overt infection [41][42][43] . For analyses of the differences in the number of DWV infections and severe DWV infections the Chi-squared test or the Fisher's exact test was used depending on the sample size. Differences in infection loads between groups were analyzed by one-way ANOVA. A possible disproportional effect of individual hives on the number of infections was analyzed using the Mann-Whitney test.

Data availability
Rough and trimmed data of honey bee breeding for suppressed in ovo virus infection, together with data on the kinship of the tested queens (pedigree) generated during the current study are available in the figshare repository: https ://doi.org/10.6084/m9.figsh are.81709 25.
Received: 10 February 2020; Accepted: 10 August 2020 (2) h 2 = σ 2 q σ 2 q + σ 2 e Table 3. Distribution of the data of the 560 Carnica queens. There were 16 queens that provided two samples and one queen with three samples. The virus score for those cases was the average of those of the different samples. The scores for TVS were the sum of those of the individual viruses. The results in Tables 1 and 2 resulted from analyses in which the presence of viruses was scored 0 or 1, where all scores larger than 0 are put to 1.