Attack of the dark clones the genetics of reproductive and color traits of South African honey bees (Apis mellifera spp.).

The traits of two subspecies of western honey bees, Apis mellifera scutellata and A.m. capensis, endemic to the Republic of South Africa (RSA), are of biological and commercial relevance. Nevertheless, the genetic basis of important phenotypes found in these subspecies remains poorly understood. We performed a genome wide association study on three traits of biological relevance in 234 A.m. capensis, 73 A.m. scutellata and 158 hybrid individuals. Thirteen markers were significantly associated to at least one trait (P ≤ 4.28 × 10-6): one for ovariole number, four for scutellar plate and eight for tergite color. We discovered two possible causative variants associated to the respective phenotypes: a deletion in GB46429 or Ebony (NC_007070.3:g.14101325G>del) (R69Efs*85) and a nonsense on GB54634 (NC_007076.3:g.4492792A>G;p.Tyr128*) causing a premature stop, substantially shortening the predicted protein. The mutant genotypes are significantly associated to phenotypes in A.m. capensis. Loss-of-function of Ebony can cause accumulation of circulating dopamine, and increased dopamine levels correlate to ovary development in queenless workers and pheromone production. Allelic association (P = 1.824 x 10-5) of NC_007076.3:g.4492792A>G;p.Tyr128* to ovariole number warrants further investigation into function and expression of the GB54634 gene. Our results highlight genetic components of relevant production/conservation behavioral phenotypes in honey bees.

Introduction

Modern western honey bees (Apis mellifera) show substantial genetic and phenotypic variation across their extensive geographic range [1]. They occur naturally in Europe, the Middle East, western Asia, and Africa, where the species is composed of between 25–35 subspecies [2-4]. This bee has been spread outside its native range to the Americas, Australia, New Zealand, and other locations globally, largely due its ability to produce honey and its use as the principal pollinator of a variety of agricultural crops.

Two subspecies of western honey bees, A.m. scutellata and A.m. capensis, are among those endemic to the Republic of South Africa (RSA) [5]. Apis mellifera scutellata is a light-colored phenotype and is adapted for survival in hot and arid climates in central and southern Africa [6]. It also displays behavioral traits that many beekeepers outside its native range consider undesirable. These include excessive swarming (colony-level reproduction), absconding (complete nest abandonment), usurpation (swarm takeover of another colony) and heightened defensiveness [6-8]. This honey bee subspecies was introduced into Brazil in the 1950’s in an effort to improve the Brazilian beekeeping industry [9]. It hybridized with local stocks of European-derived honey bees, becoming known as “Africanized” or “killer” bees. They are now considered invasive throughout South America, Central America and southern regions of North America [7, 10].

Apis mellifera capensis is a darker colored honey bee subspecies found in the Fynbos region of RSA, where the climate is Mediterranean with rainy winters. In contrast to A.m. scutellata, this bee can act as a social parasite, given its workers can reproduce via thelytoky [6], a type of parthenogenesis in which female offspring can result from unfertilized eggs. This trait allows some worker bees to develop into pseudoqueens with semi-developed spermathecae, that remain unused, and a larger-than-normal number of ovarioles [11-15]. These worker bees, then, can fly into neighboring hives and replace the queens contained within, becoming the reproductive in the nest [16]. Interestingly, colonies headed by A.m. capensis workers are doomed, as laying workers cannot maintain the egg output of that of a normal queen. The colonies eventually dwindle and die, resulting in the ‘capensis calamity’ that has plagued the South African beekeeping industry in the past [17].

Despite the perceived drawbacks associated with these bees outside their native range, beekeepers in RSA keep both subspecies for management purposes. Nevertheless, the potential movement of both bee subspecies beyond where they currently occur remains a concern of beekeepers and regulatory officials in many areas globally. These concerns have led to the search for better methods to identify both bee subspecies and their hybrids quickly and reliably. Recently developed techniques based on the reduction of genome complexity, such as Genotyping by Sequencing (GBS), have the potential to provide a large number of SNPs in understudied genomes, enabling genetic diagnostics for monitoring these two subspecies [18]. Despite genomic studies on various honey bee subspecies, the genetic basis of important phenotypes found in A.m. scutellata and A.m. capensis remain poorly understood, though progress has been made with the thelytoky trait [19-25]. We have the opportunity to fill this gap given recent work [26] that used traditional morphometric techniques to identify populations of both bees from samples collected in RSA.

In the present study, we performed a genome wide association study (GWAS) on three traits (number of ovarioles, tergite and scutellar plate color) measured in 464 A.m. capensis, A.m. scutellata and hybrid individuals (S1 Table) [26]. These same bees had been examined previously using GBS [18]. Apis mellifera capensis is known to be darker and have a greater number of ovarioles per ovary than does A.m. scutellata. Accordingly, the GWAS allowed us to determine what chromosomal regions are most associated with these phenotypic traits. The detected associations provide improved understanding of the genetic basis of phenotypic and behavioral differentiation between A.m. capensis and A.m. scutellata from RSA.

Results

GWAS associates traits mainly to two chromosomes

The GWAS resulted in significant associations to markers on chromosomes LG1, LG2, LG7, LG9 and LG10. Thirteen markers were significantly associated to at least one trait (P ≤ 4.28 x 10−6): one for ovariole number, four for scutellar plate and eight for tergite color. A total of 10 genes are annotated in candidate regions determined by markers within r ≥ 0.2 to the most significant marker, and adjacent genes (Fig 1 and Table 1).

Fig 1

Manhattan and QQ plots of the respective genome wide association study for a. tergite color ranked-transformed; b. scutellar plate color ranked-transformed; and c. ovariole number rank-transformed.

Respective annotated genes within the shared regions in chromosomes LG1 and LG7, as well as genes possessing non-synonymous variants (in bold), are also shown. The red line represents the Bonferroni corrected threshold value of P ≤ 4.28 x 10−6, and markers above this line are significantly correlated to the respective trait.

Table 1

Genome wide association study traits, significant markers, respective chromosome (Chr) location, number of base pairs, statistical information, and within region/nearby annotated genes.

Trait	Marker ID	Chr	Base Pair	P-value	Genes
Ovariole Number	S1_108729877*	7	4497718	1.8241 x 10–7	GB54634
					Tk
Scutellar Plate	S1_14074192**	1	14074192	1.12149 x 10–9	GB46427 GB46500 GB46429
	S1_14077754**	1	14077754	1.03746 x 10–7
	S1_108729877*	7	4497718	1.03692 x 10–7	GB54634
					Tk
	S1_108786716***	7	4554557	1.83841 x 10–7	–
Tergite Color	S1_14074192**	1	14074192	1.19236 x 10–13	GB46427 GB46500 GB46429
	S1_14077754**	1	14077754	1.79334 x 10–9
	S1_14080360**	1	14080360	5.6368 x 10–8
	S1_15042195	1	15042195	6.60343 x 10–7	GB52133
	S1_34570644	2	4677136	8.12408 x 10–7	–
	S1_108786716***	7	4554557	6.28334 x 10–8	–
	S1_132740677	9	1742429	5.37694 x 10–7	GB43750 GB43751 GB43755 Tpx-4
	S1_147958987	10	5840186	1.25779 x 10–6	–

*An asterisk (*) represents the marker sharing the same candidate region.

–A minus (–) means no annotated genes occurred in the candidate region.

A frameshift and a nonsense mutation are associated to color and ovariole number

Functional inspection of annotated genes within each candidate region indicated two genes with coding variants. The likely candidate gene for tergite and scutellar plate color is GB46429, mycosubtilin synthase subunit C, also known as Ebony, a non-ribosomal peptide synthetase, which also has sequence similarities to microbial enzymes [27]. This gene shares 46.99% (EnsemblMetazoa release 103, LOC409109) [28] of its sequence with the Drosophila melanogaster Ebony gene. A deletion identified by the GBS pipeline in GB46429 (NC_007070.3:g.14101325G>del) (R69Efs*85) leads to an early stop codon and truncates the normal amino acid sequence from the predicted 860aa to only 85 amino acids.

A single variant was found for ovariole number within the coding region of GB54634. The nonsense SNP (NC_007076.3:g.4492792A>G;p.Tyr128*) causes a premature stop, shortening the protein by two of the six predicted exons (45% of the protein sequence) (Fig 2).

Fig 2

Predicted protein structure (Phyre2) for both wild type and discovered variants of genes.

a.GB46429 (Ebony), correlated to both scutellar plate and tergite color phenotypes., and b. GB54634, correlated to ovariole number.

The distribution of causative variants demonstrates that the mutant form is significantly associated to phenotypes in A.m. capensis, while the wildtype locus is associated to A.m. scutellata phenotypes (Fig 3). No coding variants were discovered in our GBS dataset for the other annotated genes within each candidate region; yet these could hold biological effects of interest for the honey bee.

Fig 3

Allele distribution for variants discovered in GB54634 and GB46429 (Ebony), as well as respective color phenotypes.

a. Ruttner [2] ranking for tergite color, also applied to scutellar plate phenotyping [19]; b. Allelic distribution of the NC_007070.3:g.14101325G>del;p.R69Efs*85 variant for tergite color and c. scutellar plate color; d. Allelic distribution of the NC_007076.3:g.4492792A>G;p.Tyr128* variant for ovariole number; and e. Individuals from the Apis mellifera scutellata (above) and A.m. capensis (below) representing the variation in color [29].

Discussion

Color variation in honey bees may have diverse biological implications [30]. For example, Gloger’s rule states that coloration changes according to environmental effects, and species tend to be darker in hot and humid environments [31]. Yet, this rule might not apply to the present case, as the A.m. scutellata individuals were collected from, on average, warm semi-arid zones, while the A.m. capensis or hybrid samples came from cooler, Mediterranean or cool subtropical zones; yet, A.m. scutellata had significantly lighter phenotypes both in tergite and scutellar plate color [26]. Additionally, previous thelytoky genome mapping efforts pointed to a locus near GB46429 [24]. However further inspection into expression and gene function demonstrated that this gene has no apparent effect on the mode of parthenogenesis in the honey bee, but segregates according to subspecies and color [32].

In Drosophila melanogaster, the orthologous gene to GB46429 is Ebony (named after the mutant phenotype): darker Drosophila flies have lower expression, while lighter individuals have normal to high expression of Ebony [33]. Variants in Ebony also contribute to diverse phenotypic variations including behavioral, neurologic, locomotor, and visual ability [34, 35]. Some Drosophila Ebony mutants’ electroretinograms lacked the on- and off-transients of light response [36, 37]. Most importantly, Ebony participates in dopaminergic neuron function, metabolizing dopamine into N-β-alanyl dopamine (NBAD) [38]. We discovered a single nonsense variant in GB46429 (Ebony) significantly associated to color phenotypes in both honey bee subspecies and hybrids of the two. Furthermore, this variant severely impacts the predicted protein structure and may lead to loss-of-function of this protein. Consistent with our findings, loss-of-function Ebony mutants in Drosophila accumulate circulating dopamine, which is then directed to other pathways [39]. In the honey bee, increased dopamine levels correlate to ovary development in queenless workers, as the queen mandibular pheromone (QMP) regulates dopamine pathways in the worker bees [40, 41]. In A.m. capensis, pheromonal dominance allows for parasitic behavior, even in the presence of an A.m. scutellata queen [42]. We postulate that the mutant GB46429 causes a darker pigmentation phenotype and may play a role in dopaminergic pathways and parasitic behavior in A.m. capensis. This gene’s contribution to behavioral and reproductive traits in honey bees is worthy of further investigation.

Previous work evaluating quantitative trait loci (QTLs) impacting the number of ovarioles in honey bees resulted in a significant QTL on LG11 [43]. Although our GWAs did not associate any markers on LG11 to ovariole number, this difference in findings could be due to population genetic differences as the LG11 QTL resulted from Africanized Honey Bees (AHB) collected in Arizona, USA, compared to European Honey Bee samples collected from US commercial colonies.

Unfortunately, there is little information of the function and expression of the GB54634 gene in honey bees even though we found the significant (P = 1.824 x 10−5) allelic association of the NC_007076.3:g.4492792A>G;p.Tyr128* variant to ovariole number (Fig 3). The GB54634 gene was tagged by genomic sweeps associated to social parasitic behavior [44] and A.m. capensis versus A.m. scutellata differentiation [18], although candidate genes for thelytoky phenotype recently reported do not implicate GB54634 in this specific phenotype [24, 25]. Additionally, this uncharacterized protein (LOC725260 isoform X1) does not seem to differ in expression and splicing in the presence or absence of queen pheromones [45]. Yet, this expression analysis was conducted in an uncharacterized A. mellifera subspecies; thus, findings could be different for A. m. capensis and A. m. scutellata. Given the correlation reported here, further investigation into association of this variant to social parasitic traits such as the number of ovarioles, as well as possible pleiotropic effects, warrants additional exploration and biological characterization of GB54634.

Although we did not discover coding sequence variants for the other genes within candidate regions, biological functions related to A. m. capensis phenotypes may be of interest for future analysis. For instance, candidate regions for the color traits GB43750 (prefoldin subunit 5) are located within a haplotype associated to high altitude adaptation in A. m. scutellata [46].

TK (prepro-AmTRP or tachykinin) was also found in the candidate region associated to ovariole number in the GWAS. Previously implicated in female-related behavior, the expression levels of prepro-AmTRP are present only in the brain of female bees (queens and workers) and show lower expression levels according to labor division (lower in younger/nurse bees, higher in queens and foragers) [47]. The tachykinin neuropeptide also controls metabolic and desiccation responses in Drosophila [48, 49] and is related to aggression in other insects, such as the Leaf-Cutting Ant Acromyrmex echinatior [50]. Other AmTRP neuropeptides are implicated in the defensive behavior of Africanized honey bees [51].

Several genomic regions are likely involved in ovariole number, a social parasitism-related phenotype of A.m. capensis colonies [44]. The GB46427 gene (LOC409096) within the ovariole number LG1 candidate region is implicated in parasitism behavior and was deemed the thelytoky gene [24], also demonstrating Log2-fold differential expression of 3.24 to 4.68 between thelytokous A.m. capensis and arrhenotokous A.m. scutellata [24, 44]. A non-synonymous variant (p.Thr400Ile) likely responsible for this differential expression was suggested as the sole change responsible for thelytoky in worker bees [24]. Our GBS dataset did not possess any variants within the coding region of this gene; thus, we could not evaluate the phenotypic repercussions. Furthermore, GB46500 (LOC724495 or Ethr) is also in linkage with GB46427 [24]. In Drosophila, lower levels of the hormones transcribed by Ethr halt oogenesis and ovulation during nutritional or heat stress [52]. Therefore, its effects on honey bee social parasitism might be of biological relevance, though we could not find coding variants for this gene.

Conclusions

We associated genomic regions with important biological phenotypes as tergite color, scutellar plate color, and ovariole number within A.m. capensis and A.m. scutellata populations from RSA. Among the 28 candidate genes identified, Ebony, within the tergite color candidate region on chrLG1, possessed a variant predicted to alter protein structure significantly. Furthermore, non-functional variants of Ebony impacting pigmentation are well-documented in other insect species. Although the candidate variant correlated to ovariole number is in an uncharacterized gene, further investigations into its function are warranted given its biological implications. Our results help pave the way for the development of marker-assisted selection and diagnostic genetic differentiation in the honey bee and highlight potentially production/conservation relevant pleiotropic behavioral phenotypes.

Material and methods

Honey bee samples

The samples included 464 adult worker honey bees collected from RSA in 2013 and 2014. The samples were collected from managed colonies of A. mellifera with permission granted by the owner beekeepers (see Acknowledgements). Location data for the samples, including GIS coordinates, can be found in S1 File. Combined morphometrics, SNP, microsatellite and mitochondrial DNA data were used to determine that 73 bees were A.m. scutellata, 234 were A.m. capensis and 158 were hybrids of the two subspecies [18, 26, 53, 54]. Phenotyping methods as described in [19] determined morphometric phenotypes that significantly differed between the two subspecies of honey bees. We utilized the following traits in a GWA: number of ovarioles, pigmentation of abdominal tergite (A3) and pigmentation of the scutellar plate (Fig 4 and S1 File). The distribution for these quantitative phenotypes within the 464 samples was not normal; thus, we normalized the data prior to the GWAs using a Rank normalization on JMP®, Version 15 (SAS Institute Inc., Cary, NC, 1989–2019).

Fig 4

Distribution of morphometric phenotypic traits per subspecies, representing a. Ovariole Number quantile, b. Tergite Color quantile and c. Scutellar Plate color ranking.

Red represents increased number of ovarioles (a) or lighter phenotypes (b and c), while blue represents lower number of ovarioles (a) and darker phenotypes (b and c). The visual distribution seems to correlate to the subspecies or hybrid geographical distribution.

Genotyping and SNP QC

DNA extraction, library construction, sequencing and quality control criteria were conducted by The Genomic Diversity Facility at Cornell University. The GBS methods were previously described [18], and resulted in an average of 70,475 SNPs per individual sample. We filtered GBS SNPs (coded as major/minor allele) using VCFtools version 0.1.15 [55] and the following criteria: (1) no more than two alleles, neither of which was a gap allele, (2) a minor allele frequency (MAF) of at least 5%, (3) no more than 92% missing data, (4) mapped to one of the 16 assembled A. mellifera chromosomes in the Amel4.5 assembly [56] and (5) with an index of panmixia (FIT) of at least −0.2. After quality control, 20,006 SNPs were left. We then imputed missing genotypes for the 20,006 loci using Beagle 4.1 [57] with a window and overlap of 500 and 50 sites, respectively. After imputation, the SNPs were again filtered for MAF of at least 5%, resulting in a total of 11,656 SNPs retained per individual. The resulting VCF file was converted to PLINK format/binary ped format with the—recode—make-bed command in PLINK version 1.90b3.39 [58].

Genome wide association study

We performed a GWAS using a mixed linear model (MLM) analysis with the interrogated SNPs falling on the same chromosome as the given candidate SNP excluded from the genetic relationship matrix calculation (—mlma-loco) in GCTA ver. 1.25.2 [59]. A genetic relationship matrix (GRM) was included in the MLM analysis to compensate for population structure within the sample. We utilized a Bonferroni corrected threshold of P ≤ 1.429 × 10−6 as the significance cutoff based on 11,656 SNPs tested and the three traits analyzed (α = 0.05). We visualized GWA results in JMP®, Version 15 (SAS Institute Inc., Cary, NC, 1989–2019).

Identification of candidate genes and functional variants

Markers above Bonferroni correction were inspected for supporting linkage (r) in PLINK (—chr [ChromosomeNumber]—r2—ld-snp [MarkerID]—ld-window-r2 0.00—ld-window 100000). Loci with a r ≥ 0.2 to the lowest p-value SNP defined the boundaries of candidate regions considered for further analysis [60, 61]. We also evaluated genes adjacent to each candidate region, determined using the NCBI/GenBank annotation GCF_000002195.4 /GCA_000002195.1 (Amel 4.5) [56]. Gene function was also reported based on its homology to functionally characterized genes from the A. mellifera genome (Amel 4.5) using the EnsemblMetazoa database (release 103) [28] and a comprehensive scientific literature search on other Hymenoptera order members [62, 63].

Visual inspection of genomic regions for polymorphisms within coding regions was performed on the unfiltered, not imputed, GBS generated. vcf file, aligned, and uploaded to NCBI Apis mellifera 4.5 (accession number GCF_000002195.4), coded as major/minor allele. For candidate mutations, we evaluated protein impact using Phyre2, modeling both the wild type and the sequence containing mutation(s) [64]. Allelic association of causative polymorphisms to traits was performed on JMP®, Version 15 (SAS Institute Inc., Cary, NC, 1989–2019) using ANOVA, with the significance threshold set to P ≤ 0.00833 based on multiple tests per allele (0.05/6).

Phenotypic information, geographical coordinates and candidate variant genotypes for samples used in this study.

Ovary number value, scutellar plate and tergite color scores and respective rank transformations, as well as respective combined probable Apis mellifera subspecies ID and candidate variant genotypes per sample.

(XLSX)

Click here for additional data file.

Candidate variant distribution of alleles per subspecies.

(PDF)

Click here for additional data file.

Sample Apis mellifera subspecies assignment per source information.

Hybrid = A cross between A.m. scutellata and A.m. capensis. NA = bee samples from that location were not included in the respective analysis. The “Combined Probable ID” is inferred from the most common identification (ID) made across the four referenced studies and it parallels the identifications assigned using SNPs.

(PDF)

Click here for additional data file.

Transfer Alert

This paper was transferred from another journal. As a result, its full editorial history (including decision letters, peer reviews and author responses) may not be present.

8 Oct 2021

PONE-D-21-27869Attack of the Dark Clones The genetics of reproductive and color traits of South African honey bees (Apis mellifera spp.)PLOS ONE

Dear Dr. Patterson Rosa,

Thank you for submitting your manuscript to PLOS ONE. The reviewers and I found the topic and results of this study interesting but there are several concerns that need to be addressed. Most importantly, the major concerns raised by the second reviewer about samples, data availability, and the impact of imputation are certainly critical to address, in addition to the other concerns of both reviewers. Thus, while I decided to recommend "minor revisions" because the changes in the manuscript itself would not be extensive, this should certainly not be mistaken as a decision that publication can be guaranteed at this point.

Please submit your revised manuscript by Nov 22 2021 11:59PM. If you will need more time than this to complete your revisions, please reply to this message or contact the journal office at plosone@plos.org. When you're ready to submit your revision, log on to https://www.editorialmanager.com/pone/ and select the 'Submissions Needing Revision' folder to locate your manuscript file.

Please include the following items when submitting your revised manuscript:If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter. Guidelines for resubmitting your figure files are available below the reviewer comments at the end of this letter.

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). You should upload this letter as a separate file labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. You should upload this as a separate file labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. You should upload this as a separate file labeled 'Manuscript'.

If applicable, we recommend that you deposit your laboratory protocols in protocols.io to enhance the reproducibility of your results. Protocols.io assigns your protocol its own identifier (DOI) so that it can be cited independently in the future. For instructions see: https://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols. Additionally, PLOS ONE offers an option for publishing peer-reviewed Lab Protocol articles, which describe protocols hosted on protocols.io. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols.

We look forward to receiving your revised manuscript.

Kind regards,

Olav Rueppell

Academic Editor

PLOS ONE

Journal requirements:

When submitting your revision, we need you to address these additional requirements.

1. Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at

https://journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and

https://journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf

2. In your Methods section, please provide additional details regarding the source of the bees used in your study. If they was collected please include geographic coordinates of your field collection site if available and provide the permits you obtained for the work. Please ensure you have included the full name of the authority that approved the field site access and, if no permits were required, a brief statement explaining why. If this was purchased then please include the names of the purchasing sources (e.g., stores, markets, suppliers), if available, as well as any further details about the purchased items (e.g., lot number, source origin, description of appearance) to ensure reproducibility of the analyses. For more information regarding PLOS' policy on materials sharing and reporting, see https://journals.plos.org/plosone/s/materials-and-software-sharing#loc-sharing-materials.

3. Thank you for stating the following in the Acknowledgments Section of your manuscript:

“This work was supported through various cooperative agreements provided by the United States Department of Agriculture, Animal and Plant Health Inspection Service and by the Florida Department of Agriculture and Consumer Services through the guidance of the Honey Bee Technical Council. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.”

We note that you have provided additional information within the Acknowledgements Section that is not currently declared in your Funding Statement. Please note that funding information should not appear in the Acknowledgments section or other areas of your manuscript. We will only publish funding information present in the Funding Statement section of the online submission form.

Please remove any funding-related text from the manuscript and let us know how you would like to update your Funding Statement. Currently, your Funding Statement reads as follows:

“The author(s) received no specific funding for this work.”

Please include your amended statements within your cover letter; we will change the online submission form on your behalf.

4. Please include captions for your Supporting Information files at the end of your manuscript, and update any in-text citations to match accordingly. Please see our Supporting Information guidelines for more information: http://journals.plos.org/plosone/s/supporting-information.

5. Please review your reference list to ensure that it is complete and correct. If you have cited papers that have been retracted, please include the rationale for doing so in the manuscript text, or remove these references and replace them with relevant current references. Any changes to the reference list should be mentioned in the rebuttal letter that accompanies your revised manuscript. If you need to cite a retracted article, indicate the article’s retracted status in the References list and also include a citation and full reference for the retraction notice.

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Yes

Reviewer #2: Yes

**********

2. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

Reviewer #2: No

**********

3. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: No

**********

4. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

**********

5. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: This is an interesting paper on a genome wide association study in Apis mellifera capensis, Apis mellifera scutellata and hybrid individuals looking for associations to the number of ovarioles and two color-related traits: tergite and scutellar plate color. The authors discovered two possible genetic variations associated to these traits, a deletion in the ebony gene and a premature stop on gene GB54634.The experiments seem to be performed accurately and are well presented. The subject is put into the context of similar studies and is discussed appropriately. In my opinion, the manuscript would make a useful addition to the literature. Only a few minor suggestions to further improve the manuscript came into my mind.

Specific comments:

Minor issues:

line 102, “Bustamente et al. (2020)”: Enter the reference number as in the references list! (21)? See also line 249.

line 108, “Eimanifar et al. (2018a)”: Reference number! (18)?

line 137 or elsewhere: When the Ebony protein is introduced, it should be mentioned that Ebony is a non-ribosomal peptide synthetase, which also has sequence similarities to microbial enzymes (Hovemann et al., 1998).

lines 144/145, line 147, “Figure 2”: Check Figure numbers! Protein structures are given in Figure 3.

line 153, line 157, “Figure 3”: Check Figure numbers! Description matches Figure 2.

lines 158/159, “Ruttner (1988)” and “Bustamente et al., (2020)”: Enter reference numbers as in the references list! (2) and (21), respectively?

lines 177-178, Variants in Ebony also contribute to diverse phenotypic variations including … visual ability: In fact, it was shown some time ago that the electroretinogram (ERG) of certain Drosophila ebony mutants lacked the on- and off-transients of the light response (Hotta and Benzer, 1969; Heisenberg, 1971).

lines 178-180, “Ebony participates in dopaminergic neuron function, metabolizing dopamine into N-β-alanyl dopamine (NBAD)”: That is correct, but not the whole truth. Drosophila ebony mutants also have reduced histamine content in the head and are unable to convert histamine into its β-alanine conjugate, carcinine (Borycz et al., 2002). In vitro assays with heterologously expressed and purified Drosophila Ebony protein have shown later that the activation and binding of β-alanine occurs in a peptide-synthetase like manner. Furthermore, enzymatic activity has been observed not only with dopamine but also with histamine and other biogenic amines as substrates (Richardt et al., 2003).

line 184, line 227, “drosophila”: Drosophila should be capitalized and italicized.

lines 192/193, “… there is little information of the function and expression of the GB54634 gene in honey bees”: Is anything known about the function or the expression pattern of the homologous protein in the model organism Drosophila? Is there information about this on FlyBase or FlyAtlas?

line 224, “Aumer et al., 2019”: Reference number! (24)?

line 272, “Browning and Browning 2016”: Not in references list!

line 406, “HÄRtel”: Type “Härtel”.

Reviewer #2: The authors have re-analyzed one of their previous data sets to identify genes associated with body colour and ovary activation in African honey bees. I appreciate the approach and enjoyed reading the paper. I do have some concerns that need to be addressed prior to publication.

Major comments:

1) Do the authors have (or need) a permit to use the samples they have collected? I noticed the authors are in the US but the samples come from South Africa.

2) Data Accessibility. The authors need to upload their raw genomic data (fastq) to a repository. The results can't be replicated otherwise. The authors uploaded a vcf file to dryad but, as far as I can tell, the phenotypic data is not uploaded and with the vcf . This should be included in the raw data upload.

3) Imputation. I think the authors need to discuss and test the impacts of imputation on their data set. I have specific points on this below but the authors imply that up to 92% of the genotypes at a site were imputed. It would be worthwhile to convince the reader that imputation didn't impact your overall findings and/or the extent to which results relied on imputation.

Specific Comments

Line 101 - perhaps not the case, maybe reference previous studies by Oldroyd and Moritz groups?

Line 117 - previous work identifying QTLs for ovary activation found them on LG 11 (Linksvayer, Page, and others). Could you highlight why you're not finding an association here? Perhaps in your discussion? It's likely because of different populations being studied but this is still something that should be highlighted.

Line 117 - it would be unclear to the non-bee expert which subspecies is which color and/or has what ovariole count. Perhaps this can be elaborated on in a figure? Additionally, it would be worthwhile to see the distribution of associated SNPs within each subspecies (with and without imputation). Were these sites also found to have high Fst in previous, referenced, studies?

Line 117 - You don't elaborate much on the overlap among the traits. Do you have a quantitative relationship between body color index and ovary number? Does this relationship also associate with the associated GWAS sites.

Line 127 - your significance cut off is not explained in detail here nor in the methods. Is Bonferroni appropriate here and did you select the correct number of independant comparisons? You likely have fewer independent comparisons than you have used.

Line 165- I don't know if I would say that coloration is not a trait of interest in beekeepers. In the US, beekeepers have a long history of coloration preference in their bees (e.g. the 'three gold lines' that were preferred in Italian stocks, during early importation).

Line 168 - I'm not an ecologist, but is this Bergmann's rule? I thought that was for body size? It might be Gloger's Rule?

Line 174 - Is (24) the correct reference? Did the authors intended this one: doi.org/10.1093/molbev/msz100

Line 269 - It seems striking to allow a site to have 92% of the data and be imputed. I am curious to know the effect of imputation on this data set?

Line 271 - the authors don't explain FIT nor why they are trimming based on it.

Line 274 - 11,656 SNPs retained per individual seems like very low coverage for a genome with such a high recombination rate and very sparse for imputation given the high recombination rate. With so few SNPs per individual, how do the authors extract associated sites across all individuals (Table 1)?

Figure 4 has no description. Figure 2 and 3 descriptions seem to be swapped.

**********

6. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email PLOS at figures@plos.org. Please note that Supporting Information files do not need this step.

3 Nov 2021

Reviewer #1: This is an interesting paper on a genome wide association study in Apis mellifera capensis, Apis mellifera scutellata and hybrid individuals looking for associations to the number of ovarioles and two color-related traits: tergite and scutellar plate color. The authors discovered two possible genetic variations associated to these traits, a deletion in the ebony gene and a premature stop on gene GB54634.The experiments seem to be performed accurately and are well presented. The subject is put into the context of similar studies and is discussed appropriately. In my opinion, the manuscript would make a useful addition to the literature. Only a few minor suggestions to further improve the manuscript came into my mind.

R: We would like to thank reviewer #1 for the thoughtful and much appreciated comments throughout the review process.

Specific comments:

Minor issues:

line 102, “Bustamente et al. (2020)”: Enter the reference number as in the references list! (21)? See also line 249.R: We addressed this comment as suggested.

line 108, “Eimanifar et al. (2018a)”: Reference number! (18)?

R: We addressed this comment as suggested.

line 137 or elsewhere: When the Ebony protein is introduced, it should be mentioned that Ebony is a non-ribosomal peptide synthetase, which also has sequence similarities to microbial enzymes (Hovemann et al., 1998).

R: We have added the information, thank you for your suggestion. This is a very good reference.

lines 144/145, line 147, “Figure 2”: Check Figure numbers! Protein structures are given in Figure 3.

line 153, line 157, “Figure 3”: Check Figure numbers! Description matches Figure 2.

R: Thank you, we have fixed the figure labels accordingly.

lines 158/159, “Ruttner (1988)” and “Bustamente et al., (2020)”: Enter reference numbers as in the references list! (2) and (21), respectively?

R: Edited as suggested

lines 177-178, Variants in Ebony also contribute to diverse phenotypic variations including … visual ability: In fact, it was shown some time ago that the electroretinogram (ERG) of certain Drosophila ebony mutants lacked the on- and off-transients of the light response (Hotta and Benzer, 1969; Heisenberg, 1971).

R: Thank you for the suggestion, we have added it to the text – this is interesting and possibly relevant for the manuscript.

lines 178-180, “Ebony participates in dopaminergic neuron function, metabolizing dopamine into N-β-alanyl dopamine (NBAD)”: That is correct, but not the whole truth. Drosophila ebony mutants also have reduced histamine content in the head and are unable to convert histamine into its β-alanine conjugate, carcinine (Borycz et al., 2002). In vitro assays with heterologously expressed and purified Drosophila Ebony protein have shown later that the activation and binding of β-alanine occurs in a peptide-synthetase like manner. Furthermore, enzymatic activity has been observed not only with dopamine but also with histamine and other biogenic amines as substrates (Richardt et al., 2003).

R: Thank you for the information you provided! We would like to include it in the text; yet, to maintain the logical flow and make the information presented in the manuscript more direct, we decided to prioritize the influence of Ebony on dopamine, as this is a crucial part of the behavioral modifications cited by other authors.

line 184, line 227, “drosophila”: Drosophila should be capitalized and italicized.

R: We addressed this comment as suggested.

lines 192/193, “… there is little information of the function and expression of the GB54634 gene in honey bees”: Is anything known about the function or the expression pattern of the homologous protein in the model organism Drosophila? Is there information about this on FlyBase or FlyAtlas?

R: Unfortunately, it is the case for all other model organisms (and non-model Insecta) we have researched. It is an uncharacterized gene with no further information that could provide us any clues on its function https://useast.ensembl.org/Drosophila_melanogaster/Gene/Summary?db=core;g=FBgn0034808;r=2R:22993669-22998884;t=FBtr0071988

http://flybase.org/reports/FBgn0034808.html

line 224, “Aumer et al., 2019”: Reference number! (24)?

R: We addressed this comment as suggested.

line 272, “Browning and Browning 2016”: Not in references list!

R: We addressed this comment as suggested.

line 406, “HÄRtel”: Type “Härtel”.

R: We addressed this comment as suggested.

Reviewer #2: The authors have re-analyzed one of their previous data sets to identify genes associated with body colour and ovary activation in African honey bees. I appreciate the approach and enjoyed reading the paper. I do have some concerns that need to be addressed prior to publication.

R: We would first like to thank reviewer #2 for the thoughtful comments. !

Major comments:

1) Do the authors have (or need) a permit to use the samples they have collected? I noticed the authors are in the US but the samples come from South Africa.

R: A permit was not needed given we sampled managed colonies, with permission from the beekeeper owners. (See also the more detailed response to this question among the queries from Reviewer 1).

2) Data Accessibility. The authors need to upload their raw genomic data (fastq) to a repository. The results can't be replicated otherwise. The authors uploaded a vcf file to dryad but, as far as I can tell, the phenotypic data is not uploaded and with the vcf . This should be included in the raw data upload.

R: Thank you for noting this. We have attached the phenotype file with the manuscript. Unfortunately, for GBS studies, the file generated by the service provider (Cornell University Life Sciences Core Facility) at the time was a VCF, and not a series of fastqs representing each individual. This is the reason we uploaded the vcf. We sincerely hope this is sufficient, as we have also used the vcf (PLINK can convert VCFs to a .bed format for analysis) for our own analysis.

3) Imputation. I think the authors need to discuss and test the impacts of imputation on their data set. I have specific points on this below but the authors imply that up to 92% of the genotypes at a site were imputed. It would be worthwhile to convince the reader that imputation didn't impact your overall findings and/or the extent to which results relied on imputation.

R: Thank you for your concern with imputation. Yet, we believe the reviewer might have misinterpreted the materials and methods, as the value referenced of 92% is the threshold for missing data (missing calls or genotypes) for each SNP. This means each SNP that passed had at least 92% or more samples successfully genotyped, and this step is prior to imputation.We will address this further below.

Specific Comments

Line 101 - perhaps not the case, maybe reference previous studies by Oldroyd and Moritz groups?

We agree with the statement. We modified the text to read:

Despite genomic studies on various honey bee subspecies, the genetic basis of important phenotypes found in A.m. scutellata and A.m. capensis remain poorly understood, though progress has been made with the thelytoky trait [19-25].

Line 117 - previous work identifying QTLs for ovary activation found them on LG 11 (Linksvayer, Page, and others). Could you highlight why you're not finding an association here? Perhaps in your discussion? It's likely because of different populations being studied but this is still something that should be highlighted.

R: Thank you for your comment. We have added a comment on this study in the discussion.

“Previous work evaluating quantitative trait loci (QTLs) impacting the number of ovarioles in honey bees resulted in a significant QTL on LG11 [43]. Although our GWAs did not associate any markers on LG11 to ovariole number, this difference in findings could be due to population genetic differences as the LG11 QTL resulted from Africanized Honey Bees (AHB) collected in Arizona, USA, compared to European Honey Bee samples collected from US commercial colonies.”

Line 117 - it would be unclear to the non-bee expert which subspecies is which color and/or has what ovariole count. Perhaps this can be elaborated on in a figure? Additionally, it would be worthwhile to see the distribution of associated SNPs within each subspecies (with and without imputation). Were these sites also found to have high Fst in previous, referenced, studies?

R: Thank you for this comment – In the manuscript Introduction, we present information on the coloration and form of reproduction for both subspecies. Figures 3e and 4b/c demonstrate the variation in coloration, and Figure 4a, the variation in ovariole number by subspecies; and we have included the respective genotypes per individual sample in the supplemental file 1, as well as a S2 file with the allelic distribution per subspecies. As for the distribution of associated SNPs , polymorphisms within coding regions were discovered and analyzed based on data from the unfiltered (not imputed) GBS generated .vcf file (Material and Methods). . We did not conduct a population analysis in this study, as this was previously done and published in Eimanifar, A., Brooks, S., Bustamante, T. et al. Population genomics and morphometric assignment of western honey bees (Apis mellifera L.) in the Republic of South Africa. BMC Genomics 19, 615 (2018). https://doi.org/10.1186/s12864-018-4998-x

A full description of the Fst analysis can be found in in the referenced manuscript.

Line 117 - You don't elaborate much on the overlap among the traits. Do you have a quantitative relationship between body color index and ovary number? Does this relationship also associate with the associated GWAS sites.

R: We have performed a correlation analysis between the color traits and the ovary number, and both demonstrate a significant negative correlation (lighter individuals have fewer ovarioles; P<.0001); Yet as for the associated variants, the correlation was only significant between color and the NC_007070.3:g.14101325G>del;p.R69Efs*85 variant ( P<.0001); and between ovary number and the NC_007076.3:g.4492792A>G;p.Tyr128* variant (P=0.0003). We can add this analysis as a supplemental file if the reviewer deems necessary.

Line 127 - your significance cut off is not explained in detail here nor in the methods. Is Bonferroni appropriate here and did you select the correct number of independant comparisons? You likely have fewer independent comparisons than you have used.

R: Thank you pointing this out; we had a typo in the text. We had entered the log value for Bonferroni in the text. Bonferroni is a strict multiple testing cut frequently utilized in conservative studies with a large number of variables as in genome wide association analysis. We utilized the number of traits (3) and markers (11,656) to correct for a value of alpha of 0.05/34968 = 1.429 x 10-6 .

Line 165- I don't know if I would say that coloration is not a trait of interest in beekeepers. In the US, beekeepers have a long history of coloration preference in their bees (e.g. the 'three gold lines' that were preferred in Italian stocks, during early importation).

We agree. We deleted this text from the manuscript.

Line 168 - I'm not an ecologist, but is this Bergmann's rule? I thought that was for body size? It might be Gloger's Rule?

R: Thank you for pointing it out, we have corrected the text accordingly.

Line 174 - Is (24) the correct reference? Did the authors intended this one: doi.org/10.1093/molbev/msz100

R: Thank you for pointing it out, we have added this reference. The original reference also mentions the information in the paragraph, yet we appreciate the different perspective of the suggested paper.

Line 269 - It seems striking to allow a site to have 92% of the data and be imputed. I am curious to know the effect of imputation on this data set?

R: Thank you for your concern with imputation. Yet, we believe the reviewer might have misinterpreted the materials and methods, as the value referenced of 92% is the threshold for missing data (missing calls or genotypes) for each SNP. The number of imputed SNPs was 1600 SNPs; and these went through a second round of QC with a MAF threshold of 5%, leaving about 11,656 SNPs that were used for the GWAs.

Line 271 - the authors don't explain FIT nor why they are trimming based on it.

R: According to the findings of three published studies performed by Eimanifar et al. (References No. 18, 47, and 48), the null hypothesis of panmixia could not be rejected because the random distribution of genotypes was determined within multiple localities of the A.m. capensis population where no discernible geographic zone was recognized between two subspecies. Although the SNP data in reference no. 18 support the existence of a distinct population within A.m. scutellate, the magnitude of the panmictic gene pool across the entire geographical sampling area should not be underestimated. As a result, A.m. capensis had a large population size, with admixture evidence in the sampling areas. Therefore, we applied the panmixia index with the lowest effect to retain the most appropriate sequencings per individual for downstream analyses.

Line 274 - 11,656 SNPs retained per individual seems like very low coverage for a genome with such a high recombination rate and very sparse for imputation given the high recombination rate. With so few SNPs per individual, how do the authors extract associated sites across all individuals (Table 1)?

R: Indeed, the coverage is lower than ideal, but what was feasible given the GBS approach as originally budgeted for a smaller scale marker set suitable for the population work conducted in Eimanifar et al 2018. 11,656 SNPs is after a second round of QC and trimming. The GBS approach initially aimed to produce about 70k markers across the genome of these two subspecies. However, coverage using GBS is not at all uniform, and after applying stringent quality filters for the GWAs (32998 lost during filtering for MAF and 50 lost due to missingness) only 11,656 SNP remained). This coverage is about 1 SNP per 19,818bp. Genome-wide LD averaged 0.035 (r2) using this panel, and is comparable to that measured in other studies. The honey bee genome (~231 million bps) is considerable smaller than an average mammal species

Figure 4 has no description. Figure 2 and 3 descriptions seem to be swapped.

R: Thank you, we have fixed the figure labels. The description of figure 4 is located later in the Material and Methods, and the reviewer might not have seen it. It is:

Fig 4. Distribution of morphometric phenotypic traits per subspecies, representing a. Ovariole Number quantile, b. Tergite Color quantile and c. Scutellar Plate color ranking. Red represents increased number of ovarioles (a) or lighter phenotypes (b and c), while blue represents lower number of ovarioles (a) and darker phenotypes (b and c). The visual distribution seems to correlate to the subspecies or hybrid geographical distribution.

Submitted filename: Patterson Rosa et al 2021 Response to Reviewers.docx

Click here for additional data file.

18 Nov 2021

Attack of the Dark Clones: The genetics of reproductive and color traits of South African honey bees (Apis mellifera spp.)

PONE-D-21-27869R1

Dear Dr. Patterson Rosa,

We’re pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it meets all outstanding technical requirements.

Within one week, you’ll receive an e-mail detailing the required amendments. When these have been addressed, you’ll receive a formal acceptance letter and your manuscript will be scheduled for publication.

An invoice for payment will follow shortly after the formal acceptance. To ensure an efficient process, please log into Editorial Manager at http://www.editorialmanager.com/pone/, click the 'Update My Information' link at the top of the page, and double check that your user information is up-to-date. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org.

If your institution or institutions have a press office, please notify them about your upcoming paper to help maximize its impact. If they’ll be preparing press materials, please inform our press team as soon as possible -- no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

Kind regards,

Olav Rueppell

Academic Editor

PLOS ONE

Additional Editor Comments (optional):

Reviewers' comments:

1 Dec 2021

PONE-D-21-27869R1

Attack of the Dark Clones The genetics of reproductive and color traits of South African honey bees (Apis mellifera spp.)

Dear Dr. Patterson Rosa:

I'm pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department.

If your institution or institutions have a press office, please let them know about your upcoming paper now to help maximize its impact. If they'll be preparing press materials, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org.

If we can help with anything else, please email us at plosone@plos.org.

Thank you for submitting your work to PLOS ONE and supporting open access.

Kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Dr. Olav Rueppell

Academic Editor

PLOS ONE