Exploring the non-coding regions in the mtDNA of some honey bee species and subspecies.

The sequence of the DNA contains coding and non-coding regions. The role of the non-coding regions is not known and is hypothesized to maintain the structure of the DNA. This study aimed to investigate the structure of the non-coding sequences in honey bees utilizing bioinformatics. The non-coding sequences of the mtDNA of three honey bee species Apis dorosata, Apis florea, Apis cerana, and ten subspecies of Apis mellifera were investigated. Different techniques were utilized to explore the non-coding regions of these bees including sequence analysis, phylogenetic relationships, enzymatic digestion, and statistical tests. Variations in size and sequences of nucleotides were detected in the studied species and subspecies, but with the same nucleotide abundance (i.e. nucleotides A were more than T and nucleotides G were less than C). The phylogenetic tree based on the non-coding regions was partially similar to the known phylogenetic relationships between these bees. The enzymatic digestion using four restriction enzymes confirmed the results of the phylogenetic relationships. The statistical analysis based on numerical codes for nucleotides showed the absence of significant variations between the studied bees in their sequences in a similar way to results of neutrality tests. This study suggests that the non-coding regions have the same functional role in all the studied bees regardless of the number of nucleotides, and not just to maintain the structure of the DNA. This is approximately the first study to shade lights on the non-coding regions of the mtDNA of honey bees.

Introduction

Honey bees, genus Apis, contain about seven species, and the species mellifera contains many subspecies (Arias and Sheppard, 1996, Alqarni et al., 2011, Oleksa and Tofilski, 2015). The species mellifera occurs worldwide while the other species occur mainly in Asia. Some species are wild and not cavity-nesting bees include Apis dorosata and Apis florea while other species are cavity-nesting bees include Apis cerana and Apis mellifera. The cavity-nesting bees are kept by beekeepers and have special importance to agriculture (Breeze et al., 2011, Al-Ghamdi et al., 2016). In general, all species of honey bees have similar biological aspects: the presence of one queen in the colony, thousands of workers are responsible for all activities (Southwick and Heldmaier, 1987, Abou-Shaara et al., 2017), and the main role of drones is to mate with virgin queens (Cobey, 2007, El-Niweiri and Moritz, 2011, Heidinger et al., 2014). Approximately all bee subspecies have the same thermoregulation ability (Bernd, 1979, Villa et al., 1987, Heinrich and Esch, 1994, Jones and Oldroyd, 2006, Abou-Shaara et al., 2017). These similarities indicate the presence of common genetic characteristics between them (Abou-Shaara, 2019a). These characteristics are based mainly on the coding regions in the nuclear and the mitochondrial DNA (mtDNA).

The mtDNA sequences of some species and subspecies of honey bees were identified and are available on free online resources (Eimanifar et al., 2016, Hu et al., 2016, Eimanifar et al., 2017a, Eimanifar et al., 2017b). The availability of these sequences encourages the performance of many bioinformatics studies on honey bees (Smith and Hagen, 1996, Martimianakis et al., 2011, Ratiu et al., 2016, Eimanifar et al., 2018, Abou-Shaara and Bayoumi, 2019, Abou-Shaara, 2019a). The sequences contain coding and non-coding regions. The coding regions and the whole sequence have been used to identify the phylogenetic relationships between honey bee species and subspecies (Garnery et al., 1992, Sheppard et al., 1991, Sheppard et al., 1999, Meixner et al., 2000). However, there are no available studies on the non-coding regions of mtDNA of honey bees. The exact role of the non-coding regions is not known, and it is expected that these regions only maintain the structure of the DNA.

The similarities between honey bee species and subspecies in these non-coding regions have not been investigated. Also, the ability of the sequences of these specific regions to identify the phylogenetic relationships between bee species and subspecies is still unknown. Therefore, this study aimed to analyze the non-coding regions of the mtDNA of three honey bee species Apis dorosata, Apis florea, Apis cerana, and ten subspecies of Apis mellifera. Different analytical methods were used to identify variations in these non-coding regions between studied bees and to propose their importance.

Methods

Sequences of honey bees

The available mtDNA sequences for honey bee species and subspecies, genus Apis, were downloaded from the National Center for Biotechnology Information (https://www.ncbi.nlm.nih.gov/) (Table 1). The sequences of the identified genes of these bees were deleted from the sequences before performing the analysis to explore the non-coding regions only. To do this, the sequences were arranged in excel sheets (Abou-Shaara, 2020) and then the sequences of the coding regions were removed.

Table 1

The accession number and size (bp) of 13 species and subspecies of the genus Apis used in the study.

Species/subspecies	Size (bp)	Accession number
Apis cerana	15,895	GQ162109
Apis dorsata	15,892	NC_037709
Apis florea	15,993	KC170303
Apis mellifera capensis	16,470	KX870183
Apis mellifera lamarckii	16,589	KY464958
Apis mellifera intermissa	16,336	KM458618
Apis mellifera ligustica	16,343	NC_001566
Apis mellifera meda	16,248	KY464957
Apis mellifera mellifera	16,343	KY926884
Apis mellifera monticola	16,343	MF678581
Apis mellifera scutellata	16,288	KY614238
Apis mellifera sahariensis	16,569	MF351881
Apis mellifera syriaca	15,428	KP163643

Sequence analysis

The sequences were explored in details by calculating the percentages of nucleotides A, T, C and G. Also, two parameters based on these nucleotides, AT skew and GC skew (Perna and Kocher, 1995), were calculated. The similarity percentages between one reference species (Apis cerana) and the other bee species and subspecies were calculated to detect the variations.

Phylogenetic tree

The phylogenetic tree was constructed based on the non-coding sequences of the studied bees using MEGA7 (Kumar et al., 2016). The alignment of the sequences was done using ClustalW and IUB as DNA weight matrix with transition weight of 0.5 and 30% delay divergent cutoff. The statistical method was maximum likelihood method and using the Jukes-Cantor model (Jukes and Cantor, 1969). The phylogenetic tree was constructed after using the bootstrap method with 500 bootstrap replications.

Enzymatic digestion

Four restriction enzymes: EcoRI, EcoRV, HindIII, and NdeI were used to digest the sequences into fragments using Genome Compiler 2.2.88 (http://www.genomecompiler.com) according to Abou-Shaara (2019b). The fragments resulted from the digestion for all the studied bees were compared using NEB 100 bp ladder.

mtDNA divergence

The three species were considered as a group while the ten subspecies of Apis mellifera were considered as another group. Then, mtDNA divergence between the two groups based on the non-coding regions only was identified using DnaSP software v. 6 (Rozas et al., 2017). So, comparison between the two groups in number of polymorphic sites, total number of mutations, and nucleotide diversity was done.

Statistical variations

The method presented by Abou-Shaara and Bayoumi, 2019, Abou-Shaara, 2020 was used to perform the statistical analysis. Firstly, codes were given to the nucleotides: A = 1, G = 2, C = 3, and T = 4. Then, the non-parametric test of k independent samples (Kruskal-Wallis test) was used to identify the significant differences (P ≤ 0.05) between the sequences and nucleotides of the studied bees. The SPSS v.16 was used to accomplish the statistical analysis.

Neutrality tests

The non-coding regions of bee species and subspecies were subjected to neutrality tests based using DnaSP software v 6 (Rozas et al., 2017). The Tajima's (Tajima, 1989) and Fu and Li's (Fu and Li, 1993) tests were performed.

Results and discussion

The size of mtDNA without coding sequences for the studied bees ranged from 1160 to 5579 bp, and for subspecies of Apis mellifera ranged from 2921 to 5579 bp (Table 2). The lowest sizes were to A. dorsata and A. m. syriaca with 1160 and 2921 bp, respectively. It is clear the sizes without coding sequences are not identical between the studied bees. However, the percentages of nucleotide according to type showed the same trend for all the studied bees except A. m. sahariensis (Table 2). The percentages of nucleotide A were higher than T with positive AT skew while the percentages of nucleotide G were less than C with negative GC skew. The only exception was to A. m. sahariensis as AT skew was negative and GC skew was positive. In general, all the sequences were rich in nucleotide A followed by T then G and finally C apart from the variations in the size of the mtDNA without coding sequences. All the studied bees shared from 24.47 to 25.56% of the identical nucleotides with A. cerana (Table 2). This indicates that the non-coding regions in all the studied bees have similar identical sequences apart from the number of nucleotides.

Table 2

Components of sequences of the studied bee species and subspecies. The species A. cerana was taken as a reference to calculate the similarity percentages.

Species/subspecies	Size without genes (bp)	A%	T%	G%	C%	AT skew	GC skew	Similarity with A. cerana %
A. cerana	4859	44.10	42.77	4.86	8.27	0.01	−0.26	–
A. dorsata	1160	49.05	46.21	0.86	3.88	0.02	−0.63	25.52
A. florea	4976	45.60	42.91	4.12	7.37	0.03	−0.28	25.54
A. m. capensis	5462	45.79	42.4	4.01	7.80	0.03	−0.32	25.56
A. m. lamarckii	5579	45.81	42.48	3.94	7.77	0.03	−0.32	25.56
A. m. intermissa	5299	45.52	42.37	4.11	8.00	0.03	−0.32	25.56
A. m. ligustica	5295	45.48	42.63	4.02	7.87	0.03	−0.32	25.56
A. m. meda	5208	45.28	42.64	4.13	7.95	0.02	−0.31	25.56
A. m. mellifera	5300	45.3	42.85	4.04	7.81	0.02	−0.31	25.56
A. m. monticola	5300	45.34	42.65	4.08	7.93	0.03	−0.32	25.56
A. m. scutellata	5278	45.62	42.16	4.17	8.05	0.03	−0.31	25.56
A. m. sahariensis	4054	42.79	45.95	7.96	3.3	−0.03	0.41	25.87
A. m. syriaca	2921	44.76	41.71	4.04	9.49	0.03	−0.40	24.47

Phylogenetic relationships were detected between the studied bees based on the non-coding regions of the mtDNA (Fig. 1). The subspecies of A. mellifera were grouped together except A. m. sahariensis. Some close relationships detected in this tree are similar to previous studies using the full mtDNA sequences including the relationship between A. m. lamarckii and A. m. syriaca. The absence of close relationships between the African bee subspecies and A. m. lamarckii is supported in this study and is in line with previous studies (Franck et al., 2001, Eimanifar et al., 2017a, Abou-Shaara, 2019a) while the relationships between A. m. intermissa, A. m. monticola, A. m. scutellata, and A. m. capensis were not identical to previous studies (Eimanifar et al., 2016, Eimanifar et al., 2017a). The species A. cerana, A. florea, and A. dorosata were correctly placed away from the subspecies of A. mellifera, suggesting the possibility of using the non-coding regions to identify the phylogenetic relationships between bees. The similarities between phylogenetic relationships based on the protein coding genes from previous studies and the non-coding regions in the present study support the idea that these non-coding regions are not random or only to conserve the DNA structure, but it has a specific role. This role may be not less in its importance than the coding regions.

Fig. 1

The phylogenetic tree based on the non-coding sequences of the studied bee species and subspecies.

The enzymatic digestion of the non-coding regions yielded fragments (Fig. 2). These fragments showed clear variations between bee species A. dorosata, A. florea, and A. cerana in side and the studied subspecies of A. mellifera in the other side. Similar fragments were found between A. m. ligustica and A. m. meda. The other subspecies had 3 fragments with close sizes except A. m. sahariensis and A. m. syriaca. The variations in the sizes of the fragments can be explained by the variations in number of sequences of the non-coding regions. The similarities between closely related bees based on fragments support the genetic relationships based on the phylogenetic tree. This indicates that these non-coding regions are not random and have a specific functional role. The same trend of relationships between results of enzymatic digestion analysis and the phylogenetic relationships was found in previous studies.

Fig. 2

The enzymatic digestion of the non-coding sequences of the studied bees using Genome Compiler 2.2.88. A: A. cerana, B: A. dorsata, C: A. florea, D: A. m. capensis, E: A. m. lamarckii, F: A. m. intermissa, G: A. m. ligustica, H: A. m. meda, I: A. m. mellifera, J: A. m. monticola, K: A. m. scutellata, L: A. m. sahariensis, and M: A. m. syriaca.

The number of polymorphic sites, total number of mutations, and nucleotide diversity of three honey bee species as a group and the ten subspecies of A. mellifera as another group are shown in (Table 3). It is clear that the non-coding regions of the mtDNA of honey bee subspecies had more polymorphic sites and mutations than the bee species. However, the nucleotide diversity was higher in bee species than subspecies only by 0.08. The two groups shared 829 mutations. This analysis showed the presence of high similarities between the non-coding regions of bee species and subspecies regardless of the high number of sequences (10 sequences) in bee subspecies group than 3 sequences in bee species group.

Table 3

mtDNA divergence between three honey bee species (Apis cerana, Apis dorsata, and Apis florea) and the ten subspecies of Apis mellifera.

Parameters	Honey bee species	Honey bee subspecies
Number of sequences	3	10
Number of polymorphic sites	946	1121
Total number of mutations	1084	1800
Average number of nucleotide differences	k: 676.667	k: 588.378
Nucleotide diversity	Pi(1): 0.58333	Pi(1): 0.50722
Shared mutations	829

Statistical variations

No significant difference was found between the studied bees in the whole sequences without the coding regions regardless of variations in number and type of nucleotides (Table 4). This can be explained by the significant variations in the component of the sequence. For example, A. dorsata had the highest value in nucleotides T while A. m. syriaca had the highest value in nucleotides C; therefore there was a balance between them in the overall mean based on the whole sequence. The same situation was found in the other studied bees. Significant differences between bee species/subspecies in values of nucleotide A, T, G, and C were detected (Table 4). The mean values of nucleotides A were higher than nucleotide T while nucleotides C were higher than nucleotides G for all the studied bees except A. m. sahariensis. The statistical analysis used in this study can greatly assist in comparing sequences as shown from previous study (Abou-Shaara, 2020). This analysis suggests the high similarities between the studied bees because no significant differences were detected in their non-coding sequences as well as sharing the same trend of mean values of nucleotides (A was higher than T while C was higher than G).

Table 4

Statistical analysis of the non-coding regions of the mtDNA of the studied bee species and subspecies. The Chi-Square and p value respectively from the Kruskal-Wallis test are shown between brackets.

Species/subspecies	Sequences	Nucleotides
		A	T	G	C
A. cerana	2.50 ± 0.02	0.44 ± 0.01	0.43 ± 0.01	0.05 ± 0.003	0.08 ± 0.004
A. dorsata	2.47 ± 0.04	0.49 ± 0.01	0.46 ± 0.01	0.01 ± 0.003	0.04 ± 0.006
A. florea	2.48 ± 0.02	0.46 ± 0.01	0.43 ± 0.01	0.04 ± 0.003	0.07 ± 0.002
A. m. capensis	2.47 ± 0.02	0.46 ± 0.01	0.42 ± 0.01	0.04 ± 0.003	0.08 ± 0.004
A. m. lamarckii	2.47 ± 0.02	0.46 ± 0.01	0.42 ± 0.01	0.04 ± 0.003	0.08 ± 0.004
A. m. intermissa	2.47 ± 0.02	0.46 ± 0.01	0.42 ± 0.01	0.04 ± 0.003	0.08 ± 0.004
A. m. ligustica	2.48 ± 0.02	0.45 ± 0.01	0.43 ± 0.01	0.04 ± 0.003	0.08 ± 0.004
A. m. meda	2.48 ± 0.02	0.45 ± 0.01	0.43 ± 0.01	0.04 ± 0.003	0.08 ± 0.004
A. m. mellifera	2.48 ± 0.02	0.45 ± 0.01	0.43 ± 0.01	0.04 ± 0.003	0.08 ± 0.004
A. m. monticola	2.48 ± 0.02	0.45 ± 0.01	0.43 ± 0.01	0.04 ± 0.003	0.08 ± 0.004
A. m. scutellata	2.47 ± 0.02	0.46 ± 0.01	0.42 ± 0.01	0.04 ± 0.003	0.08 ± 0.004
A. m. sahariensis	2.52 ± 0.02	0.43 ± 0.01	0.46 ± 0.01	0.08 ± 0.004	0.03 ± 0.003
A. m. syriaca	2.48 ± 0.02	0.45 ± 0.01	0.42 ± 0.01	0.04 ± 0.004	0.09 ± 0.005
Kruskal-Wallis test	(13.11, 0.36)	(21.25, 0.04)	(25.95, 0.01)	(167.18, 0.00)	(153.90, 0.00)

Neutrality tests

No significant differences were detected in the genetic variations within the non-coding mtDNA sequences of bee species/subspecies. The Tajima's test: for synonymous sites (Tajima's D(Syn): −0.46216, not significant, P > 0.10) and for nonSynonymous sites (Tajima's D(NonSyn): 0.58981, not significant, P > 0.10). Fu and Li's tests: (Fu and Li's D* test statistic: 0.19200, not significant, P > 0.10) and (Fu and Li's F* test statistic: 0.07851, not significant, P > 0.10). These findings support the absence of significant differences between the non-coding regions of the studied bees as shown from the statistical variations based on Kruskal-Wallis test.

Conclusion:

This study suggests that the non-coding sequences have special and similar role in all the studied bees regardless of number and type of nucleotides, and not just to maintain the DNA structure. This suggestion was inferred from the nucleotide percentages, phylogenetic tree, enzymatic digestion, and the statistical tests. Nucleotides A were higher than nucleotides T while nucleotides C were higher than nucleotides G approximately in all the studied bee species/subspecies. The phylogenetic relationships based on the non-coding sequences were partially in line with the previous studies using the coding sequences or the whole sequences. Also, fragments resulted from enzymatic digestion supported the results of the phylogenetic relationships. Moreover, the statistical tests showed the absence of significant variations between species/subspecies in their non-coding sequences. The suggestion that these non-coding regions are only to maintain the DNA structure was not supported by this study because the phylogenetic analysis and the enzymatic digestion yielded informative results. Also, nucleotide percentages and the statistical analysis showed that the sequences of the none-coding regions were not random but were similar in the studied bees although the variations in the number and type of nucleotides. This study encourages researchers to conduct additional studies to explore the specific roles of the non-coding sequences in bees.