Genome sequence variation in two subspecies of western honeybee, A. m. carnica and A. m. ligustica

Populations of western honeybee (Apis mellifera) show differences in morphology, physiology and behaviour as a result of adaptation to various habitats. A. m. carnica, which inhabits the South-East and Central Europe, and A. m. ligustica, which is endemic on Apennine peninsula, represent 2 closely related honeybee subspecies living in the neighbouring climatic regions. In the current study, 3,655,618 polymorphisms were identified from the whole genome sequences of 37 individual drone genomes, from A. m. carnica (n=27) and A. m. ligustica (n=10). The analysis revealed variation in genes involved in biological pathways associated with energy production and conversion, cell cycle and cytokinesis.


Introduction
The genus Apis consists of at least 9 honeybee species which are all except western honeybee (Apis mellifera) native to southern and eastern Asia (Engel, 1999). As the Asian species of honeybee live in mild tropical climate, they have no need to gather large food stores because they can meet their needs throughout the year (Ruttner, 1988). The western honeybee, which is kept by most beekeepers, is native to Europe, Africa and part of Asia, west of the Caspian sea (Webster, 2008). However, due to human activities, it has become naturalized through much of the rest of the world. Populations of honeybee show differences in morphology, physiology and behaviour as a result of adaptation to various habitats (Le Conte et al., 2008;De la Rúa et al., 2009). There are more than 26 distinct Western honeybee subspecies which can be assigned to 5 evolutionary branches: A (South and Central Africa), C (Central and Southeast Europe), M (Western and North Europe), O (Near and Middle East) and Y (Ethiopia) (Cridland et al., 2017). The honeybees of group C, A. m. carnica, A. m. ligustica, A. m. caucasica and A. m. macedonica, are variable in behaviour and colour and are adapted to various climates (Ruttner, 1988).
A. m. carnica, commonly known as the Carniolan bee, which inhabits the South-East and Central Europe, and A. m. ligustica, commonly known as the Italian bee, inhabiting Apennine peninsula, represent 2 closely related honeybee subspecies living in neighbouring climatic regions (Ruttner, 1988) with evidence of admixture in the natural hybridization zone (Dall'Olio et al., 2007). A. m. ligustica bees are yellowish-brown and exhibit low swarming tendencies, while A. m. carnica are darker in colour, less aggressive, resistant to disease and good house cleaners (Caron et al., 2013). Studies on metabolic rate and thermal limits of both subspecies suggest that A. m. ligustica populations can tolerate a high temperature exposure better than A. m. carnica populations (Kovac et al., 2014). There are also variations in the waggle dance between the 2 subspecies (Boch, 1957;Bozic et al., 2016). Several experiments showed no differences regarding learning and memory retention in both subspecies (Hoefer et al., 1975;Iqbal et al., 2019). Also, there were no differences between queens of both subspecies morphometrically or in spermathecal expression of antioxidantencoding genes (Catalase, GLOD4L, SOD1, TXN2, TXNRD1 and Vitellogenin). However, there was a trend for a higher level of GSTD1 gene expression in the spermathecae of virgin Italian queens compared to virgin Carniolan queens (Gonzalez et al., 2018). Furthermore, there were no statistical differences in spermathecal sperm counts or sperm viability between Carniolan and Italian mated queens (Gonzalez et al., 2018). Proteomic analysis of royal jelly from both subspecies revealed differences in the heterogeneity of the major royal jelly proteins, in particular MRJP3 (Li et al., 2007). The comparison of enzyme activities of Carniolan and Italian honey bees, both, uninhibited and after incubation with synergists, gave no significant differences for studied enzymes, which leads to the conclusion that the 2 populations are biochemically similar and comparable (Todeschini et al., 2017). However, in an earlier study, Rinkevich et al. (2015) tested the variation in insecticide susceptibility among honeybee subspecies and there was a large variation in sensitivity to several insecticides between A. m. ligustica and A. m. carnica. The largest variation in sensitivity was observed in neonicotinoid bioassays, where mutations or post-transcriptional modifications of nicotinic acetylcholine receptors can affect neonicotinoid sensitivity. The significant variance was noticed also in pyrethroid toxicity, which indicates variation in P450 metabolism or voltage-gated sodium channel associated with pyrethroid resistance (Rinkevich et al., 2015).
Whole genome sequences of 10 individual genomes of A. m. ligustica drones and 27 individual genomes of A. m. carnica drones were applied in the current study to identify sequence variation in subspecies. The identified polymorphisms at the whole genome level among both subspecies could be the first step for their differentiation at the molecular level.

DNA extraction and sequencing
Eight drones of A. m. carnica were collected in Slovenia in the pupae/nymph stage. The genomic DNA was extracted from the head and thorax of the frozen larvae cut into small pieces with a scalpel. The DNA extraction was performed using E.Z.N.A. Tissue DNA Kit (Omega Bio-Tek, U.S.A.) according to manufacturer's recommendations. The concentration of genomic DNA for DNA sequencing was adjusted to 50 ng/µl. Pair-end sequencing was performed on Illumina NovaSeq 6000 platform (Genewizz, Leipzig, Germany).

Gene ontology and biological pathways analysis
Bioinformatics tools g:Profiler (Raudvere et al., 2019) and EnrichmentMap (Isserlin et al., 2014) plugin for Cytoscape (Shannon et al., 2003) were used for functional annotation and enrichment analysis of a filtered set of genes. Functional annotation and enrichment analysis were based on the latest honeybee reference genome assembly (HAv3.1). All known honeybee genes were used as a background gene set in the gene enrichment analysis.

Results and discussion
The search for polymorphisms in 37 whole genome sequences of drones belonging to 2 Wallberg et al. (2014). This could be the evidence for natural selection diversifying both subspecies. In 47 out of 64 genes with missense SNPs there was one missense SNP per gene, and in 13 genes there were 2 missense SNPs. However, in 4 genes we identified 3 or more missense SNPs (Table 1). The gene with 7 missense SNPs, LOC100576610, is mitochondrial ribonuclease P catalytic subunit (mtRNase P; Figure 1). The evolutionary rate for this gene is high (2.08) in the superfamily Apoidea and lower (1.69) at the level of the class Insecta according to OrthoDB v10.1 (Kriventseva et al., 2019). The mtRNase P has been identified as an important factor in RNA processing in a wide number of eukaryotic taxa including yeasts, plants and animals (Gissi et al., 2008). Mitochondrial DNA of the majority of metazoan species encodes only a tiny proportion of proteins (<5%) required for a normal mitochondrial function. Proteins encoded by mitochondrial DNA are components of the enzyme complexes which are responsible for oxidative phosphorylation. Mitochondrial encoded proteins are synthesized by mitochondrial translation machinery, where tRNA plays a crucial role. For a normal tRNA function, cleavage of polycistronic transcripts and specific posttranscriptional modifications of the primary transcripts at both ends are required (Jackman and Alfonzo, 2013). Recently, a specific RNase P (PRORP) has been identified in human mitochondria, being responsible for the 5'-end cleavage of pre-tRNA (Holzmann et al., 2008). In addition to PRORP, the mtRNase P complex contains two additional riboprotein complexes, RNase P Protein 1 and RNase Protein2, both required for selective methylation of tRNA (Vilardo et al., 2012). The PRORP is characterized by the short signal peptide, required for the import of the protein to the mitochondria, repetitive region responsible for RNA binding and the active domain, which is metallonuclease and is responsible for RNA cleaving. The analogue to the human PRORP protein is Mulder protein in Drosophila melanogaster (Sen et al., 2016). We compared the coding region for RNase P of A. m. carnica and A. m. ligustica and identified 7 missense mutations in the coding region suggesting that RNase P is a good candidate locus for molecular discrimination of both subspecies.
Functional enrichment analysis of 64 genes with at least one missense variant revealed enrichment of GO terms ( Figure 2) and KEGG pathways associated with energy production and conversion. Enriched KEGG pathways are pyruvate metabolism, fatty acid metabolism, fatty acid elongation, fatty acid degradation, valine, leucine and isoleucine degradation, citrate cycle and glycolysis/gluconeogenesis.

Figure 2
Network representing enriched pathways based on GO terms. Edges represent genes in common between 2 pathways (nodes). The node colour corresponds to the significance of the geneset. Sources GO (BP, MF, CC) Genes in intersections of enriched biological pathways are ral GTPase-activating protein subunit beta (LOC412130), rab GTPase-binding effector protein 1 (LOC726432), regulator complex protein LAMTOR1 (LOC725789), pyruvate dehydrogenase E1 component subunit alpha (LOC551103) and trifunctional enzyme subunit beta, mitochondrial (LOC551775).
In our material we found 2 missense mutations in the RalGAPbeta subunit between A. m. carnica and A. m. ligustica, which due to involvement of ral GTPase in different cellular functions may contribute to physiological differences between the 2 subspecies. Members of the RAS-GTPase superfamily were found to have decreased expression in the solid honeybee larval tissues during the development (Chan et al., 2008). The Ral-GTPase is regulated by guanine nucleotide exchange factors and GTPase-activating proteins (GAP) which represent a heterodimeric complex composed of alpha and beta subunit (Shirakawa et al., 2009). RalGAP beta is recruited to the mitotic spindle and intercellular bridges in dividing cells and contributes to the spindle assembly checkpoint. Cells with RalGAPbeta deficiency show altered meta-to anaphase transition in humans (Personnic et al., 2014).
The pyruvate dehydrogenase (PDH) complex is a combination of 3 catalytic components, E1 and E2, which generate acetyl-coenzyme A and FAD/NAD+ dependent E3, which performs redox recycling (Perham, 1991). The PDH catalyses the oxidative decarboxylation of pyruvate and produces acetly CoA, CO2 and NADH(H + ). The PDH complex occupies a key position in the oxidation of glucose by linking gycolysis to the cycle of the tricarboxylic acid and maintains the glucose homeostasis (Patel et al., 2014). The PDH catalytic component E1 catalyses 2 consecutive steps, the decarboxylation of pyruvate and the reductive acetylation of the lipoyl groups. The gene encoding PDH belongs to the 3 duplicated glycolysis/gluconeogenesis genes in the honeybee genome (Kunieda et al., 2006).
In our material we found 4 point mutations in the catalytic component E1 of the honeybee PDH. The missense mutations could have an effect on the efficiency of the oxidative decarboxylation of pyruvate also in honeybee and so affect the efficiency of the energy metabolism. Since the activity of the E1 component affects the activity of whole PDH complex, the different genetic variants of the E1 component may play a role in adaptive traits related to energy metabolism. In some cases the deficiency of E1 subunit can cause the demyelination, resulting in neurological symptoms (Singhi et al., 2013).

Conclusions
Analysis of genomic sequences of A. m. carnica and A. m. ligustica revealed differences in the genetic background between the 2 subspecies which may be associated with phenotypic differences in metabolic rate and significant differences in insecticide sensitivity. Analysis of genetic variations in the protein coding part of the genome revealed polymorphisms fixed for different alleles in sampled populations. Among those loci, there was high proportion of missense SNPs, which might be an evidence for natural selection. Functional analysis of genes containing missense SNPs revealed pathways associated with energy metabolism, cell cycle regulation and cytokinesis.