Understanding the genetic basis of species differences is a major goal in evolutionary biology. Pigmentation divergence between Drosophila species often involves genetic changes in pigmentation candidate genes that pattern the body and wings, but it remains unclear how these changes affect pigmentation evolution in multiple body parts between the same diverging species. Drosophila elegans and D. gunungcola show pigmentation differences in the thorax, legs, and wings, with D. elegans exhibiting male-specific wing spots and D. gunungcola lacking wing spots with intensely dark thoraces and legs. Here, we performed QTL mapping to identify the genetic architecture of these differences. We find a large effect QTL on the X chromosome for all three body parts. QTL on Muller Element E were found for thorax pigmentation in both backcrosses but were only marginally significant in one backcross for the legs and wings. Consistent with this observation, we isolated the effects of the Muller Element E QTL by introgressing D. gunungcola alleles into a D. elegans genetic background and found that D. gunungcola alleles linked near the pigmentation candidate gene ebony caused intense darkening of the thorax, minimal darkening of legs, and minimal shrinking of wing spots. D. elegans ebony mutants showed changes in pigmentation consistent with Ebony having different effects on pigmentation in different tissues. Our results suggest that multiple genes have evolved differential effects on pigmentation levels in different body regions.

Fly stocks

Drosophila elegans HK (Hong Kong) and D. gunungcola SK (Sukarami) species stocks were a gift from John True (Stony Brook University). Stock maintenance and food recipes are described in Massey et al. (2020). In brief, stocks were maintained at 23ºC on a 12 h light-dark cycle. At the third instar larval stage of development, adults were transferred onto new food, and Fisherbrand filter paper (cat# 09-790-2A) was added to the larval vials to facilitate pupariation.

Generating hybrid progeny

Male and female D. elegans HK flies will reproduce with male and female D. gunungcola SK in the laboratory to produce fertile F₁ hybrid female and sterile F₁ hybrid male offspring (Yeh et al., 2006). Creating these F₁ hybrids in populations large enough for genetic analysis, however, is difficult. Detailed methods are described in Massey et al. (2020). In brief, both D. elegans HK and D. gunungcola SK species stocks were expanded to establish populations of more than 10,000 flies. Next, virgin males and females from each species were placed in heterospecific crosses (D. elegans HK males with D. gunungcola SK females and vice versa) in groups of ten males and females to generate fertile F₁ hybrid female offspring. Dozens of these crosses were set to create ~120 F₁ hybrid female offspring. Since F₁ hybrid males are sterile (Yeh et al., 2006; Yeh and True, 2014), the F₁ hybrid females were used to generate two backcross populations for QTL mapping. Briefly, for the D. elegans HK backcross population, ~60 F₁ hybrid females were crossed in the same vial with ~60 D. elegans HK males and transferred onto new food every two weeks for ~2.5 months, resulting in 724 recombinant individuals. For the D. gunungcola SK backcross population, ~60 F₁ hybrid females were crossed in the same vial with ~60 D. gunungcola SK males and transferred onto new food every two weeks for ~2.5 months, resulting in 241 recombinant individuals.

Pigmentation quantification

For thorax pigmentation QTL analysis, male recombinants from each backcross population were organized into three thorax pigmentation classes. The lightest recombinants showed thorax pigmentation intensities similar to D. elegans HK and were given a score of 0; recombinants with intermediate thorax pigmentation intensities were given a score of 1; and recombinants with dark thorax pigmentation intensities similar to D. gunungcola SK were given a score of 2 (Figure 3A). To quantify the effects of the Muller Element E introgression region on thorax pigmentation (Figure 4B,C), individuals were placed thorax-side up on Scotch double sided sticky tape on glass microscope slides (Fisherbrand) (cat# 12-550-15) and imaged at the same exposure using a Canon EOS Rebel T6 camera mounted to a Canon MP-E 65 mm macro lens equipped with a ring light. The images were then imported into ImageJ software (version 1.50i) (Wayne Rasband, National Institutes of Health, USA; http://rsbweb.nih.gov/ij/) and converted to 32 bit grayscale. Using the “straight line segment” tool, a line was drawn between the anterior scutellar bristles (Figure 4B white dashed box) to measure the mean grayscale value of the cuticle.

For leg pigmentation QTL analysis, methods were identical to the thorax procedures above except recombinants were organized into light (0), intermediate (1), and dark (2) classes based on the pigmentation intensity of the medial side of their right hindleg femur (Figure 3B). This region of the leg was chosen because it contains few cuticular bristles that could obscure pigmentation. To quantify the effects of the Muller Element E introgression region on leg pigmentation (Figure 4B,E), the same methods as above were used except right medial hindleg images were captured. In ImageJ, the “polygon selections” tool was used to draw a polygon selection around the medial side of the right hindleg femur (Figure 4B red transparency selection) to measure the mean grayscale value of the cuticle. Data shown in Figure 4E were inverted so that higher values indicated darker pigmentation.

For wing spot pigmentation QTL analysis, methods are described in Massey et al., (2020). Briefly, right wings from male recombinants from each backcross population were imaged, and spots were quantified in ImageJ using the “polygon selections” tool to quantify wing spot size relative to wing size (Figure 4B red transparency selection). The same procedure was used to quantify the effects of the Muller Element E introgression region (Figure 4D).

Library preparation, sequencing, and genome assembly

Detailed methods for preparing and sequencing the genomic DNA (gDNA) libraries for the D. elegans HK and D. gunungcola SK backcross populations and advanced introgression line were described in Massey et al., (2020). Briefly, gDNA was extracted from male recombinants by homogenizing individuals singly in each well of a 96-well plate (Corning, cat# 3879). Each recombinant gDNA sample was then barcoded with unique adaptors, pooled into a single multiplexed sequencing library, size selected, and sequenced in a single lane of Illumina HiSeq by the Janelia Quantitative Genomics Team. Methods for assembling the D. elegans HK and D. gunungcola SK genomes to facilitate marker generation were described in Massey et al. (2020).

Marker generation with Multiplexed Shotgun Genotyping

Chromosome ancestry “genotypes” for the backcross progeny and introgression line were estimated with two Multiplexed Shotgun Genotyping (MSG) (Andolfatto et al., 2011) libraries, following methods described in Cande et al., (2012). Briefly, reads generated from the Illumina backcross sequencing library were mapped to the assembled D. elegans HK and D. gunungcola SK parental genomes to estimate chromosome ancestry for each backcross individual. We generated 3,425 and 3,121 markers for the D. elegans HK and D. gunungcola SK backcrosses, respectively, for QTL analysis [markers, phenotypes, and procedures for QTL mapping are deposited on Dryad (doi:10.5061/dryad.gb5mkkwm5)]. PDFs of chromosomal breakpoints for each recombinant are available here: https://deepblue.lib.umich.edu/data/concern/data_sets/j098zb17n?locale=en.

QTL analysis

QTL analysis was performed using R/qtl (Broman and Sen, 2009) in R for Mac version 3.3.3. (R Core Team 2018). We imported ancestry data for both backcross populations into R/qtl using a custom script (https://github.com/dstern/read_cross_msg). This script directly imports the conditional probability estimates produced by the Hidden Markov Model (HMM) of MSG (described in detail in Andolfatto et al., 2011). We performed genome scans with a single QTL model using the “scanone” function of R/qtl and Haley-Knott regression (Haley and Knott, 1992) for thorax, leg, and wing pigmentation. For QTL mapping using the D. elegans HK backcross population, we excluded 18 and 20 individuals for thorax and leg pigmentation, respectively, because fly samples were either too poor to image or sequencing reads were too shallow to map. For the D. gunungcola SK backcross population, we excluded 12 and 9 individuals for thorax and leg pigmentation, respectively, for the same reasons. For wing pigmentation QTL analysis, we previously (Massey et al., 2020) reported on an ~400 kb fine-mapped region on the X chromosome explaining the majority of variation for wing spot size. We also reported on QTL underlying variation in wing spot size with spotless recombinants removed from the analysis (Massey et al., 2020). Here, in Figure 3D, we show QTL underlying this variation to emphasize the role of Muller Element E in spot size divergence. Significance of QTL peaks at α = 0.01 was determined by performing 1000 permutations of the data.

Introgression of black body color alleles from D. gunungcola SK into D. elegans HK

To isolate individual QTL underlying body color differences between D. elegans HK and D. gunungcola SK, F₁ hybrid females and D. elegans HK backcross recombinants were generated using the methods described above. Next, dark black/brown female recombinants were repeatedly backcrossed with D. elegans HK males for four generations (BC3-BC6). Finally, we generated a single BC6 dark black homozygous introgression line that we then genotyped using MSG (Figure 4A,B).

Creating an ebony null allele in Drosophila elegans HK via CRISPR-Cas9 genome editing

Using methods described in Bassett et al. 2013, we in vitro transcribed (MEGAscript T7 Transcription Kit, Invitrogen) three single guide RNAs (sgRNAs) (Supplemental File S1) with target sequences designed based on conserved sites between D. elegans HK and D. gunungcola SK in exon 2 of the ebony gene. After transcription, sgRNAs were purified using an RNA Clean and Concentrator 5 kit (Zymo Research), eluted with nuclease-free water, and quantified using a Qubit RNA BR Assay Kit (Thermo Fisher Scientific). Next, mature (>2 weeks old) D. elegans HK males and females were transferred to 60 mm embryo lay cages (GENESEE Part Number: 59-100) on top of 60 mm grape plates (3% agar + 25% grape juice + 0.3% sucrose) at high densities (>300 flies) after brief CO₂ anesthesia. After CO₂ anesthesia, mature, mated D. elegans HK females will often dispel an embryo from their abdomen where it sticks briefly to the anus. Tapping the grape plate + embryo cage down on a hard surface 10 times causes the embryos to stick to the grape agar. Flies were then transferred back into food vials and embryos were lined up on glass cover slips taped to a glass microscope slide. For CRISPR/Cas9 injections into D. elegans HK embryos, Cas9 protein (PNA Bio #CP01), phenol red, and all three sgRNAs were mixed together at 400 ng/µl, 0.05%, and 100 ng/µl final injection concentrations, respectively. All CRISPR injections were performed in-house, using previously described methods (Miller et al, 2002). Finally, we screened for germline mutants based on body pigmentation and confirmed loss of Ebony protein by western blot (Supplemental Figure S1). To the best of our knowledge, these are the first gene editing experiments to succeed in D. elegans. We attempted the same experiments in D. gunungcola SK, but failed to recover any mutants.

Western blotting

Western blot methods were followed similar to Wittkopp et al. (2002). In brief, for each replicate per genotype, four newly eclosed (within 60 min) male flies were homogenized in 100 µl of 125 mM Tris pH 6.8, 6% SDS and centrifuged for 15 min. The supernatant was then transferred to a new protein low-bind Eppendorf tube with an equal amount of 2X Laemmli sample buffer [4% SDS, 20% glycerol, 120 mM Tris-Cl (pH 6.8)], boiled for 10 min, and stored at -80ºC. Before gel electrophoresis, samples were thawed at room temperature for 30 min, and 20 µl of each sample was loaded into individual wells of an Invitrogen NuPAGE 4-12% Bis-Tris Gel. The gel was run at 175 V for 60 min, washed, and transferred to an Invitrogen iBlot 2 PVDF Mini Stack Kit. The mini stack was run on an iBlot 2 (ThermoFisher, Catalog Number: IB21001) to perform western blotting transfer, and blocked using an Invitrogen Western Breeze anti-rabbit kit (Catalog Number: WB7106). After two washes with diH₂O, samples were incubated in 1:400 rabbit anti-Ebony (Wittkopp et al., 2002) overnight at 4ºC. Finally, samples were washed using the Western Breeze wash solution, incubated in secondary antibody solution alk-phos. (Catalog Number: WP20007) with conjugated anti-rabbit for 30 min, washed for 2 min with diH₂O, and prepared for chromogenic detection.

Statistics

Statistical tests were performed in R for Mac version 3.3.3 (R Core Team 2018). ANOVAs were performed with post hoc Tukey HSD for pairwise comparisons adjusted for multiple comparisons. See “QTL analysis” methods for statistical tests used during QTL mapping.