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.