Eye can see clearly now: Identifying the locus associated with a white‐eye mutation in honey bees ( Apis mellifera )

TL;DR

This study identifies a genetic region in honey bees linked to white-eye mutations using community science and genome analysis.

Contribution

The study discovers a specific eight-base pair region in a non-coding RNA associated with white-eye pigmentation in honey bees.

Findings

A genome-wide association scan identified a 52 kb region on chromosome 11 linked to white-eye pigmentation.

An eight-base pair region with SNPs and a deletion in a non-coding RNA is strongly associated with the white-eye phenotype.

The mutation likely affects an uncharacterized non-coding RNA (LOC100578475) involved in eye pigmentation.

Abstract

In insects, eye pigmentation is vital for various adaptive functions, including foraging, mating and predator avoidance. Due to its visible and often heritable variation, altered eye pigmentation in mutants provides an excellent model for studying biosynthetic pathways and identifying new genes involved in pigmentation. Eye mutants can also be valuable in science communication as they offer simplified examples to help the public understand complex genetic concepts. In this study, we used a community science‐based approach to identify the inheritance pattern and mutation(s) responsible for white‐eye pigmentation in honey bees. With the help of the beekeeping community, we identified a honey bee queen that produced a fraction of haploid sons (called drones) with white‐eyes. As the queen was wild‐type, we predicted that the mutation is most likely recessive to wild‐type. Using genome‐wide association and differentiation scans of wild‐type and white‐eyed drone brothers, we identified a single elevated region (52 kb) of chromosome 11. This region contains four non‐coding RNA (ncRNA) and one protein‐coding gene. We identified an eight‐base pair region with two SNPs and a four‐nucleotide deletion that are likely responsible for the phenotype. The mutation likely affects the expression and/or downstream effects of an uncharacterized ncRNA (LOC100578475). Our efforts highlight the value of community‐based science in novel gene discovery. We hope this serves, not only as a new example for the burgeoning field of honey bee functional genomics, but also as a teaching tool for both researchers and educators.

We performed a community science investigation with the beekeepers nationwide to find phenotypic mutants.Alongside a beekeeper in Kentucky, we performed a GWAS to identify a region associated with a white‐eye phenotype.We identified an eight bp region in a non‐coding RNA strongly associated with white eyes.

We performed a community science investigation with the beekeepers nationwide to find phenotypic mutants.

Alongside a beekeeper in Kentucky, we performed a GWAS to identify a region associated with a white‐eye phenotype.

We identified an eight bp region in a non‐coding RNA strongly associated with white eyes.

Figures

(a) Lateral and anterior views of an exemplar white‐eye mutant. (b) Wild‐type drone collected from the white‐eye mutant colony (Image credit Izaak Gilchrist).

Genetic mapping of white eye. (a) Manhattan plot of 100 individual brothers. Orthologs of Drosophila eye pigmentation genes are listed across the chromosomes. A peak at chromosome 11 is the only region of significant association with the phenotype. (b) A differentiation scan (F ST) across the same region of genetic differentiation between pools of white‐eyed versus wild‐type brothers (white‐eye F STmax = 0.83). (c) The fourth exon of the candidate locus (LOC100578475) contains two SNPs (bold/red) and a four‐nucleotide deletion in white‐eyed drones. This eight‐base‐pair region of variation occurs in 43 reads of white‐eye and does not occur in wild‐type reads. Lowercase letters of the transcript signify nucleotides adjacent to the transcript, and uppercase letters signify nucleotides within the transcript.

Full text

INTRODUCTION

Insects display an impressive array of colours and patterns, largely due to pigments and genetic mechanisms that differ from those in vertebrates (Futahashi & Osanai‐Futahashi, 2021; Stavenga, 2002). Insect eyes are typically dark, but may exhibit a variety of colours and patterns, with the colour determined by structural colouration (e.g., Bernard & Miller, 1968) and screening pigments in the pigment cells. Colour variation within species can result from mutations that inhibit pigment production or deposition (Stavenga, 2002). Such mutations are rare because they can reduce an insect's ability to correctly filter light (Stavenga, 2002), they can impact pigment cell activation (Ferreiro et al., 2018), and can increase the chance of damage to retinal cells (Insausti et al., 2013).

The rarity of eye mutations has not deterred investigation. Eye mutations have been important subjects in genetics for over 100 years (e.g., Morgan, 1911). Eye mutants often serve in genetic studies as markers used to track specific genetic regions in model species (Hales et al., 2015), and as markers in early genetic maps (e.g., Hiroyoshi, 1961; Rothenbuhler & Kerr, 1968). The first linkage map of Drosophila used the frequency of offspring with combinations of eye, body and wing mutations to determine distances between genes on chromosomes and identify genetic linkages (Sturtevant, 1913). This technique has been replicated outside of model systems. Eye mutations found in honey bees were used to construct the first (incomplete) linkage map for the species in 1968 (Rothenbuhler & Kerr, 1968). Beyond their use as markers, eye colour mutations can be valuable scientifically. They can provide insights into how biological pigments are synthesized (Dustmann, 1968; Kim et al., 2013; Summers et al., 1982; Xu et al., 2020). For example, in honey bees and Drosophila, the chemical pathway from tryptophan to ommochrome has been mapped with eye colour mutants (Dustmann, 1968). Recent work on the well‐characterized white locus (CG2759) of Drosophila has further explored how mutations in white can impact sexual reproduction (Xiao et al., 2017), serve as a model system for human retina neuropathies (Ferreiro et al., 2018), and provide the basis for ‘controls’ for functional genomic screens even outside of Drosophila (e.g., Heu et al., 2022). Ultimately, eye mutants and the genes underpinning them have historically informed and generated novel biological insights from the basic to the applied.

Most work to date on insect eye colour variation has come from the genetic model system, Drosophila. While this body of work has provided significant insights, it inevitably misses the unique biology of other species, particularly those outside the Drosophila genus. In the case of the honey bees (Apis mellifera), for example, eye colour mutations have been noted in the English scientific literature since at least 1931 (Beljavsky, 1931). It is often presumed that the same genes underpin eye colour variation in honey bees as in Drosophila (e.g., Dustmann, 1968). However, little progress has been made in characterizing the identified phenotypes since the pioneering work of Dustmann (1968). To better understand the underlying genetics of eye colour variation and eye pigmentation biology, we need to explore systems outside of Drosophila. The advantages of working with eye mutations in honey bees over some other systems include the large community of interested beekeepers globally and haplodiploidy. On the former, there are at least 250,000 beekeepers in the US keeping over 2.5 million colonies (NASS, 2024). This is incredible person‐power to draw from for natural observation or directed experimentation (Cooper, 2016). On the latter, haploid honey bee drones (males) are a direct meiotic product of the queen and, therefore, will not gain any paternally inherited alleles, while a given worker (diploid and fertilized female offspring of a queen) will inherit 50% of the queen's alleles and 100% of the father's alleles. Having access to a haploid state is highly advantageous for detecting recessive alleles: they will all be expressed in drones that inherit them. Access to haploids for genome sequencing also allows for higher genotype and phasing accuracy with fewer reads (Wragg et al., 2016). The honey bee is additionally useful for genetic discovery because of its small and contiguous genome build (Wallberg et al., 2019).

Here, we used a community science approach with the beekeepers of the US to locate eye mutants and characterize their genetics. We were able to identify one mutant‐producing colony whose wild‐type queen created white‐eyed and wild‐type drones. Working with the owner of this colony, we assessed patterns of inheritance and undertook a large‐scale genomic scan to identify the mutation(s) responsible. The results of this work provide insight into the function of a previously uncharacterized non‐coding RNA.

RESULTS AND DISCUSSION

Emergence ratios and inheritance of eye colour variants

Throughout 2021 and 2022, we engaged with the beekeeping community throughout the United States to identify colonies containing abherrant (mutant) drones. A beekeeper in Kentucky (PP) responded to our callout with a colony producing drones with white eyes (Figure 1a). The colony's queen and her diploid offspring had wild‐type eye colour (Figure 1b), while a fraction of her haploid drones were white‐eyed. Given the rarity of eye colour mutations and the queen's wild‐type phenotype, we hypothesized that this was the result of a single recessive mutation carried by the queen. To test this hypothesis, we induced the queen to produce a cohort of drones (see Experimental Procedures) and recorded their phenotypes at emergence. We observed 825 white‐eyed and 850 wild‐type drones. This was not significantly different from 1:1 (Fisher's exact test p > 0.05), suggesting this phenotype is underpinned by a single recessive allele.

(a) Lateral and anterior views of an exemplar white‐eye mutant. (b) Wild‐type drone collected from the white‐eye mutant colony (Image credit Izaak Gilchrist).

White eyes are likely caused by a single mutation on chromosome 11

We used two independent genomic mapping methods to identify the mutation associated with the white‐eyed phenotype. First, we performed a Genome Wide Association Study (GWAS). We performed skim‐seq (Li et al., 2021) following methods developed for honey bees (Yang et al., 2011). After testing for association with a mixed‐linear model across 30,125 unlinked loci, we identified one significant peak (−log_10_(P) > 8) associated with the white‐eye phenotype spanning approximately 52 kb on chromosome 11 (13,342,830 to 13,395,032; Figure 2a, Table S1). We next sought to identify sites fixed in frequency between the white‐eye and wild‐type drones using a pool‐seq approach (N = 60 each) (see Experimental Procedures). We identified a single highly differentiated peak between the two phenotypes occurring in the same region of chromosome 11 identified above (Figure S1; Figure 2b).

Genetic mapping of white eye. (a) Manhattan plot of 100 individual brothers. Orthologs of Drosophila eye pigmentation genes are listed across the chromosomes. A peak at chromosome 11 is the only region of significant association with the phenotype. (b) A differentiation scan (F ST) across the same region of genetic differentiation between pools of white‐eyed versus wild‐type brothers (white‐eye F STmax = 0.83). (c) The fourth exon of the candidate locus (LOC100578475) contains two SNPs (bold/red) and a four‐nucleotide deletion in white‐eyed drones. This eight‐base‐pair region of variation occurs in 43 reads of white‐eye and does not occur in wild‐type reads. Lowercase letters of the transcript signify nucleotides adjacent to the transcript, and uppercase letters signify nucleotides within the transcript.

There are no known orthologous genes from the ommochrome pathway of insects in the candidate region (Figure 2a). It contains only five annotated genes (Table S2): one protein‐coding gene (LOC410351‐ Delta) and four uncharacterized ncRNAs (Figure 2c). The authors are unaware of pigmentation defects caused by mutations in Delta, and we found no fixed sites between white‐eye and wild‐type drones in Delta. Most documented mutants in Delta are lethal, likely because of the critical role Delta plays in the Notch pathway, which regulates cell fate decisions (e.g., Vässin & Campos‐Ortega, 1987). In the retina, altering Delta expression can impact the overall morphology of a fly's eye, inclusive of pigment cells, but this is apparent at the gross anatomical level, and there has been no documented pigment deficit (Parks et al., 1995).

We identified eight fully fixed sites between phenotypes in the candidate region (Table S3). Five are present in samples of honey bees deposited on NCBI and in wild‐type samples previously sequenced by the senior author (Harpur et al., 2014), suggesting to us that these mutations are not associated with the novel phenotype. Three sites were found to be fixed in white‐eye samples and not found to be in wild‐type drones nor in any honey bee sample sequenced and publically available to date (Table S3; Figure 2c). These three sites (T > C, T > A, and a four bp deletion) are within eight bp of each other and all are found n the coding region of LOC100578475 (an uncharacterized ncRNA; Figure 2c).

Non‐coding RNA can make up the bulk of a cell's transcriptome (Uchida & Adams, 2019) but their functional roles are largely yet to be uncovered (Chen & Kim, 2024). There are examples of ncRNA playing critical roles in pigmentation pathways in insects (Livraghi et al., 2024). Previous work in Apis cerana has shown some tentative evidence of ncRNA expression being involved in body pigmentation (Abdelmawla et al., 2023). Additionally, honey bee drones express a retinal‐specific ncRNA called Nb‐1 which may have a role in visual development (Tadano et al., 2024). Unfortunately, very little is known about LOC100578475 nor the three neighbouring transcripts in our associated region. None have been functionally characterized, though they may conceivably regulate eye pigmentation processes and prevent omochrome deposition or production. Mutations that inhibit the ommochrome pathway lead to a white or off‐white eye colour in honey bees (Dustmann, 1968). This is unlike the red in Drosophila because honey bees do not produce drosopterins, the red pigments produced by a tandem pathway (Kim et al., 2013; Ziegler & Harmsen, 1970). Honey bees do produce pteridines in their eye cells, but pteridines apparently have no role in eye pigmentation (Dustmann, 1968). Future studies could consider targeting LOC100578475 for its functional role in eye pigmentation to help develop our understanding of the role of ncRNA in more diverse biological processes.

CONCLUSIONS

We have identified an eight‐base‐pair region where a recessive mutation likely causes white‐eye colouration in honey bees. We are unable to pinpoint the exact mechanism behind its effect and are left with only a few reasonable hypotheses. Primarily, we suspect that the novel mutation causes the misregulation of a non‐coding RNA, which impacts the synthesis or deposition of ommochromes. Further research is needed to explore this possibility.

Our discovery serves as a valuable example for the growing field of honey bee genomics. While the genome was first sequenced in 2006 (Weinstock et al., 2006), few studies have been able to make direct connections between a specific mutation and a phenotype. Beyond the scientific value, we think this case offers a tangible example that stakeholders, especially beekeepers, can relate to. We see this as a useful example to present for extension and outreach on the value of community science. Eye colour mutations have been reported for over a century, and although some basic understanding of their inheritance has been established, much remains to be learned. By working alongside beekeepers to investigate this example, we believe we have developed an insightful case study that is both scientifically valuable and practically relevant.

EXPERIMENTAL PROCEDURES

Finding mutant‐producing colonies and quantifying emergence ratios

In summer 2021 until winter 2022, we solicited the beekeeping community nationwide through social media, in‐person talks, and word‐of‐mouth with a request for information on colonies with ‘mutant drones’ (Figure 1). We were contacted by three beekeepers from Georgia, Kentucky (PP) and Minnesota who had mutant‐producing wild‐type queens. The Minnesotan queen did not survive transport to our facility in Indiana, and the Georgian queen was paint‐marked and introduced to a standard single‐deep colony and maintained until her death in summer 2023. The focal colony of this study was maintained by its owner (PP) using standard beekeeping practice except where noted below. All collection and counting were performed by PP. To generate white‐eye and wild‐type drones for sequencing (below) and to assess emergence ratios, we compelled the queen to lay drones (Büchler et al., 2024). In brief, the queen was placed on a drone foundation (frame with large‐sized cells) and she was then ‘caged’ inside a queen excluder to prevent her from leaving that comb. She was left for 24 h on the comb to lay eggs. After 24 h, the queen was freed and placed on another drone frame with the original frame left in an excluder to prevent her from laying additional eggs into it. We repeated this three additional times, creating four drone frames. We allowed the eggs to hatch, develop, and emerge over the next 24 days. On the day of emergence, the phenotype of each emerged drone was noted and then they were placed into specimen containers containing 70% ethanol and shipped to Purdue University for DNA extraction.

Sample size and whole genome sequencing

We performed both pooled and individual sequencing of white‐eye and wild‐type drones to identify the causal mutation(s). We determined the minimum sample size mathematically. As the white‐eye phenotype appears to be underpinned by a single mutation, our goal was to identify that mutational difference between wild‐type and white‐eyed haploid brothers. We first presume that there are no more than three million heterozygous sites in a diploid honey bee queen (Harpur et al., 2014). This is an extreme assumption as it represents heterozygosity in the honey bee's native range and of its most diverse population (Harpur et al., 2014). It is likely not representative of the commercial honey bee populations in the US. Regardless, under this conservative assumption, our sample size must then be the number of haploid samples in two groups (e.g., white‐eye and wild‐type) that we need to sequence such that we do not observe a mutational difference between them by chance across three million sites. If a queen is heterozygous at a locus, then there is a 50% probability that a son will inherit one or the other copy. The probability of randomly selecting N sons all with the same copy is 12N. Therefore, the probability of drawing two pools of sons, each of size N, with opposing alleles (one pool is all wild‐type and one is all white‐eyed, i.e., fixed) is, PN=2122N. If the queen has H heterozygous sites, the expected number of sites that should be fixed between pools can be written as EH=H×PN. As an example, assuming H=3,000,000 (as above), then sampling N=25 drones in each pool (50 total drones), P25=21250=1.776×10−15 and EH=5.33×10−9≈0.

Skim‐seq

Individual sequencing was performed on 50 white‐eye and 50 wild‐type drones each (100 total drones). Sequencing (Skim‐seq) was performed by Gencove. Library preparation was performed with Bead‐Linked Transposome using the Illumina Nextera DNA Flex Library Preparation kit. Libraries were quantified with fluorometry, pooled and assessed for quality with TapeStation (Agilent). Low‐pass sequencing (to ~5× depth) was performed on an Illumina NovaSeq 6000 and run as pair‐ended reads. Reads were assembled and aligned to the honey bee genome (Amel_HAv3.1; Wallberg et al., 2019) with BWA (Li & Durbin, 2009). Genotype imputation was performed by Gencove using an adapted version of the software introduced by Rubinacci et al. (2023) which applies the Li and Stephens (2003) method to estimate diploid genotype probabilities for each sample using a reference panel of phased, high‐depth genomes sourced broadly from across the honey bee native and introduced ranges. Genotypes called with less than 0.99 posterior probability (the probability that a given genotype call is correct, given the observed data and the reference panel used for imputation) were considered as missing data.

Pool‐seq

Pooled sequencing was performed on two individual groups of 60 white‐eye and 60 wild‐type drones. The thoraces of drones within each pool were ground together in a mortar and pestle and the resulting tissue homogenate was used for DNA extraction. DNA was extracted via a standard phenol/chloroform protocol (Martinson et al., 2011). Pooled samples (N = 2; 60 wild‐type in one pool; 60 white‐eye in the other) were sent for library preparation and sequencing. Samples were sequenced on an Illumina NovaSeq 6000 to 100× sequence depth, with 150 bp pair‐end reads and a 350 bp insert size (Novogene, CA). The pooled‐seq reads were aligned to the most recent honey bee genome (Amel_HAv3.1; (Wallberg et al., 2019)) using NextGenMap 0.5.0 (Sedlazeck et al., 2013). Aligned reads were then formatted, sorted and indexed with SAMTOOLS version 1.16 (Li et al., 2009).

Long‐read sequencing

DNA from an individual drone of each phenotype was sent for library preparation and long‐read sequencing (PacBio sequencing; Life Technologies, Carlsbad, CA, USA). DNA samples were quantified using a Qubit 2.0 Fluorometer. The SMRTbell library for PacBio Sequel was constructed using the SMRTbell Express Template Prep Kit 2.0 (PacBio, Menlo Park, CA, USA). We bound the library to polymerase using the Sequel II Binding Kit 2.0 (PacBio) and loaded it onto PacBio Sequel II using 8 M SMRT cells. We aligned the long‐read data (wild‐type and white‐eye) to the honey bee reference genome (Amel_HAv3.1) with BWA‐SW (Li & Durbin, 2010). We then followed up the alignment with four rounds of polishing using the short reads from wild‐type with Pilon v1.22 (Walker et al., 2014).

Genome‐wide association

All alignments were initially visualized using Integrative Genomics Viewer to ensure they appeared sound and to confirm associated sites (below) segregated fully (Robinson et al., 2011). For the 100 individually sequenced samples, the VCF file was converted to PLINK format and then run through the following pipeline: First, sites were filtered using PLINK (Purcell et al., 2007) to include only bi‐allelic SNPs that passed the base‐call confidence threshold and contained two or more minor alleles across the sample. Sites within 4000 bp of each other were arbitrarily thinned to aid in computation time and remove linked sites. Next, GCTA (Yang et al., 2011) was used to construct genomic relationship matrices and run a mixed‐linear model association between each polymorphic site and the binary mutant phenotype, accounting for genomic relatedness. After Bonferroni adjustment for multiple‐hypothesis testing, the p‐value of the association of each site to the phenotype was visualized on a Manhattan plot.

We repeated a GWA analysis using pool‐seq. Sequences were converted to mpileup file types using SAMTOOLS (Li et al., 2009). We then measured sites of differentiation between brothers and phenotypes using aligned sequences from each sample type using POPOOLATION2 (Arnold et al., 2005; Kofler et al., 2011). The fixation index (F ST) was used to measure the differentiation for every SNP (window size = 1). BED files were then sorted by chromosome and visualized in R (R version 4.2.2 through RStudio 2023) (R Core Team, 2021).

Mutation scanning

Alignments of the pool‐seq and long‐read seq (to ensure no large structural variation was also segregating) for wild‐type and white‐eye were loaded into the Integrative Genomics Viewer to identify mutations fixed in white‐eye and absent in wild‐type (Robinson et al., 2011). We identified eight fixed mutations (Table S3). Five of the eight mutations are found in previously sequenced honey bees, and therefore, unlikely to be causal to the white‐eye phenotype (NCBI EVARefSNP releases 6 and 7; Harpur et al., 2014; Sayers et al., 2019).

AUTHOR CONTRIBUTIONS

Riley R. Shultz: Formal analysis; writing – review and editing; data curation; writing – original draft; methodology; validation; investigation; project administration. Dylan K. Ryals: Writing – review and editing; formal analysis; data curation; investigation; software. Phillip Patterson: Writing – review and editing; data curation; methodology; investigation. Jonathan M. Nixon: Writing – review and editing; data curation. Izaak R. Gilchrist: Writing – review and editing; data curation. Brock A. Harpur: Writing – original draft; project administration; writing – review and editing; supervision; resources; funding acquisition; validation; methodology; investigation.

FUNDING INFORMATION

This work was funded by small donations from beekeepers, an FFAR New Investigator Award and USDA AFRI to BAH.

CONFLICT OF INTEREST STATEMENT

The authors declare that they have no competing interests.

Supporting information

Figure S1. Manhattan plot of the F _ ST _ visualizing the region of genetic differentiation between two pools of brothers (white‐eye and wild‐type) across the entire genome.

Data S1. Supporting Information.