The evolution of sexual traits often involves correlated changes in morphology and behavior. For example, in Drosophila, divergent mating displays are often accompanied by divergent pigment patterns. To better understand how such traits co-evolve, we investigated the genetic basis of correlated divergence in wing pigmentation and mating display between the sibling species Drosophila elegans and D. gunungcola. Drosophila elegans males have an area of black pigment on their wings known as a wing spot and appear to display this spot to females by extending their wings laterally during courtship. By contrast, D. gunungcola lost both of these traits. Using Multiplexed Shotgun Genotyping (MSG), we identified a ~440 kb region on the X chromosome that behaves like a genetic switch controlling the presence or absence of male-specific wing spots. This region includes the candidate gene optomotor-blind (omb), which plays a critical role in patterning the Drosophila wing. The genetic basis of divergent wing display is more complex, with at least two loci on the X chromosome and two loci on autosomes contributing to its evolution. Introgressing the X-linked region affecting wing spot development from D. gunungcola into D. elegans reduced pigmentation in the wing spots but did not affect the wing display, indicating that these are genetically separable traits. Consistent with this observation, broader sampling of wild D. gunungcola populations confirmed the wing spot and wing display are evolving independently: some D. gunungcola males performed wing displays similar to D. elegans despite lacking wing spots. These data suggest that correlated selection pressures rather than physical linkage or pleiotropy are responsible for the coevolution of these morphological and behavioral traits. They also suggest that the change in morphology evolved prior to the change in behavior.

Fly stocks

The D. elegans HK (Hong Kong) and D. gunungcola SK (Sukarami) lines used in this study were a gift from John True (Stony Brook University). Species stocks were kept on a 12 h light-dark cycle at 23ºC on a University of Michigan “R food” diet containing molasses (http://lab-express.com/flyfoodsupplies.htm#rfood) (Wirtz and Semey, 1982). Maintaining these species on R food at high densities (50-100 flies per vial) allowed for the parental population to build up to thousands of flies to collect hundreds of virgins for interspecific crosses (see below). Neither D. elegans nor D. gunungcola pupate on the sides of the vial, so adults were flipped out when 3^rd instar L3 larvae developed and Fisherbrand filter paper (cat# 09-790-2A) was added to the food to create pupation space.

Generating hybrid progeny

Virgin males and females of D. elegans HK and D. gunungcola SK (the same lines used previously in Yeh et al., 2006; Yeh and True, 2014) were isolated upon eclosion and stored in groups of ten for one week on University of Michigan “M food”, which is the standard cornmeal diet from the Bloomington Drosophila Stock Center (https://bdsc.indiana.edu/information/recipes/bloomfood.html) with 20% higher agar content. Virgin males from D. elegans were crossed to virgin females from D. gunungcola, and virgin males from D. gunungcola were crossed to virgin females from D. elegans in groups of ten males and ten females to generate fertile F₁ female and sterile F₁ male hybrids. These crosses took ~3-4 weeks to produce hybrid progeny. The switch from R food to M food for interspecific crosses was necessary, because R food tended to accumulate condensation and bacterial growth much faster than M food when few flies occupied a vial. Since crossing D. elegans and D. gunungcola to generate F₁ hybrids tends to take several more weeks than within species crosses, the switch to M food diet allowed for maximum breeding time and the development of dozens of hybrid progeny. Once hybrid females eclosed from both interspecific cross directions, they were pooled into the same vial and aged for ten days. We did not keep track of F₁ hybrid female maternity, because previous work (Yeh and True, 2014) found no effect of F₁ hybrid maternity on trait means for wing spots and wing display in backcross populations. Multiple high-density groups of ~60 F₁ hybrid females were then backcrossed to ~60 virgin male D. elegans flies in individual vials on M food diet to create the D. elegans backcross recombinant population (724 individuals). To create the D. gunungcola backcross recombinant population (241 individuals), groups of ~60 F₁ hybrid females were backcrossed to ~60 virgin male D. gunungcola flies in individual vials on M food diet; this backcross was less successful at producing recombinant progeny than the D. elegans backcross direction.

Behavioral assays

Virgin D. elegans females were isolated upon eclosion, aged 10-20 days, and stored in groups of 30-40 for courtship assays. F₁ hybrid and recombinant backcross males were isolated individually in M food vials using CO₂ upon eclosion for at least 5 days before each courtship assay. For each assay, a single individual male was gently aspirated into a custom built 70 mm diameter bowl arena that matches the specifications in Simon and Dickinson (2010). Next, a single virgin D. elegans female was aspirated into the chamber and videotaped for the next 20 min, using a Canon VIXIA HF R500 camcorder mounted to Manfrotto (MKCOMPACTACN-BK) aluminum tripods. Videos were recorded between 09:00 and 16:00 at 23ºC. D. elegans virgin females were used in all courtship assays in case any D. elegans female cues were necessary to elicit male wing display behavior. After each assay, both the male and female were aspirated back into an M food vial and left for up to 5 days, after which each male was frozen in individual 1.5 mL Eppendorf tubes for wing spot quantification (see Quantification of wing spots), genomic DNA (gDNA) extraction, and sequencing (see Library preparation and sequencing). All courtship videos (~900 total) are available here: https://deepblue.lib.umich.edu/data/concern/data_sets/j098zb17n?locale=en.

Quantification of wing display behavior

F₁ hybrid and recombinant males from both backcross directions performed variable wing display behaviors during courtship as described previously (Yeh et al., 2006; Yeh and True, 2014). To generate quantitative measurements of wing display variation between individuals, each courtship video was played using QuickTime (version 10.4) (Apple Inc., Cupertino, CA) software in a MacOS environment and digital screenshots were manually taken for each wing display bout, defined as a bilateral wing extension performed near the female (Supplementary Figure S1). Next, for each individual fly, wing display screenshots were compared to each other to identify the maximum wing display bout per fly, defined by comparing the distance between the tips of each wing relative to the center of the fly. These maximum wing display screenshots were then imported into ImageJ software (version 1.50i) (Wayne Rasband, National Institutes of Health, USA; http://rsbweb.nih.gov/ij/) to manually measure the “Maximum wing display angle” for F₁ hybrid and recombinant males. In ImageJ, each screenshot image was inverted using the “Find Edges” function to enhance the contrast between the arena background and the edges of the fly wings (Supplementary Figure S1). Next, the “Polygon Selections” tool was used to fit an ellipse around the fly body using the “Fit Ellipse” function (Supplementary Figure S1). A Macros function (Supplementary File S1) was then used to generate major and minor axes inside the ellipse to identify the center of the fly body (Supplementary Figure S1). Finally, the “Angle Tool” was used to measure the “Maximum wing display angle” centering the vertex at the intersection of the major and minor axes and extended from wing tip to wing tip (Supplementary Figure S1). “Maximum wing display angle” varied between ~50º and ~220º between backcross recombinant individuals. 

Quantification of wing spots

Since wing spots fully form ~24 h after eclosion in D. elegans, all parental male D. elegans, D. gunungcola, F₁ hybrids, and backcross recombinants were aged at least 7 days before being frozen at -20C in 1.5 mL Eppendorf tubes. Next, using a 20 Gauge stainless steel syringe tip (Techcon) (cat# TE720100PK) the right wing of each fly was cut away from the thorax and placed on a glass microscope slide (Fisherbrand) (cat# 12-550-15) to image using either a Leica MZFLIII stereoscope equipped with a Leica DC480 microscope camera or a Canon EOS Rebel T6 camera equipped with a Canon MP-E 65 mm macro lens. Each camera was calibrated using an OMAX 0.1 mm slide micrometer to define pixel density in ImageJ software. JPEG images of wings were imported into ImageJ to measure wing spot size relative to total wing area (wing spot size / total wing area). We quantified wing spot size, rather than wing spot intensity, because we aimed to map previously identified wing spot size QTL (Yeh et al., 2006; Yeh and True 2014) and their relationship with wing display behavior. Total wing area (wing length x wing width) was approximated using length and width proxies following methods described in Yeh and True (2014). Using the “Polygon Selections” tool, the margins of black pigmentation defining each “Wing spot size” was traced and the polygon area quantified in mm² using the “Measure” function. “Wing spot size” varied between 0 mm² (spotless) and 0.15 mm² between recombinant individuals.

Library preparation and sequencing

We estimated chromosome ancestry “genotypes” for 724 D. elegans backcross progeny and 241 D. gunungcola backcross progeny with a single Multiplexed Shotgun Genotyping (MSG) (Andolfatto et al., 2011) library using 965 barcoded adaptors following methods described in Cande et al., (2012). In brief, to extract gDNA from all male backcross individuals, single flies were placed into individual wells of 96-well (Corning, cat# 3879) plates containing a single steel grinding bead in each well (Qiagen, cat# 69989). Eleven plates in total were prepared for 965 individual gDNA extractions. gDNA was isolated and purified using the solid tissue extraction procedure from a Quick-DNA 96 Kit (Zymo, cat# D3012) and a paint shaker to homogenize tissue. gDNA was tagmented using a hyperactive version of Tn5 transposase charged with annealed adaptor oligos following the methods described in Picelli et al. (2014). Unique barcoded adaptor sequences were ligated to each sample of tagmented gDNA with 14 cycles of PCR using OneTaq 2x Master Mix (NEB, cat# M0482S), and all samples were pooled into a single multiplexed sequencing library. Agencourt AMPure XP beads (Beckman Coulter, cat# A63881) were used to size select ~150-800 bp fragments and eluted in 35 uL of molecular grade water (Corning, cat# MT46000CI). The library was quantified by qPCR and sequenced in a single lane of Illumina HiSeq by the Janelia Quantitative Genomics Team.

In addition to generating the backcross sequencing library, both D. elegans HK and D. gunungcola SK parental species were sequenced at 20x coverage using an Illumina MiSeq Reagent Kit (v.3, 600 cycle PE) to facilitate genome assembly. In brief, gDNA was extracted using a Quick-DNA Microprep Kit (Zymo, cat# D4074) from 10 pooled females for each species and quantified on a Qubit 2.0 (Invitrogen). These samples were sent to the University of Michigan DNA Sequencing Core to prepare 300 bp paired-end libraries, which were quantified by qPCR and sequenced in a single lane of Illumina MiSeq. 

Genome assembly

In brief, Illumina reads from all 965 backcross recombinants were used to perform MSG on the Baylor College of Medicine D. elegans genome assembly (accession number: GCA_000224195.2). Using custom scripts in R and Python (https://github.com/masseyj/elegans), the recombination fraction between the Baylor and MSG contigs was calculated and plotted to manually tabulate joins and splits between newly assembled contigs. These new contigs were then used to assemble approximately chromosome length scaffolds in D. elegans (accession number: WVIB00000000) and partially assembled scaffolds in D. gunungcola (accession number: WTSR00000000).

Marker generation with Multiplexed Shotgun Genotyping

Following methods described previously (Andolfatto et al., 2011; Cande et al., 2012), we used the MSG software pipeline (https://github.com/JaneliaSciComp/msg/tree/master/instructions) to perform data parsing and chromosome ancestry estimation to generate markers for quantitative trait locus (QTL) analysis. In brief, using data from the Illumina backcross sequencing library (see Supplementary File S2 for the number of reads per individual), we mapped reads to the assembled D. elegans and D. gunungcola parental genomes to estimate chromosome ancestry for each backcross individual. We generated 3,425 and 3,121 markers for the D. elegans and D. gunungcola backcrosses, respectively (Supplementary Files S3, S4), for QTL analysis. 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 et al., 2003) in R for Mac version 3.3.3 (R Core Team 2018) in a MacOS environment. Ancestry data for both backcross directions were imported into R/qtl using a custom script (https://github.com/dstern/read_cross_msg), which directly imports the conditional probability estimates produced by the Hidden Markov Model (HMM) of MSG (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 “Wing spot size” and “Maximum wing display angle”. Note, for “Wing spot size”, 68 and 42 recombinants from the D. elegans and D. gunungcola backcross populations, respectively, were excluded from the QTL mapping because their wings were too damaged to quantify spot variation. Similarly, for “Maximum wing display angle”, 314 and 94 recombinants from the D. elegans and D. gunungcola backcross populations, respectively, were excluded from the QTL mapping because these males did not perform any courtship behavior during the assay. Significance of QTL peaks at α = 0.01 was determined by performing 1000 permutations of the data. Effect sizes for each QTL peak were individually estimated by comparing the mean “Wing spot size” or “Maximum wing display angle” between individuals that inherited either D. elegans or D. gunungcola alleles at each QTL peak position.

Since we detected multiple QTL peaks on separate chromosomes for “Maximum wing display angle”, we tested for the presence of epistatic interactions using two methods: First, we performed two- and three-way ANOVAs comparing the effect of each QTL peak in multiple QTL peak genetic backgrounds and found no evidence of an interaction. For two-way ANOVAs, we tested for any statistically significant interactions for max wing display angles between two different QTL peaks in the D. elegans backcross. For three-way ANOVAs, we tested for any statistically significant interactions for max wing display angles between three different QTL peaks in the D. gunungcola backcross. Second, we performed genome-wide pairwise tests using the “scantwo” function of R/qtl and Haley-Knott regression to test for non-additive interactions across all markers; LOD significance thresholds at α = 0.05, 0.01, and 0.001 were determined by performing 1000 permutations of the data for each model (Supplementary Figure S2, Supplementary Tables S1,S2).

Annotating the wing spot QTL interval

To annotate genes within the ~440 Kbp fine-mapped wing spot locus, we performed nucleotide BLAST (BLASTn) (Johnson et al., 2008) searches against the D. melanogaster genome (taxid: 7227) using ~10 Kbp windows of assembled D. elegans chromosomal regions spanning the wing spot QTL interval. Using the “GBrowse” tool on Flybase (Thurmond et al., 2018), we mapped regions of microsynteny to identify the orientation of each gene and exported the respective D. melanogaster coding region (CDS) FASTA sequences to align with the D. elegans X chromosome.

In situ hybridization

Fly genomic DNA (gDNA) was extracted from ten homogenized D. elegans and D. gunungcola females using a Quick-DNA Microprep Kit (Zymo, cat# D3021). The following forward and reverse primers were designed and synthesized by Integrated DNA Technologies (IDT) to PCR amplify 321 bp DNA templates targeting exon 5 of the omb locus in D. elegans: 5’-GCTGAGGATCCATTCGCTAGATTTG-3’ and 5’-GTTGTTGGAACTAGAGTTGTTGGTG-3’, and D. gunungcola: 5’- GCTGAGGATCCATTCGCTAGATTTG-3’ and 5’-GTTGTTGGAACTGGAGTTGTTGGTG-3’. Reverse primers were designed beginning with a T7 RNA polymerase binding sequence (TAATACGACTCACTATAG) to facilitate in vitro transcription. Raw PCR products were then used to generate digoxigenin-labeled RNA probes using a T7 RNA in vitro transcription kit (Promega / Life Technologies). RNA was ethanol precipitated and resuspended in water to analyze on a Nanodrop. Each probe was stored at -20ºC in 50% formamide before in situ hybridization.

All tissues underwent primary dissection in PBS, fixed for 30 mins in 4% PFA, washed 3X in PBT and underwent secondary dissection in PBT, were then washed 2X in MeOH, and 2X in EtOH before being stored at -20C. Male D. elegans and D. gunungcola L3 wing discs were dissected first to validate that our omb probes detected an mRNA expression pattern similar to D. melanogaster (Grimm and Pflugfelder, 1996; Supplementary Figure S3). Next, pupal wings were dissected at 30 and 48 h after pupal formation (APF) to probe for omb mRNA. To prepare pupal wings, appropriately staged pupae underwent a primary dissection: were cut in half along the anterior-posterior axis using Astra Platinum Double Edge Razor Blades, and fat body was washed out of the pupal casing using a pipette and PBS prior to fixation. After fixation, pupal wings underwent a secondary dissection to pull off the cuticle surrounding each wing and then washed using the procedure described above. Finally, in situ hybridization was carried out as previously described (Vincent et al., 2019). Briefly, we used an InsituPro VSi robot to rehydrate in PBT, fix in PBT with 4% PFA, and prehybridize in hybridization buffer for 1 hr at 65˚C. Samples were then incubated with probe for 16 h at 65˚C before washing with hybridization buffer and PBT. Samples were blocked in PBT with 1% bovine serum albumin (PBT+BSA) for 2 hours. Samples were then incubated with anti-digoxigenin Fab fragments conjugated to alkaline phosphatase (Roche) diluted 1:6000 in PBT+BSA. After additional washes, color reactions were performed by incubating samples with NBT and BCIP (Promega) until purple stain could be detected under a dissecting microscope. Samples were mounted in glycerol on microscope slides coated with poly-L-lysine and imaged at 10X magnification on a Leica DFC450C camera. 

Generating advanced recombinant introgressions on the X chromosome

To try to isolate the QTL effects for “Wing spot size” and “Maximum wing display angle” localized to the X chromosome according to the D. elegans backcross experiment, F₁ hybrid females were generated using the procedures described above. F₁ hybrid females were then backcrossed to D. elegans males, and backcross males lacking wing spots were isolated to measure “Maximum wing display angles” during courtship as described above. This procedure was repeated for seven generations to generate BC3-BC9 backcross individuals: backcross females were backcrossed en masse to D. elegans males, and BC3 backcross males lacking wing spots were isolated to measure “Maximum wing display angles” during courtship with D. elegans virgins (and so on to BC9). At each generation, an attempt was made to create stable introgression lines of advanced recombinant males lacking wing spots, but all failed to produce offspring, suggesting that D. gunungcola X-linked loci might also contain hybrid sterility factors. After seven generations of backcrossing, gDNA from all backcross males lacking wing spots was extracted and sequenced for MSG as described above. Backcross males lacking wing spots from BC4-BC9 were homozygous for D. elegans genomic regions across all autosomes but varied for the amount of D. gunungcola genome regions on the X chromosome.

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

In the D. gunungcola backcross, QTL mapping for wing spot size revealed QTL peaks linked to Muller Element C and E when spotless recombinants were excluded from the analysis (Supplementary Figure S4; Supplementary Table S3). The Muller Element E QTL peak is located near the ebony gene, which appears to contribute to variation in body color between D. elegans and D. gunungcola (unpublished data). We therefore reasoned that introgressing dark body color from D. gunungcola into D. elegans would introgress the Muller Element E QTL peak underlying wing spot size differences. After six generations of backcrossing dark brown female recombinants with D. elegans males, we crossed dark brown male and female recombinants together to create black offspring homozygous for the introgressed region. We then performed MSG on a single, dark black introgression line and found that it was homozygous for ~1.5 Mb of D. gunungcola alleles linked near the Muller Element E QTL peak (Supplementary Figure S4A,B).

Observing and collecting wild D. gunungcola and D. elegans in Indonesia

Throughout early July 2018, D. elegans and D. gunungcola were recorded performing courtship in East Java, Indonesia on Ipomoea sp. and  Brugmansia sp. flowers using Canon VIXIA HF R500 camcorders mounted to Manfrotto (MKCOMPACTACN-BK) aluminum tripods. Both species were observed in sympatry on flowers near Coban Rondo Waterfall in Batu, Batu City, East Java, Indonesia (-7.884985, 112.477311). After video recording courtship, males and females were captured using a mouth pipette and gently aspirated into glass vials containing standard fly media (glucose, corn meal, yeast extract, and agar). Isofemale lines of D. gunungcola from Bumiaji District (Batu City, East Java Province, Indonesia) were established in the laboratory on standard fly media at 24°C temperature. We quantified (see Methods- Quantification of wing display behavior), to the best of our knowledge, the first recorded observations of D. gunungcola wing displays on flowers in the field and in the laboratory (Supplementary Figure S11; Videos 8, 9, and 11). To confirm species identification of D. gunungcola and D. elegans from the field sites mentioned above, we dissected and imaged male genitalia and compared with the laboratory strains (D. gunungcola SK and D. elegans HK) used in this study and described previously (Sultana et al., 1999; Kopp and True, 2002) (Supplementary Figure S9). The distal paramere [also called the pregonite (Rice et al., 2019)] was especially diagnostic of species identity (Supplementary Figure S9). We also performed low coverage sequencing of the new D. gunungcola strains’ genomes from Coban Rondo (see Methods-Library preparation and sequencing) and aligned coding sequences from the omb locus with the D. gunungcola SK lab strain (Supplementary Figure S10). A nonsynonymous coding change that distinguished the laboratory D. gunungcola SK strain from D. elegans HK also distinguished the new D. gunungcola Coban Rondo strain from D. elegans HK, matching the D. gunungcola SK sequence (Supplementary Figure S10).   

Statistics

Statistical tests were performed in R for Mac version 3.3.3 (R Core Team 2018) using Student’s t-test (two-tailed) to test for statistically significant effects of pairwise comparisons of continuous data with normally distributed error terms. For tests comparing more than two groups, 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.