DOI: 10.5061/dryad.gb5mkkwm5
title: Data from: Co-evolving wing spots and mating displays are genetically separable traits in Drosophila

description: 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.
description: Fly stocks The <i>D. elegans HK
 </i>(Hong Kong)<i> </i>and <i>D. gunungcola SK
 </i>(Sukarami)<i> </i>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 <i>D. elegans </i>nor <i>D. gunungcola </i>pupate on the sides of the vial, so adults were flipped out when 3<sup>rd</sup> 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 <i>D. elegans HK</i> and <i>D. gunungcola SK </i>(the same lines used previously in Yeh <i>et al</i>., 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 <i>D. elegans </i>were crossed to virgin females from <i>D. gunungcola</i>, and virgin males from <i>D. gunungcola </i>were crossed to virgin females from <i>D. elegans</i> in groups of ten males and ten females to generate fertile F<sub>1</sub> female and sterile F<sub>1</sub> 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 <i>D. elegans </i>and <i>D. gunungcola</i> to generate F<sub>1</sub> 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<sub>1</sub> hybrid female maternity, because previous work (Yeh and True, 2014) found no effect of F<sub>1</sub> hybrid maternity on trait means for wing spots and wing display in backcross populations. Multiple high-density groups of ~60 F<sub>1</sub> hybrid females were then backcrossed to ~60 virgin male <i>D. elegans </i>flies in individual vials on M food diet to create the <i>D. elegans</i> backcross recombinant population (724 individuals). To create the <i>D. gunungcola</i> backcross recombinant population (241 individuals), groups of ~60 F<sub>1</sub> hybrid females were backcrossed to ~60 virgin male <i>D. gunungcola </i>flies in individual vials on M food diet; this backcross was less successful at producing recombinant progeny than the <i>D. elegans </i>backcross direction. Behavioral assays Virgin <i>D. elegans </i>females were isolated upon eclosion, aged 10-20 days, and stored in groups of 30-40 for courtship assays. F<sub>1</sub> hybrid and recombinant backcross males were isolated individually in M food vials using CO<sub>2</sub> 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 <i>D. elegans </i>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. <i>D. elegans </i>virgin females were used in all courtship assays in case any <i>D. elegans</i> female cues were necessary to elicit male wing display behavior.<i> </i>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<i> </i>Quantification of wing spots)<i>, </i>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<sub>1</sub> hybrid and recombinant males from both backcross directions performed variable wing display behaviors during courtship as described previously (Yeh <i>et al.</i>, 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<sub>1</sub> 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 <i>D. elegans</i>, all parental male <i>D. elegans</i>, <i>D. gunungcola</i>, F<sub>1</sub> 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 <i>et al</i>., 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<sup>2</sup> using the “Measure” function. “Wing spot size” varied between 0 mm<sup>2</sup> (spotless) and 0.15 mm<sup>2</sup> between recombinant individuals. Library preparation and sequencing We estimated chromosome ancestry “genotypes” for 724 <i>D. elegans </i>backcross progeny and 241 <i>D. gunungcola </i>backcross progeny with a single Multiplexed Shotgun Genotyping (MSG) (Andolfatto <i>et al</i>., 2011) library using 965 barcoded adaptors following methods described in Cande <i>et al</i>., (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 <i>et al</i>. (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 <i>D. elegans HK </i>and <i>D. gunungcola SK </i>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 <i>D. elegans </i>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 <i>D. elegans </i>(accession number: WVIB00000000) and partially assembled scaffolds in <i>D. gunungcola </i>(accession number: WTSR00000000). Marker generation with Multiplexed Shotgun Genotyping Following methods described previously (Andolfatto <i>et al</i>., 2011; Cande <i>et al</i>., 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 <i>D. elegans </i>and <i>D. gunungcola </i>parental genomes to estimate chromosome ancestry for each backcross individual. We generated 3,425 and 3,121 markers for the <i>D. elegans </i>and <i>D. gunungcola </i>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 <i>et al</i>., 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 <i>et al</i>., 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 <i>D. elegans </i>and <i>D. gunungcola </i>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 <i>D. elegans </i>and <i>D. gunungcola </i>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 <i>D. elegans </i>or <i>D. gunungcola </i>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 <i>D. elegans </i>backcross. For three-way ANOVAs, we tested for any statistically significant interactions for max wing display angles between three different QTL peaks in the <i>D. gunungcola </i>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 <i>et al</i>., 2008) searches against the <i>D. melanogaster </i>genome (taxid: 7227) using ~10 Kbp windows of assembled <i>D. elegans </i>chromosomal regions spanning the wing spot QTL interval. Using the “GBrowse” tool on Flybase (Thurmond <i>et al</i>., 2018), we mapped regions of microsynteny to identify the orientation of each gene and exported the respective <i>D. melanogaster </i>coding region (CDS) FASTA sequences to align with the <i>D. elegans </i>X chromosome. In situ hybridization Fly genomic DNA (gDNA) was extracted from ten homogenized <i>D. elegans </i>and <i>D. gunungcola </i>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 <i>omb</i> locus in <i>D. elegans</i>: 5’-GCTGAGGATCCATTCGCTAGATTTG-3’ and 5’-GTTGTTGGAACTAGAGTTGTTGGTG-3’, and <i>D. gunungcola</i>: 5’- GCTGAGGATCCATTCGCTAGATTTG-3’ and 5’-GTTGTTGGAACTGGAGTTGTTGGTG-3’. Reverse primers were designed beginning with a T7 RNA polymerase binding sequence (TAATACGACTCACTATAG) to facilitate <i>in vitro</i> transcription. Raw PCR products were then used to generate digoxigenin-labeled RNA probes using a T7 RNA <i>in vitro</i> 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 <i>in situ </i>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 <i>D. elegans </i>and <i>D. gunungcola </i>L3 wing discs were dissected first to validate that our <i>omb </i>probes<i> </i>detected an mRNA expression pattern similar to <i>D. melanogaster </i>(Grimm and Pflugfelder, 1996; Supplementary Figure S3). Next, pupal wings were dissected at 30 and 48 h after pupal formation (APF) to probe for <i>omb </i>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, <i>in situ </i>hybridization was carried out as previously described (Vincent <i>et al</i>., 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 <i>D. elegans </i>backcross experiment, F<sub>1</sub> hybrid females were generated using the procedures described above. F<sub>1</sub> hybrid females were then backcrossed to <i>D. elegans </i>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 <i>D. elegans</i> males, and BC3 backcross males lacking wing spots were isolated to measure “Maximum wing display angles” during courtship with <i>D. elegans </i>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 <i>D. gunungcola</i> 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 <i>D. elegans</i> genomic regions across all autosomes but varied for the amount of <i>D. gunungcola </i>genome regions on the X chromosome. Introgression of black body color alleles from D. gunungcola<i> into </i>D. elegans In the <i>D. gunungcola </i>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 <i>ebony </i>gene, which appears to contribute to variation in body color between <i>D. elegans </i>and <i>D. gunungcola</i> (unpublished data). We therefore reasoned that introgressing dark body color from <i>D. gunungcola</i> into <i>D. elegans </i>would introgress the Muller Element E QTL peak underlying wing spot size differences. After six generations of backcrossing dark brown female recombinants with <i>D. elegans </i>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 <i>D. gunungcola </i>alleles linked near the Muller Element E QTL peak (Supplementary Figure S4A,B). Observing and collecting wild <i>D. gunungcola</i> and <i>D. elegans</i> in Indonesia Throughout early July 2018, <i>D. elegans </i>and <i>D. gunungcola </i>were recorded performing<i> </i>courtship in East Java, Indonesia on <i>Ipomoea sp. </i>and<i> </i> <i>Brugmansia sp.</i> 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 <i>D. gunungcola </i>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-<i> </i>Quantification of wing display behavior), to the best of our knowledge, the first recorded observations of <i>D. gunungcola </i>wing displays on flowers in the field and in the laboratory (Supplementary Figure S11; Videos 8, 9, and 11). To confirm species identification of <i>D. gunungcola </i>and <i>D. elegans </i>from the field sites mentioned above, we dissected and imaged male genitalia and compared with the laboratory strains (<i>D. gunungcola SK </i>and <i>D. elegans HK</i>) used in this study and described previously (Sultana <i>et al</i>., 1999; Kopp and True, 2002) (Supplementary Figure S9). The distal paramere [also called the pregonite (Rice <i>et al</i>., 2019)] was especially diagnostic of species identity (Supplementary Figure S9). We also performed low coverage sequencing of the new <i>D. gunungcola </i>strains’ genomes from Coban Rondo (see Methods-Library preparation and sequencing) and aligned coding sequences from the <i>omb</i> locus with the <i>D. gunungcola SK </i>lab strain (Supplementary Figure S10). A nonsynonymous coding change that distinguished the laboratory <i>D. gunungcola</i> <i>SK </i>strain<i> </i>from <i>D. elegans HK </i>also distinguished the new <i>D. gunungcola </i>Coban Rondo strain from <i>D. elegans </i>HK, matching the <i>D. gunungcola SK</i> 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 <i>post hoc</i> Tukey HSD for pairwise comparisons adjusted for multiple comparisons. See “QTL analysis” methods for statistical tests used during QTL mapping.