Discussion

We demonstrated that colony genotypes derived from pooled honey bee workers can successfully be applied to ascertain high-resolution population structures, including the computation of FROHand the detection of breed-specific selection signatures. To date, mostly drone genomes were used to assess the genetic diversity of honey bees, as their haploid nature facilitates cost-efficient whole-genome sequencing (Parejo et al., 2016). However, based on the haploid data structure, which exhibits systematic homozygosity, it is likely to overestimate genetic relationships and subsequently inbreeding, compared to other livestock species (Wragg et al., 2016). Another disadvantage of honey bee drones is that they only explain part of the genetic diversity, as multiple patrilines are involved in the formation of honey bee colonies (Estoup, Solignac, & Cornuet, 1994; Neumann, Moritz, & van Praagh, 1999; Tarpy, Nielsen, & Nielsen, 2004).
The applied population structure analyses clearly differentiated MEL from CAR, whilst the FST = 0.26 between the two subspecies was lower compared to our previous findings using drone genomes (FST = 0.36) (Parejo et al., 2016). This reduced genetic difference might be the result of highly admixed MEL colonies, which simultaneously highlights the challenges to conserve native honey bees due to random mating. Therefore, the high-resolution population network illustrates that a successful honey bee conservation programme requires an appropriate management tool including a legal framework, a suitable geographical isolated location, and ancestry informative marker testing, like the conservation strategy of MEL in the Canton Glarus, the only area without admixed colonies (CS_CH and P1). However, in our view the strong gene flow between the two subpopulations can have a negative impact on the in-situ conservation as foreign genetic variants might be introduced to the CS_CH gene pool by the mating station.
Compared to CS_CH, the origin of CS_FR colonies has only been sporadically assessed in the past based on wing vein measurements, which simultaneously explains the highly observed diversity of the colonies, whereas four colonies showed a high genetic relatedness with CS_CH. The population structure of CS_FR and the random mating events observed in some SL_CH patrilines indicate that current applied conservation strategies including the geographical locations are not suitable for in-situ conservation. Ex-situ conservation by means of artificial insemination (Cobey, Tarpy, & Woyke, 2013), could be a more efficient alternative to maintain the gene pool of native honey bees.
In spite of the small population size, MEL showed significantly lower FROH than CAR. This result corresponds to the current applied selection strategies of both subspecies. Within MEL, SL_CH are carefully selected to contribute to the local genetic diversity, whereas CAR can be considered as a highly specialized transboundary honey bee breed. This breed characteristic of CAR was also evident in the high-resolution population network, which successfully ascertained the numerous substructures within MEL, but failed to clearly separate CAR colonies according to their geographical origin. The ROH results, according to the observed population structure (Table 2), revealed that the population admixture of MEL was another reason for the low FROH. It also simultaneously indicated a direct relationship between admixture and ROH length in honey bees. The effect of population admixture on homozygosity patterns has already been documented in other populations, such as cattle (Purfield et al., 2012) and goats (Bertolini et al., 2018). The population admixture also had an effect on the relationship between FROH and FPED, whereas the improved value (R2 = 0.10) was significantly lower than commonly observed values in livestock, such as goats (R2 = 0.27 – 0.65) (Burren et al., 2016). Compared to other livestock populations, the paternal origin of honey bees is not known, as honey bee queens naturally mate in flight with 10 – 20 drones (polyandrous mating system) (Estoup et al., 1994; Neumann et al., 1999; Tarpy et al., 2004). Hence, the paternal origin must be estimated by restricting paternal origins to the drone-producing colonies located at the mating station. However, the proportion of foreign drones contributing to the mating remains unknown. Therefore, FPED of colonies from insufficiently isolated mating stations (with higher admixture proportions) are overestimated, while a low pedigree completeness results in lower FPEDcompared to FROH. To improve the pedigree quality of honeybees, we suggest confirming the parental origin with a marker-based parentage analysis.
Within CAR-specific homozygosity islands, we identified several genes that are directly associated with the current applied selection traits, including increased productivity, as well as reduced defensive behaviour and swarming drive (Bouga et al., 2011; Guichard et al., 2021). Based on highly selected A. m. ligustica lineages, it has already been demonstrated that RpL35 controls royal jelly production and larval growth (Ararso et al., 2018). Furthermore, the differential expression of Ndufs1 (Guo et al., 2019) and Gpdh (Seehuus, Taylor, Petersen, & Aamodt, 2013) may also increase foraging behaviour, and consequently, productivity. The reduced aggressiveness in CAR (Guichard, Neuditschko, Fried, Soland, & Dainat, 2019; Friedrich Ruttner, 1988) may be associated to the private homozygosity island containing the gene 5-ht7, coding for serotonin receptor 7, as higher serotonin levels increase the likelihood of bee stings, the ultimate defence mechanism of the colony (Nouvian et al., 2018).
The longest private homozygosity island of CAR included several cuticular protein-coding genes (CPR3 and CPR4 in particular). These proteins may be involved in the CAR specific morphotype of broader hairy stripes (Figure 1), as they effect the thickness and colouring of the exoskeleton (Costa et al., 2016; Soares et al., 2013). Another gene, Chmp1, is known to influence the veining pattern in Drosophila (Valentine, Hogan, & Collier, 2014), which might explain the morphological differences in vein patterns between the two subspecies (Bouga et al., 2011; Friedrich Ruttner, 1988).
The genes chr-BP and Uvop , embedded in private islands for CAR and MEL respectively, are both involved in the resistance to ultraviolet (UV) exposure, but reveal different adaptive mechanisms due to the ancestral geographical origin of the two subspecies. The gene crh-BP was shown to be upregulated in honey bees in response to UV exposure and heat stress (Even, Devaud, & Barron, 2012). CAR-specific homozygosity in that gene could therefore indicate local adaptation to more constant sun exposure and higher temperatures of Southern Europe. The homozygous state of the Uvop gene in MEL is associated with retinal development and the circadian rhythm (Lichtenstein, Grübel, & Spaethe, 2018), which may enable MEL to deal with seasonally more variable sun exposure. For example, within fish certain Malawi Cichlids also show differential expression of opsin genes depending on their photic environment (Parry et al., 2005). Diurnal mammal species also produce different quantities of ultraviolet-sensitive pigments depending on their ecological niche (Emerling, Huynh, Nguyen, Meredith, & Springer, 2015). Furthermore, it has recently been demonstrated that genes involved in the response to UV exposure are associated with the local adaption of horse breeds (Grilz-Seger, Neuditschko, et al., 2019).
Several characterised genes in a homozygous state for CAR and MEL shared functions associated with stress response (ATP5G2 , Nmdar1, Tpx-4, GstS1 (Alburaki, Karim, Lamour, Adamczyk, & Stewart, 2019; Watts, Williams, Nithianantharajah, & Claudianos, 2018; Yan, Jia, Gao, Guo, & Xu, 2013)), DNA integrity (Uzip , WRNexo (Ding et al., 2011; Rossi, Ghosh, & Bohr, 2010)) and immunity (CPR2, SP34, Pban, CTL5 (Badaoui et al., 2017; Lin et al., 2020; McDonnell et al., 2013; Zou, Lopez, Kanost, Evans, & Jiang, 2006)). However, at the current stage of research, it is not clear whether the homozygosity state of these genes has a positive or negative effect on the aforementioned functions. Therefore, fine-tuned gene expression studies are required to assess the selection direction within the two subspecies.
In summary, we have described a number of novel aspects to investigate the genetic diversity of honey bees that are of potential interest. Firstly, the application of colony genotypes derived from pooled honey bee workers to ascertain fine-scale population structures. Secondly, the identification of ROH segments from pool-seq data to compute genomic inbreeding of honey bee colonies. Finally, the assessment of breed-specific selection signatures by means of ROH islands. Therefore, we believe that ROH derived from pool-seq data will be of invaluable benefit to investigate complex population structures in honey bees and other insects, whereas the hard-call genotype threshold and ROH setting should be further investigated.