tartan underlies the evolution of Drosophila male genital morphology.

Male genital structures are among the most rapidly evolving morphological traits and are often the only features that can distinguish closely related species. This process is thought to be driven by sexual selection and may reinforce species separation. However, while the genetic bases of many phenotypic differences have been identified, we still lack knowledge about the genes underlying evolutionary differences in male genital organs and organ size more generally. The claspers (surstyli) are periphallic structures that play an important role in copulation in insects. Here, we show that divergence in clasper size and bristle number between Drosophila mauritiana and Drosophila simulans is caused by evolutionary changes in tartan (trn), which encodes a transmembrane leucine-rich repeat domain protein that mediates cell-cell interactions and affinity. There are no fixed amino acid differences in trn between D. mauritiana and D. simulans, but differences in the expression of this gene in developing genitalia suggest that cis-regulatory changes in trn underlie the evolution of clasper morphology in these species. Finally, analyses of reciprocal hemizygotes that are genetically identical, except for the species from which the functional allele of trn originates, determined that the trn allele of D. mauritiana specifies larger claspers with more bristles than the allele of D. simulans Therefore, we have identified a gene underlying evolutionary change in the size of a male genital organ, which will help to better understand not only the rapid diversification of these structures, but also the regulation and evolution of organ size more broadly.

The morphology of male genitalia can differ dramatically even between very closely related animal species (1). In Drosophila mauritiana males, for example, the size, shape, and bristle morphology of the claspers (surstyli), posterior lobes (epandrial posterior lobes), and anal plates (cerci) are strikingly different from those of its sister species Drosophila simulans and Drosophila sechellia (Fig. 1). Moreover, these differences have evolved in only the last 240,000 y since these species last shared a common ancestor (2) (Fig. 1).

Fig. 1.

Divergence in periphallic structures in the D. simulans clade and its relationship to the outgroup D. melanogaster (2). (A) Schematic representation of the male analia and external genitalia (posterior view). Posterior lobes are illustrated as dissected away on the right-hand side, in order to facilitate visualization of the claspers (outlined in orange), which are typically covered by the posterior lobes. While the shape and size of the posterior lobes is species-specific, the claspers and anal plates are very similar between D. simulans and D. sechellia, which are smaller and have less bristles than those of D. mauritiana and D. melanogaster. In addition, the clasper bristles of D. mauritiana are shorter and thicker than those of the other three species (19, 20, 48). (B–E) Scanning electron micrographs of external male genitalia (B and C) and dissected claspers (D and E) of Dmau D1 and Dsim w501, respectively. (Scale bars, 50 μm.)

As in other animal groups (1, 3–5), interspecific differences in the morphology of genital structures are thought to have been driven by sexual selection (6). However, the mechanisms [female choice, sperm competition, or sexual antagonism (5)], and their contribution to reproductive isolation between populations and species, have been difficult to address and resolve theoretically (7–9) and experimentally (10, 11). Genetic manipulation of the evolved loci would allow us to test directly the effect of male genital divergence on mating behavior and reproductive fitness and, therefore, facilitate the empirical study of these questions (12, 13). Although quantitative mapping studies of morphological differences in male genitalia between species of the D. simulans clade were first carried out more than three decades ago (14–21), the genetic bases of male genital divergence between these species has remained elusive. This is due, at least in part, to the large number of loci found to contribute to variation in size and shape of these structures (18, 19, 21).

The claspers are periphallic structures with an essential role in grasping and proprioception of the female and in securing genital coupling (12, 22–27). Previously, we found that multiple loci contribute to divergence in clasper size and bristle number between D. simulans and D. mauritiana (19). Here, we report the identification of one of these loci, tartan (trn). In addition, mapping and functional experiments strongly suggest that cis-regulatory changes in this gene underlie differences in clasper morphology between these two species.

Results and Discussion

Previously, we identified two regions on the left arm of chromosome 3 that contribute to differences in clasper size and bristle number between D. mauritiana and D. simulans (19). Here, we have generated recombinant introgression lines (ILs) between the D. mauritiana D1 (Dmau D1) and D. simulans w501 (Dsim w501) strains (Dataset S1) to increase the resolution of one of these regions, C2, from ∼3.5 Mb (24) to 177 kb. This interval explains about 16.3% of the difference in clasper size (and 37.9% of clasper bristle number) between the two parental strains (Fig. 2 and Datasets S2 and ). The claspers of lines that are homozygous for introgressed D. mauritiana DNA in C2 are significantly larger than those of natural strains of D. simulans (P < 0.001, ). The change in clasper size caused by differences in C2 is therefore outside the range of variation in clasper size in D. simulans, suggesting that C2 underlies interspecific divergence between D. mauritiana and D. simulans and not merely intraspecific polymorphism in clasper size in either or both of these species.

Fig. 2.

High-resolution mapping of differences in clasper morphology between Dsim w501 and Dmau D1. (A) Introgression line breakpoints on chromosome arm 3L define the 177-kb region C2 (gray, orange, and white boxes indicate DNA regions from Dsim w501, Dmau D1, or not verified, respectively). Coordinates are given in megabases with respect to the D. simulans genome (Flybase R2.02). This region contains eight protein coding genes including trn and is flanked by CG11279 and caps. (B) Introgression lines containing region C2 from Dmau D1 (IL43 and IL16.14) contribute 37.9% of the difference in bristles (Upper graph) and 16.3% of the clasper size difference (Lower graph) of this strain compared to Dsim w501 (Dataset S2). IL43 and IL16.14 differed significantly from IL16.30 and IL82 in clasper bristle number and in clasper area (P < 0.001). Asterisks indicate significance comparisons where P < 0.001 (Dataset S2). Shading in the pictures underneath the Lower graph indicates the area measured at the distal end of the claspers in lines containing Dsim w501 (gray) or Dmau D1 (orange) regions for C2. Boxes indicate the range, upper and lower quartiles, and median for each sample. (Scale bar, 20 μm.)

C2 contains eight protein-coding genes with orthologs in Drosophila melanogaster. RNA-Sequencing (RNA-Seq) data suggests that only one of these genes, trn, is expressed in the terminalia of D. simulans and D. mauritiana when the difference in clasper morphology develops between these two species (Dataset S3 and ). However, if the causative gene has a spatially restricted pattern of expression it may not have been detected in the RNA-Seq. Therefore, we knocked down the expression of all genes in the candidate region (with the exception of CG34429, for which there was no available UAS line) using RNA interference (RNAi) in D. melanogaster to test if these positional candidates are involved in clasper development (Dataset S4). In addition, we knocked down CG11279 and capricious (caps), a gene closely related to trn and that functionally overlaps with trn in some contexts (28–34). CG11279 and caps flank C2, but their cis-regulatory sequences may still be within this region (Fig. 2). We found that while knockdown of trn significantly reduced the size of the claspers (Dataset S4 and ), RNAi against the other nine genes tested, including caps, had no effect on clasper morphology in D. melanogaster (Dataset S4). Note that trn RNAi had no effect on the posterior lobes consistent with region C2 only affecting the claspers (Dataset S4).

trn encodes a leucine-rich repeat transmembrane protein (28, 30, 32, 33, 35, 36), and it is thought that its main function is to confer differences in affinity between cells and mediate their correct allocation to compartments in developing tissues such as the nervous system, trachea, eyes, wings, and legs (28, 30, 32, 35, 37, 38). Intriguingly, changes in trn expression can affect the allocation of cells between compartments, cause misspecification of compartmental boundaries, and even result in invasive movements of cells across such boundaries (33, 35, 38).

Our RNA-Seq data indicates that trn is more highly expressed in D. simulans during early terminalia development but is subsequently up-regulated in D. mauritiana at a later stage (Dataset S3). However, these data correspond to the sum of all of the expression domains of trn throughout the terminalia at each of these stages and may conceal more subtle localized expression differences between these species in specific tissues like the developing claspers. Therefore, we investigated the spatial pattern of trn expression throughout terminalia development using mRNA in situ hybridization (ISH) in Dmau D1 and Dsim w501 (Fig. 3 and ). Concomitantly, we observed a 4-h difference in the timing of terminalia development between the two strains used (Fig. 3 and ). We found that during early pupal stages trn is more highly expressed in Dsim w501 compared to Dmau D1 at the center of the terminalia, from where the internal genital structures will develop, which may explain the overall higher expression of trn in D. simulans at 30 h after puparium formation (hAPF) according to the RNA-Seq data (Fig. 3 and Dataset S3). However, during later stages, the expression of trn is detected in a wider domain and persists for longer at the base of the developing claspers of Dmau D1 compared to Dsim w501 (Fig. 3 ). This is consistent with higher expression of trn in D. mauritiana detected in the RNA-Seq data at ∼50 hAPF (Dataset S3). These results are also consistent with the RNAi results in D. melanogaster where knockdown of trn results in the loss of trn expression at the base of the claspers () and the development of smaller claspers (). Together, these results suggest that the higher and/or more persistent expression of the trn allele relative to the trn allele in the developing claspers is at least partially responsible for the larger claspers in D. mauritiana.

Fig. 3.

The spatial and temporal expression of trn differs in the developing claspers of D. simulans and D. mauritiana. Expression shown at 4-h intervals hAPF in Dsim w501 (A) and Dmau D1 (B). (C) Illustration of the developing structures at each morphological stage (). Black arrowheads indicate expression at the base of the developing claspers. CL, clasper; IG, internal genitalia. (D) Analysis of trn expression domain at the base of the developing clasper at stage 5. trnsim males, Dsim w501 and IL 16.30, exhibit significantly smaller expression domains than trnmau males, IL 43 and Dmau D1 (all comparisons in trn expression domain between lines are significant [P < 0.001], except for those indicated by ns (nonsignificant), see also Dataset S5). Boxes show the range, upper and lower quartiles, and the median for each sample. Representative trn expression at the base of the claspers is shown on the right-hand side. (E) The proportion of males with trn expression at the base of the clasper at 48 hAPF (between stages 5 and 6) in IL 16.30 and IL 43. IL43 males exhibit on average 51.9% more trn expression at the base of the claspers than IL 16.30 males (Dataset S5).

Quantitative analysis of trn ISH confirmed that males containing trn (Dmau D1 and IL43) exhibit a larger expression domain at the base of the developing claspers at stage 5 (50 hAPF for Dmau D1 and 46 hAPF for Dsim w501 and ILs) than those containing trn (Dsim w501 and IL16.30) (Fig. 3 and Dataset S5). Moreover, although at stage 6 IL43 and IL16.30 seem to recapitulate the pattern observed in Dsim w501 (i.e., trn expression no longer detected), we found that just before this, between stages 5 and 6 (48 hAPF in these ILs and D. simulans, ), there was variability in the presence of trn expression at the base of the developing claspers: expression was observed in 21% of IL16.30 males (i.e., males with trn) but in 74% of IL43 males (i.e., males with trn) (Fig. 3 and Dataset S5). These data further support the hypothesis that spatial and/or temporal divergence in the expression of trn underlies differences in clasper size between D. simulans and D. mauritiana.

We also carried out ISH for CG11279 and caps (which are both also expressed in the terminalia, Dataset S3) and CG34429 (which we were unable to knock down in D. melanogaster). This showed that, unlike trn, these genes are either not expressed in the developing genitalia or at least not in a pattern consistent with a role in clasper development and evolution (). For example, although caps expression in the male genitalia is generally similar to that of trn, caps transcripts were never detected at the base of the developing claspers (Fig. 3 and ).

Although our expression analyses of developing claspers suggest that cis-regulatory changes in trn are likely to contribute to differences in clasper morphology between D. mauritiana and D. simulans, we also found a total of 22 nucleotide differences in the coding sequence of trn between our mapped strains, Dmau D1 and Dsim w501. Only three of these differences are nonsynonymous and none is fixed between the two species (Dataset S6). In addition, a comparison of clasper size between strains of D. simulans and D. mauritiana with different combinations of amino acids at these three sites suggests that none of these substitutions is sufficient to explain the contribution of trn to the difference in clasper size between the species ( and Dataset S2). However, although the clasper size of the two mapped strains is also well within the range of their species ( and Dataset S2), we cannot rule out that one or more of these three amino acid substitutions may contribute to the difference in clasper size between the two strains used in this study.

To confirm that sequence divergence in trn contributes to the difference in clasper morphology between Dmau D1 and Dsim w501, we used CRISPR/Cas9 to make null alleles of D. simulans trn (in Dsim w501) and D. mauritiana trn (in IL43, see Fig. 2 and ). We then generated reciprocal hemizygotes for trn i.e., genetically identical male flies that differ only in whether they have a functional copy of trn from D. mauritiana or D. simulans (Fig. 4) (39). Comparison of the claspers between male reciprocal hemizygotes of trn shows that flies with a functional D. mauritiana trn allele have significantly larger claspers (P < 0.05) with more bristles (P < 0.001) than those with a functional D. simulans trn allele (Fig. 4 and Dataset S2). This confirms that, consistent with the effects of the introgressions containing trn (Fig. 2), D. mauritiana trn has evolved to confer larger claspers than D. simulans trn.

Fig. 4.

Reciprocal hemizygotes of trn show that this locus contributes to evolutionary differences in male clasper morphology. (A) Schematic at the top illustrates the third chromosome of the reciprocal hemizygotes carrying a functional allele of trn from only Dmau D1 (Left) and Dsim w501. We found a significant difference in their clasper area (F(3, 61) = 7.012, P < 0.001) and clasper bristle number (F(3, 83) = 26.29, P < 0.001), shown in the boxplots underneath. Flies with a functional trn allele from D. mauritiana (IL431.1 and IL431.3) have significantly larger claspers (*P < 0.05) with more bristles (***P < 0.001) than those with a functional D. simulans trn allele, w5011.1 and w5012.1 (Dataset S2). Boxes show the range, upper and lower quartiles, and the median for each sample. (B) Evolutionary changes increased the spatial domain and temporal expression of trn during clasper development in D. mauritiana have led to larger claspers with more bristles in this species compared to D. simulans. Orange and gray shading indicate broad and narrow expression of trn at the base of the developing claspers in D. mauritiana and D. simulans, respectively. The correspondingly colored clocks indicate differences in the persistence of this expression domain. (Scale bar, 20 μm.)

We have found that trn is a gene that underlies the rapid evolution in the size of a male genital organ and more generally a gene that contributes to differences in animal organ size (e.g., refs. 40–42). Many examples of phenotypic evolution, including differences in genital bristles between other Drosophila species (43), have been found to be caused by changes in the expression of transcription factors (44). However, trn is a transmembrane protein that appears to mediate differences in cell-cell contact directly through its extracellular domain, directing cells toward their correct positions via cues that are currently unknown (33, 35). Our results suggest that differences in trn expression in Drosophila are able to alter clasper size. Therefore, changes in cell affinity caused by variation in the temporal and/or spatial expression of transmembrane proteins that mediate cell interactions may represent another mechanism for the evolution of organ size. There is also some evidence, however, that trn could act as a ligand and may transduce signals, although its intracellular domain appears to be dispensable for most of its functions (28, 30, 33, 36). Therefore, further study of the function of trn and characterization of its role in organ size regulation and evolution is required.

Materials and Methods

Introgression Mapping and Phenotyping.

We increased the resolution of the previously predicted C2 region by generating recombinants between introgression line D11.01 and Dsim w501 (19) (). Flies were phenotyped and genotyped as described previously (19) using molecular markers (Dataset S7). All stocks and crosses were maintained on a standard cornmeal diet at 25 °C under a 12-h:12-h dark/light cycle unless otherwise stated. The periphalic structures and T1 legs were dissected, imaged, and measured as described in ref. 19 ( and Dataset S2). All statistical analyses were conducted in R Studio.

RNA Sequencing and Differential Expression Analysis.

Three independent biological replicates of RNA-Seq libraries were generated from abdominal tip tissue dissected from 20 to 30 males for Dsim w501 and Dmau w males at 30 and 50 hAPF (). Indexed libraries were sent to Macrogen Japan for sequencing in a single lane of HiSeq4000 (Illumina), producing 100-bp paired-end reads. Raw fastq files are deposited at DDBJ under the accession no’s. DRA006755 and DRA006758 for D. mauritiana and D. simulans, respectively. Genes were considered not expressed if reads per kilobase million (RPKM) was below 1.5. RNA-Seq analysis of genes in C2 as well as CG11279 and caps is summarized in Dataset S3.

RNAi Knockdown of C2 Candidate Genes.

We conducted RNAi knockdown of all of the genes within region C2 (with the exception CG34429 for which there was no available UAS line) in D. melanogaster using UAS-RNAi lines from Vienna Drosophila Resource Center (VDRC, www.vdrc.at) and from Bloomington Drosophila Stock Center (Dataset S8). UAS males of our candidate genes were crossed to NP6333–GAL4 driver female virgins (P(GawB)PenNP6333) (45) carrying the transgene UAS-Dicer-2 P(UAS-Dcr-2.D). Crosses for the RNAi were carried out at 25 °C. The genital morphology of the male knockdowns was compared to NP6333-GAL4, UAS-Dicer, and UAS-RNAi controls. Clasper bristle number and tibia length were measured for 16 individuals of each genotype. Differences in clasper bristle number and tibia size were assessed using a one-way ANOVA followed by a Tukey’s test (Dataset S4). For raw phenotypic data, see Dataset S2.

trn Sequence Analysis.

To assess the number of interspecific nucleotide substitutions in the coding sequence of trn, we used polymorphism data from Pool-seq data from 107 strains of D. mauritiana and from 50 strains of sub-Saharan D. simulans (46, 47) available at http://www.popoolation.at/pgt/ as well as from whole genome data for 10 strains of each species submitted to the SRA database by the University of Rochester (D. mauritiana lines: SRX135546, SRX688576, SRX688581, SRX688583, SRX688588, SRX688609, SRX688610, SRX688612, SRX688710, SRX688712; D. simulans lines: SRX497551, SRX497574, SRX497553, SRX497563, SRX497558, SRX497564, SRX497559, SRX495510, SRX495507, SRX497557). The sequence analysis is summarized in Dataset S6.

In Situ Hybridization.

We performed in situ hybridization to detect trn, CG11279, CG34429, and caps expression in the male terminalia of Dmau D1, Dsim w501, and D. melanogaster w1118 at a range of developmental time-points ().

Quantification of Temporal and Spatial trn Expression.

To investigate potential differences in the trn expression domain between introgression lines used to map the C2 interval, in situ hybridizations were carried out at 46 hAPF in Dsim w501, IL 16.30, IL 43, and 50 hAPF in Dmau D1. It is at these time points that the largest differences in trn expression can be detected between the two parental species (Fig. 3 ). Dsim w501, IL 16.30, and IL 43 are morphologically equivalent at these stages ().

Generation of Reciprocal Hemizygotes and Statistical Analysis.

We inserted 3xP3-DsRed to disrupt the reading frame of trn in Dsim w501 and IL43 using CRISPR/Cas9 (). Injections were carried by The University of Cambridge Department of Genetics Fly Facility. Transgenic Dsim w501 and IL43 males heterozygous for the mutation were then crossed to noninjected IL43 and D. simulans w501 virgin females, respectively, to generate F1 males carrying the mutation (i.e., hemizygous for trn allele).

See extended methodological details in .