Genetic architecture and functional characterization of genes underlying the rapid diversification of male external genitalia between Drosophila simulans and Drosophila mauritiana.

Male sexual characters are often among the first traits to diverge between closely related species and identifying the genetic basis of such changes can contribute to our understanding of their evolutionary history. However, little is known about the genetic architecture or the specific genes underlying the evolution of male genitalia. The morphology of the claspers, posterior lobes, and anal plates exhibit striking differences between Drosophila mauritiana and D. simulans. Using QTL and introgression-based high-resolution mapping, we identified several small regions on chromosome arms 3L and 3R that contribute to differences in these traits. However, we found that the loci underlying the evolution of clasper differences between these two species are independent from those that contribute to posterior lobe and anal plate divergence. Furthermore, while most of the loci affect each trait in the same direction and act additively, we also found evidence for epistasis between loci for clasper bristle number. In addition, we conducted an RNAi screen in D. melanogaster to investigate if positional and expression candidate genes located on chromosome 3L, are also involved in genital development. We found that six of these genes, including components of Wnt signaling and male-specific lethal 3 (msl3), regulate the development of genital traits consistent with the effects of the introgressed regions where they are located and that thus represent promising candidate genes for the evolution these traits.

MALE genitalia exhibit striking diversity in morphology and often evolve faster than other phenotypic traits between groups of animals with internal fertilization (Eberhard 1985). To understand the rapid evolution of these traits it is necessary to characterize the genetic architecture underlying differences in genital morphology. This requires assessment of the number of loci involved, how they interact, and the extent to which they affect other traits. For example, studies of quantitative traits suggest that epistasis is pervasive (reviewed in Mackay 2013), and there is experimental evidence that this may constrain the number of possible evolutionary paths that lead to trait variation (Weinreich ; Blount ; Bridgham ; Gong ; Park and Lehner 2013). However it has been shown that in certain theoretical scenarios, epistasis may actually facilitate evolution (Wagner ). Similarly, while in general, pleiotropic mutations are thought to contribute less to morphological change than mutations with more restricted effects (Stern and Orgogozo 2008; Wagner and Zhang 2011; Nunes ), it has been argued that high pleiotropy may not always limit evolution (Wang ; Hill and Zhang 2012; Zhang 2012).

Identifying the loci underlying differences in genital morphology between species will ultimately also allow us to determine whether they are dominant or recessive, and infer if they evolved preferentially through selection on standing variation, by selection on new mutations, or through a combination of both. In most evolutionary scenarios the chance and speed of fixation of a beneficial allele is greater if the allele was already present in the population than if it arose through a new mutation (Hermisson and Pennings 2005; Barrett and Schluter 2008). However, the frequency of new beneficial mutations can also increase very quickly if the mutations are dominant and the selection coefficient and/or the population is large enough.

Here, we address these questions by investigating differences in the male terminalia of two species of the Drosophila simulans clade, D. simulans and D. mauritiana. In Drosophila, genital and anal structures often vary dramatically among males of closely related species (Rizki 1951; Coyne 1983; Liu ; True ; Lachaise ; Zeng ; Kopp and True 2002; Markow and O’Grady 2005; Soto ; Richmond ; Lang ), as well as more subtly within species (Lachaise ; Soto , 2013; Andrade ; McNeil ; Schafer ).

The male analia develops from abdominal segment A10 and consists of two anal plates that are positioned either side of the anus (Sanchez and Guerrero 2001) (Figure 1A). The male genital arch develops from segment A9 and is a horseshoe shaped structure composed of two bilateral pairs of ventral projections, the lateral plate (epandrium) and the clasper (surstylus) (Sanchez and Guerrero 2001) (Figure 1A). Within the D. simulans species clade D. mauritiana males have claspers and anal plates that are larger than those of D. simulans and D. sechellia and that carry more bristles, which are generally shorter and thicker (Figure 1, B, C, F, and G and Tsacas and David 1974; True ). In addition, the posterior lobes, which are extensions of the lateral plate, have diverged in shape and size in all species of the group (Figure 1, D and E and Tsacas and David 1974; True ).

Figure 1

Drosophila male external genitalia and analia. (A) SEM image of a D. mauritiana male genital arch highlighting the right clasper (Cl), posterior lobe (PL), and anal plate (AP). Morphological divergence between D. mauritiana and D. simulans is shown, respectively, in B and C for claspers, D and E for posterior lobes, and F and G for anal plates. Bars, 100 μm (in A) and 50 μm (in B–G).

Studies of inter- and intraspecific copulation anatomy in Drosophila suggest that claspers, posterior lobes, and anal plates play an important role during copulation (Robertson 1988; Acebes ; Jagadeeshan and Singh 2006; Kamimura and Mitsumoto 2011; Yassin and Orgogozo 2013). Therefore, differences in the morphology of these traits may contribute to the delayed and incomplete genital coupling observed in interspecific crosses in comparison to conspecific crosses (Jagadeeshan and Singh 2006).

Quantitative trait locus (QTL) mapping experiments of posterior lobe differences between species of the D. simulans clade have revealed that they have a polygenic and predominantly additive basis (Liu ; Laurie ; True ; Macdonald and Goldstein 1999; Zeng ). Zeng previously found a minimum of 19 QTL underlying the difference in posterior lobe morphology between D. mauritiana and D. simulans, and nearly all of these loci acted in the same direction, suggesting a history of strong directional selection. True mapped the genetic basis of anal plates and clasper differences between D. mauritiana and D. simulans and identified several QTL for both traits. For clasper morphology they found a single QTL on chromosome 3L that explained ∼50% of the clasper bristle number difference between D. mauritiana and D. simulans (True ); however, the resolution achieved was only 20 cM. Furthermore, although it is not known if the same region affects clasper size because these authors did not map this trait, this region does overlap with QTL for both posterior lobe and anal plate traits (True ; Zeng ), which may suggest that the same or closely linked loci could underlie coordinated changes in genital morphology and therefore help explain their rapid diversification.

We carried out QTL mapping of clasper size and bristle number between D. mauritiana and D. simulans. We then used an introgression-based approach to investigate the causative regions on chromosome 3 at higher resolution. We found that all loci contributing to clasper size divergence also affected clasper bristle number but that the genetic basis of clasper divergence is independent from the genes underlying posterior lobe and anal plate evolution. Consistent with previous reports (Long ; True ; Gurganus ), we found evidence for epistasis between loci affecting clasper bristle number but not for cuticle size. In addition, we found that most D. mauritiana alleles are dominant with respect to the D. simulans alleles in regions affecting the claspers and recessive in regions affecting posterior lobes, suggesting that these two traits have experienced different evolutionary histories.

To then investigate which positional candidate genes in our introgressed regions on chromosome 3L are required for genitalia development, we conducted an RNAi screen in D. melanogaster. We tested 56 genes with reported differences in expression between species (Masly ) or sexes (Chatterjee ) and found that only 6 of these genes had an effect on the morphology of genital traits, which was consistent with the effect of the introgression lines and the expression differences. None of these genes was previously known to be involved in the regulation of genital development but our study shows that they are promising candidates for the divergence of the claspers or posterior lobes between D. simulans and D. mauritiana.

Materials and Methods

Fly strains and mapping populations

For QTL mapping, we generated a backcross population of 244 males as described in Arif , by crossing virgin D. simulans yellow (), vermillion (), (), (hereafter YVF) females (Drosophila Species Stock Center, DSSC, University of California, San Diego, Stock no.14021-0251.146), to D. mauritiana TAM16 males and backcrossing the resultant F1 virgin females to D. mauritiana TAM16 males.

The X chromosome region between and was introgressed from D. mauritiana TAM16 into D. simulans YVF by crossing virgin D. simulans YVF to D. mauritiana TAM16 males, followed by 16 generations of backcrossing of y virgin females to D. simulans YVF males (Supporting Information, Table S1). Using a similar crossing scheme, we introgressed the region on 3L, between markers D1 and Q1 (cytological bands 66B and 79F, respectively; True ), from D. mauritiana w into D. simulans w501 by backcrossing for up to 23 generations (Table S1). D1 and Q1 are P-element insertions that carry and therefore rescue eye pigmentation in D. mauritiana w flies. The Q1 insertion rescues eye pigmentation more strongly than D1, allowing us to distinguish between flies that carry only D1 from flies that carry only Q1, and both from flies carrying D1 and Q1. D1 and Q1 introgression lines are hereafter called D and Q lines. In both cases, double recombinant flies were identified and selected against in every new generation of introgression using molecular markers (Table S2). All stocks and crosses were maintained on a standard cornmeal diet at 25° under a 12-h:12-h dark/light cycle unless otherwise stated.

Phenotypic data

We excised the terminalia from each experimental male and mounted dissected claspers, posterior lobes, and anal plates in Hoyer’s medium. When possible, both left and right structures were collected and the average value was used for subsequent analysis. All images were captured using a Leica M205 stereomicroscope at ×400 magnification and a DFC300 camera. Cuticle area measurements for all three traits are shown in Figure S1. Areas were measured using either Adobe Photoshop CS5 or ImageJ. Clasper and anal plate bristle numbers were counted directly under the stereomicroscope. We also mounted and measured T1 tibias from each male to use as a proxy for body size. However, since genitalia are hypoallometric relative to the rest of the body (Coyne ; Liu ; Macdonald and Goldstein 1999; Eberhard 2009; Shingleton ; Masly ) and none of our conclusions were affected by analysis performed with a body size correction (most of the introgression lines did not differ in T1 tibia length, Table S1), we chose to present our phenotypic measurements of genitalia traits without body size correction.

We performed principal component analysis (PCA) of posterior lobe shape using singular value decomposition of elliptical Fourier coefficients calculated from posterior lobe outlines standardized for size, orientation, location, and handedness, as described in Masly (Table S3).

Scanning electron micrographs (SEMs) were taken on a Hitachi S3400N after dehydrating the terminalia excised from specimens in a graded ethanol series and critical point dried in CO2 (Tousimis Samdri-780). The dried samples were mounted on stubs that were then sputter coated with gold (Polaron coater).

Genotypic data

DNA was extracted from fly bodies, after dissection of the terminalia (see above), using a high-salt extraction protocol (Miller ). The markers used for the QTL mapping are described in Arif and the genotypes of the backcross population are given in Table S4. A mix of visible (, , and on the X chromosome and D1 and Q1 on 3L), restriction fragment length polymorphisms (RFLPs), insertion/deletion (indel) polymorphisms, and SNPs were used for the introgression mapping. Marker location and primer sequences are shown in Table S2.

We also sequenced the genomes of these introgression lines to improve the resolution and accuracy of recombination breakpoints in the D and Q introgression lines and to confirm that these regions have been introgressed from D. mauritiana to D. simulans without introgressing other regions of the genome that may also have an effect on the genitalia traits studied. The DNA of 10 pooled individuals (five males and five females) from each homozygous introgression line was extracted using the high-salt extraction protocol (Miller ). Each sample of DNA was ligated to a different barcoded adaptor sequence. All 40 samples were pooled, prepared following Cande , and sequenced in a single lane of Illumina HiSeq. Reads were mapped to the D. simulans w501 genome (Hu ), ancestries of chromosome regions were estimated with a hidden Markov model, and recombination breakpoints calculated using the multiplexed shotgun genotyping (MSG) pipeline (Andolfatto ). Files with breakpoint data are available at Dryad:doi:10.5061/dryad.c194d.

QTL mapping

QTL analyses were performed in R/qtl (Broman ; R Development Core Team 2012). The genetic map was constructed using default parameters and the Kosambi map function with a total of 34 markers on the 2nd, 3rd, and X chromosomes. We performed QTL scans using Haley–Knot (HK) regression, which is an approximation to standard interval mapping (Haley and Knott 1992), and determined genome-wide statistical significance thresholds (α = 0.05 and 0.005) for each phenotype (clasper bristle number and size) using 1000 permutations. We tested for any possible nonadditive interactions between QTL by fitting full linear models with all QTL significant at α ≤ 0.05 and all possible interactions between them, using the fitqtl function of R/qtl. To estimate effect sizes, we fit a new model with only those QTL and interactions that remained significant at α ≤ 0.05 and calculated their position (optimized using maximum likelihood) and 2-LOD support intervals. Furthermore, we estimated additive allelic effects of all significant QTL in three ways: (1) the difference in phenotype averages between homozygotes (D. mau/D. mau) and heterozygotes (D. mau/D. sim) for autosomes and hemizygotes for the X chromosome (labeled effect size in Table 1); (2) differences in no. 1 (above) standardized by half the difference between parental lines for autosomes and the entire difference in the case of the X chromosome (labeled relative effect size in Table 1); and finally (3) the percentage of phenotypic variation accounted for by the significant QTL in the backcross population (labeled variance explained in Table 1).

Table 1

QTL for clasper bristle number and size

			2 LOD support region (cM)		2 LOD support region (Chr@Mb)a		Additive allelic effectsb
Trait	Peak location (Chr@cM)	Peak significance (LOD)	From	To	From	To	Effect size (SE)	Relative effect size (%)	Variance explained (%)
Clasper bristle number	X@22.0	3.632	0	44	X@0	X@5.4	−0.56 (0.18)	11.2	6.42
	X@77.0	7.141	59	95	X@6.9	X@11.7	−0.91 (0.21)	18.2	13.14
	3L@49.0	4.283	32	99	3L@5.7	3R@1.9	−0.47 (0.17)	18.8	7.63
	X@22.0:X@77.0	1.868	—	—	—	—	−1.29 (0.44)	—	3.24
	X@77.0:3@49.0	2.389	—	—	—	—	1.44 (0.44)	—	4.16
Clasper size	2R@94.5	3.408	10	113	2L@1.2	2R@4.4	−155.76 (39.05)	17.4	5.84
	3L@56.0	2.408	9	111	3L@1.2	3R@8.1	−99.58 (40.01)	11.14	4.08
	3R@112.0	6.325	96	133	3R@1.9	3R@15.3	−208.98 (46.00)	23.4	11.20
	3L@56.0:3R@112.0	1.118	—	—	—	—	−214.81 (95.23)	—	1.87

Regions corresponding to the 2-LOD support interval based on the next marker position closest to the interval boundaries.

See main text (Materials and Methods) for details on how these different measures of effect size were calculated.

Statistical analysis

To determine if the X chromosome introgression lines had an effect on clasper bristle number, we performed a one-way ANOVA between genotypes (yf and yvf as well as vf, for introgression lines X:13.25 and X:13.27), and a two-way ANOVA between genotypes per introgression line. To determine if the 3rd chromosome introgression lines had an effect on different genitalia traits, we compared D and Q introgressions to D. simulans w501 using Dunnett’s test. To identify significant differences in effect size between introgression lines, we conducted pairwise comparisons between D introgression lines using Tukey’s honest significant difference test in conjunction with a one-way ANOVA (Table S5).

To determine the dominance relationship(s) between D. mauritiana and D. simulans alleles in the introgressed regions, we compared absolute differences between D. simulans w501 and heterozygous introgression lines (MS-SS) to the absolute differences between heterozygous and homozygous introgression lines (MM-MS) with a Wilcoxon rank-sum test. The significance (at α ≤ 0.01) of the W statistic was evaluated by random permutation, with replacement, of trait values (n = 2000). If MM-MS is not significantly different from MS-SS, D. mauritiana and D. simulans alleles are co-dominant. D. mauritiana alleles are dominant or recessive over D. simulans alleles if the difference between MM-MS is smaller or larger than MS-SS, respectively (Table S6). To test for epistatic interactions between D and Q introgressed regions, we conducted a two-way ANOVA of the clasper and posterior lobe size of w501, and D-introgression/w501, Q-introgression/w501 and D-introgression/Q-introgression heterozygous males (Table S7). All basic statistics were implemented using R (R Development Core Team 2012).

RNAi screen of candidate genes

We conducted an RNAi screen in D. melanogaster to determine if positional candidate genes that are differentially expressed either between sexes (Chatterjee ) or between species (Masly ) are required for genitalia development (Table S8). UAS-RNAi lines (Vienna Drosophila RNAi Center) for genes located in all 3L regions affecting the morphology of posterior lobes and the smallest region affecting clasper morphology were crossed to the NP6333-Gal4 driver () (Chatterjee ) also carrying a UAS-Dicer-2 transgene insertion to increase RNAi activity (Dietzl ). NP6333-Gal4 drives expression in the wing, leg, and genital discs (Stieper ) and was kindly provided by Mark Siegal (Chatterjee ). We determined the effect of knocking down positional candidate genes by comparing the genital morphology (clasper bristle number, posterior lobe size, and shape) of these flies to that of NP6333-Gal4; UAS-Dicer driver males and UAS-RNAi controls. RNAi flies and controls were reared at either 25° or 28°.

Results

The claspers of D. mauritiana males are more than twice as large and carry 25–30% more bristles than those of D. simulans (True ) (Figure 1 and Table S1). To identify the genetic loci underlying these differences and to determine if the same loci contribute to variation in clasper size and bristle number, we carried out QTL mapping of both traits using progeny from a D. simulans/D. mauritiana F1 backcross to D. mauritiana.

QTL analysis of clasper size and bristle number

Clasper bristle number is moderately positively correlated with clasper size in our mapping population (Pearson’s r = 0.24, P < 0.001), suggesting a partially common genetic basis for both traits. Interval mapping of clasper size using Haley–Knot regression identified two significant QTL on chromosome 2 (at genome-wide P < 0.05) and a further three QTL on chromosome 3 (Figure 2). After testing the fit of a full linear model that included all significant QTL and their interactions, we retained only those QTL that remained significant at P < 0.05 (see Materials and Methods). We found that clasper size variation between D. simulans and D. mauritiana is best explained by a minimum of three QTL, one at 94 cM on the 2nd chromosome and two on the 3rd chromosome at 56 cM and 112 cM. We also found a significant less-than-additive (i.e., negative) interaction between the two QTL on the 3rd chromosome (Table 1 and Figure S2), where the effect of the double homozygote is smaller than expected under simple additivity. The full model accounted for a total of 26% of the phenotypic variation in the backcross population.

Figure 2

QTL affecting clasper bristle number and size. LOD profiles from a HK regression analysis of clasper bristle number (gray curves) and clasper size (black curves). The 5% significance threshold LOD was the same for both traits (gray horizontal dashed line). The dotted lines represent the 0.5% significance threshold for clasper bristle number (gray) and clasper size (black). Ticks above and below the x-axis represent position of markers and distance (in centimorgans), respectively.

For clasper bristle number differences, we found two significant QTL on the X chromosome as well as two marginally significant peaks on 3L (Figure 2). After testing the fit of the full linear model, the estimated positions of the QTL on the X chromosome are 22 cM and 77 cM (Table 1). While the QTL at the proximal end of chromosome 3 was excluded after testing the fit of the full linear model, the QTL at 49 cM remained significant and shows a significant positive interaction with the largest effect QTL at 77 cM on the X chromosome (Table 1). In fact the effect of the D. mauritiana allele at 49 cM can only be detected if the X chromosome is D. mauritiana at 77 cM (Figure S2). The full linear model, which also includes an interaction between the two X chromosome QTL (Table 1 and Figure S2), explains 25% of the phenotypic variation in the backcross population. Given the proximity of the QTL at 49 cM for clasper bristle number and the QTL at 56 cM for clasper size, and the large 2-LOD confidence intervals of these QTL, it is possible that these effects could be caused by the same locus. Indeed, the effect of the 3L locus at 49 cM on clasper size is very similar to the effect of the 3L locus at 56 cM (not shown).

Our QTL analysis shows that clasper size and clasper bristle number differences between D. simulans and D. mauritiana are each explained by a minimum of three QTL and that one of these regions may contribute to both traits. In addition, we also found nonadditive interactions between QTL for both of these genital traits.

X chromosome introgression lines have no effect on clasper bristle number

To verify the results of our QTL analysis, we generated six introgression lines containing the X chromosomal regions underlying QTL for clasper bristle number identified above. After 16 generations of backcrossing to D. simulans, we collected male progeny and compared clasper bristle number of yf males (and vf and for introgression lines X:13.25 and X:13.27) with D. mauritiana DNA introgressed on the X chromosome (Figure S3) to that of their yvf sibling males from the same introgression line (i.e., without D. mauritiana DNA). Based on the estimated QTL effects (Table 1), standard deviation of the parental D. simulans strain, and sample size of introgression lines, our power to detect the major X chromosome QTL is only 60% for introgression lines X9.18 and X13.27, but is at least 80% for the remaining lines. Surprisingly, however, we found that none of the introgression lines has an effect on clasper bristle number (Table S1). One possible explanation for these results is that our statistical power is overestimated due to upward biases in the estimation of QTL effect sizes (Beavis 1998; Rockman 2012). Alternatively, the X chromosome effect may be dependent on the genetic background; since the QTL mapping population and introgression lines were generated by backcrossing in different parental directions, the X chromosome introgression lines lack the D. mauritiana 3rd chromosome QTL at 49 cM that interacts with the X chromosome QTL of large effect (Figure S2). We tried to verify this interaction by crossing X and 3rd chromosome introgression lines. However, perhaps expectedly, given the greater density of hybrid male sterility loci on the X chromosome than on the autosomes (Coyne and Orr 1989; True ; Tao ; Masly and Presgraves 2007), we were unable to generate homozygous X chromosome introgressions from the lines used in this study because yf males are sterile. We were also unable to obtain viable male progeny from crosses between the different X chromosome introgressions and introgressions on chromosome 3 (see below). These fertility and viability effects are probably caused by the disruption of the Winters sex-ratio system (Tao ; Kingan ) in these flies.

At least three regions on the 3rd chromosome contribute to divergence in clasper morphology

We used a D. mauritiana w strain carrying the visible markers D1 (at ∼7.5 Mb on 3L) and Q1 (on 3L near the centromere) (True ) to introgress the QTL regions affecting clasper size on 3L and 3R, respectively into D. simulans w501 (Figure 3A and Figure S1). After at least eight generations of backcrossing, all introgression lines exhibited an effect on clasper size (when compared to D. simulans w501, Figure 3B and Table S1). Interestingly, all the lines also have significantly more clasper bristles than D. simulans w501, suggesting that the same regions contribute to both clasper size and bristle number differences between these two species (Figure 3C and Table S1).

Figure 3

Mapping resolution and effects of the 3rd chromosome introgression lines. (A) The two arms of chromosome 3 are shown with D. simulans coordinates in megabases. Open triangles indicate the position of the visible markers D1 (on the left) and Q1 (on the right). The centromere position is indicated by an open oval. For each introgression line, the colors indicate DNA from D. mauritiana (solid), D. simulans (open) and either D. mauritiana or D. simulans (shaded). (I) All D and Q introgression lines have larger clasper size (B) and clasper bristle number (C) than D. simulans w501, indicating the contribution of regions C1 and C3 to the difference between species. Since introgression line D11.01 has a significantly larger effect than the other D introgression lines (Table S5) an additional region, C2, must also contain at least one gene contributing to variation in both traits. (II) Significant differences in effect size between D introgression lines (D21.43 and D21.43e: no effect, D20.32/D23.42: ∼7%, D08.04/D11.03/D20.37: 14–17%, D11.01: ∼22%, Table S1 and Table S4) facilitate mapping of posterior lobe size variation (D) to at least three regions on 3L, P1 (for D08.04, D11.03 and probably D11.01), P2 (for all D introgression lines except D21.43), and P3 (for D20.37 and D11.01). Posterior lobe shape (E) maps only to P1 and P2 because D20.37 does not differ significantly from D20.32. One additional region on chromosome arm 3R (P4) also contributes to posterior lobe size and shape differences (Q08.05 has an effect but Q08.02 does not). (III) Mapping resolution of anal plate size (F). Among the D introgression lines, only five (D08.04, D11.01, D11.03, D20.37, and Q08.05) show significantly smaller anal plate sizes than D. simulans w501. Therefore, anal plate size variation between species maps to three regions, A1–A3. For each region the black and gray bars indicate minimal and maximum regions, respectively, from known mapping resolution. Whole genome sequencing confirmed that there is no D. mauritiana DNA elsewhere in the genome of any of our introgression lines. Asterisks represent the significance level (Dunnett’s test comparing with w501): *P > 0.05; **P > 0.01; ***P > 0.001.

Comparative analysis of D and Q introgression lines allowed us to map the differences in clasper size and bristle number to one region, C1, common to the D introgressions and one region, C3, common to the Q introgressions (Figure 3A). Moreover, since introgression D11.01 has a significantly larger effect on clasper morphology than the other D lines (Table S5), this implies that there is an additional region, C2, on 3L (Figure 3A). While this analysis relies on the assumption that the effects of the loci are additive, whereas the observations could be explained by more complicated genetics, it allowed us to delimit minimal regions that should contain at least one gene underlying the trait difference.

Region C1 explains 8 and 13% of the interspecific difference in clasper size and bristle number, respectively. This region is <500 kb and harbors a maximum of 52 genes (with D. melanogaster orthologs, see Table 2). Region C2 is up to 3.5 Mb and explains 16 and 20% of the interspecific differences in clasper size and bristle number, respectively. Region C3 could be as large as 5.5 Mb and can explain 12 and 38% of the interspecific differences in clasper size and bristle number, respectively (Table 2).

Table 2

Effect of the 3rd chromosome introgressions on genital morphology

Region	Relative effect size (%)a				Maximum size (bp)	Coordinates Dsim R1.4 (FlyBase)	Maximum number of D. sim (D. mel orthologs)	Maximum number of sex biased genesb	Maximum number of D. mau biased genesc	Candidates common in b and c
	CS	CBS	PLS	APS
C1	8	13	0	0	446,500	3L:7905718.. 8352218	74 (52)	3	12	1
C2	16	20	0	0	3,486,822	3L:11304252..14791074	399 (312)	14	55	4
C3	12	38	0	0	5,448,709	3L:2059152..3R:3040233	693 (550)	10	96	0
P1/A1	0	0	9	−63	1,833,957	3L:5494884..7328841	236 (209)	5	30	3
P2	0	0	6	0	1,947,968	3L:8670985..10618953	251 (207)	2	31	0
P3/A2	0	0	7	−65	561,618	3L:10745193..11306811	92 (78)	1	13	1
P4/A3	0	0	13	−24	825,209	3R:3034107..3859316	98 (82)	4	9	2

Effect of homozygous introgressions as a percentage of the difference between parental strains.

Chatterjee et al. 2011.

Masly .

Clasper regions have no pleiotropic effects on the morphology of posterior lobes or anal plates

Although anal and genital structures develop from different compartments within the genital disc (Sanchez and Guerrero 2001), it has been shown in D. melanogaster that the development of these two primordia is coordinated (Estrada ; Gorfinkiel ). Hence, it is possible that genetic variation leading to clasper size diversification also affects the size of other traits such as the posterior lobes and anal plates. Indeed, the 3L region we have identified that affects clasper size and clasper bristle number broadly overlaps with previously identified QTL for posterior lobe and anal plate differences between D. simulans and D. mauritiana (True ). Therefore, we investigated the effects of D and Q introgression lines on these two traits and found that all but three introgression lines, D21.43, D21.43e (generated from D21.43 by an additional five generations of backcrossing to D. simulans w501) and Q08.05 have significantly smaller posterior lobes than D. simulans w501 (Figure 3D and Table S1), i.e., the D. mauritiana alleles move the average phenotype values in the direction of the average D. mauritiana parental phenotype average. Given the differences in effect size between lines (Figure 3D and Table S5), and again assuming the loci act additively, we were able to map interspecific variation in posterior lobe size to a minimum of three regions on chromosome arm 3L (P1–P3) and one region on 3R (P4) (Figure 3A). While P1 and P2 are each up to 2 Mb in length, P3 and P4 are only 560 kb and 825 kb, respectively (Table 2).

In a previous study of differences in posterior lobe morphology between D. mauritiana and D. sechellia, Masly found that none of the introgression lines had a consistent effect on both posterior lobe size and shape with respect to the parental difference. Therefore, we performed PCA of Fourier coefficients for each posterior lobe to determine if any of our introgression lines have a significant change in posterior lobe shape relative to D. simulans w501 and in the direction of D. mauritiana. While the first three principal components (PC1–PC3) explain >90% of the variation in the dataset (Table S3), only PC1 (which explains 69% of the variation) differs significantly between introgression lines and D. simulans w501 (Figure 3E and Table S8). All introgression lines with an effect on posterior lobe size also differ significantly in posterior lobe shape relative to D. simulans w501 and all show a shape change toward the phenotype of D. mauritiana (Figure 3E and Table S8). However, introgression line D20.37 does not differ significantly from D20.32; therefore, posterior lobe shape differences map to at least three regions, P1, P2, and P4 (Figure 3A).

Anal plate size maps to two regions that also affect both posterior lobe size and shape, P1 and P4, and one region (P3) that affects only posterior lobe size (Figure 3, A and F and Table S1). However, contrary to their effects on posterior lobe morphology, the D and Q introgression lines affect anal plate size in the opposite direction of the D. mauritiana parental average.

In contrast to their effects on cuticle area, none of the introgression lines affect the number of anal plate bristles with respect to D. simulans w501 (Figure 3G and Table S1). This suggests not only that anal plate bristle number is unaffected by the regions responsible for anal plate cuticle formation, but also that clasper bristle number can evolve independently from anal plate bristle number.

Overall, regions that contribute to interspecific differences in clasper morphology exhibit no pleiotropic effects on either anal plates or posterior lobes. While, one posterior lobe region, P2, does not affect the anal plates, generally all regions underlying posterior lobe differences also affect anal plate size, suggesting that the underlying causative loci affect both these traits.

Allelic interactions within and between loci

Intralocus allelic interactions can affect the rate of evolution at a given locus (Hartl and Ruvolo 2012). To determine the type of genetic interaction between D. mauritiana and D. simulans alleles within each region underlying differences in analia/genitalia traits, we compared absolute differences between D. simulans w501 and heterozygous D or Q introgression lines (MS-SS) to the absolute differences between heterozygous and homozygous D or Q introgression lines (MM-MS).

Our results suggest that for clasper size, D. mauritiana alleles are partially dominant to D. simulans alleles in all three regions (C1, C2, and C3) that affect clasper size (Table S6). In contrast, for anal plate size and posterior lobe morphology (size and shape), we found that D. simulans alleles in the different regions on chromosome 3L are partially dominant to D. mauritiana alleles. Finally, in region P4/A3, D. mauritiana alleles seem to be dominant for anal plate size but recessive for posterior lobe morphology.

As described above, we found a negative interaction between QTL contributing to clasper size on 3L and 3R (Table 1 and Figure S2). To determine if this effect is caused by dominance of D. mauritiana alleles or by epistasis between regions, we compared the clasper size of double heterozygotes, flies only heterozygous at one region, and D. simulans w501 males, our reasoning being that in the case of epistasis between loci, we should still detect a less-than-additive effect on the clasper size of double-heterozygous males. However, we found no evidence for less-than-additive interactions between regions C1 and C3 or between C1/C2 jointly and C3 (Table S7), which suggests that the interaction between the clasper QTL on the 3rd chromosome is caused by partial dominance of D. mauritiana alleles.

For clasper bristle number, we detected a significant negative interaction (i.e., smaller clasper than expected based on additive effects) between C1 alone and C3, as well as between C1/C2 together and C3 (Table S7), which suggests that there is epistasis at least between region C1 and region C3 for this trait. While we found no evidence for interactions between regions underlying posterior lobe size differences (Table S7), for posterior lobe shape we found a significant negative interaction between P1 or P2 and P4 (Table S1 and Table S7), as well as a positive interaction between P2 or P3 and P4 (Table S1 and Table S7). One of the caveats of this analysis, however, is that while a significant interaction would indicate epistasis between regions, a failure to detect it can result from the masking of the effect of D. mauritiana alleles by the dominant D. simulans alleles, and therefore, we cannot completely exclude epistasis between regions that show nonsignificant interactions (as in the case of posterior lobe shape).

Candidate genes for morphological diversification in genitalia

The small introgressed regions that we have identified still contain many genes that could potentially underlie the differences in terminal traits between D. mauritiana and D. simulans. Therefore we decided to test if positional candidate genes that exhibit differential expression in developing genitalia either between sexes (Chatterjee ) or between species (Masly ) are involved in the development of these traits. To do this, we conducted an RNAi screen in D. melanogaster using the NP6333-Gal4 driver line (Stieper ; Chatterjee ) to express double-stranded RNA (dsRNA) in the genital discs to knockdown such candidate genes located in the minimal interval (Figure 3A, blue bars and Table S8) of all chromosome 3L regions affecting the morphology of posterior lobes (P1, P2, and P3) and the smallest region affecting clasper morphology (C1).

Of 12 genes tested for region C1, we found that three of them, , , and affect the clasper development (Figure 4A and Table S8). Interestingly, and also appear to be involved in the development of the posterior lobe. While RNAi against had only subtle effects on posterior lobe size (Figure 4B and Table S8), RNAi against had a strong effect on both the shape and size of the posterior lobes (Figure 4, B and D and Table S8). In addition, we found significant effects of RNAi against two other C1 positional candidates, CG6673 and , on the shape and size of the posterior lobe, respectively (Figure 4, E and B and Table S8), although the claspers did not appear to be affected (Figure 4A and Table S8).

Figure 4

Effect of RNAi against positional candidate genes on clasper bristle number (A), posterior lobe size (B), and posterior lobe shape (C–H). Only lines with a significant difference with respect to the control (NP6333; UAS-Dicer) are shown. The results of the remaining crosses are available in Table S8.

In region P1, although 5 genes (of 17 tested) had an effect on posterior lobe size, only seemed to affect both the shape and size (Figure 4, B and F and Table S8). In region P2, we found 4 genes (of 18 tested) with an effect on posterior lobe size (, , , and ) (Figure 4B) and one gene () with an effect on the shape of this structure (Figure 4G and Table S8). Finally, in region P3, we found 2 genes (of 9 tested) with an effect on posterior lobe size ( and ) (Figure 4B) and 1 gene with an effect on posterior lobe shape () (Figure 4H and Table S8). RNAi against (P1) and (P2) had no effect on posterior lobe morphology but did cause clasper defects (Figure 4A).

In summary, of 56 positional candidate genes tested, we found that 20 of them seem to be required for the normal development either of the claspers and/or the posterior lobes. However, it is unlikely that all of these genes underlie natural variation for genital morphology. Instead, the effects detected are probably a result of the pleiotropic functions of a large proportion of the genes tested. Indeed, only 6 of these genes had effects consistent with the effect of the introgressed regions they are in and previously reported differences in expression in developing genitalia between species and sexes.

Discussion

The goal of our study was to investigate the genetic basis and architecture of interspecific differences in the morphology of anal and genital traits between D. simulans and D. mauritiana. The D. simulans species clade has a complex evolutionary history (Nunes ; Garrigan ) but the species split is thought to have occurred ∼250,000 years ago (McDermott and Kliman 2008; Garrigan ). While the size and shape of the posterior lobes has evolved in all species of the D. melanogaster group (Liu ), the morphology of the claspers and anal plates in D. mauritiana is probably derived, because in D. simulans and D. sechellia these structures are very similar in size and bristle number (True ).

The genetic basis of clasper size and clasper bristle number is coupled

We confirmed the location of a QTL for clasper bristle number on chromosome 3L found previously by True . However, in contrast to their results, the largest effect QTL for this trait in our mapping population was located on the X chromosome. We also found that at least three QTL, one on chromosome 2 and two on chromosome 3, are necessary to explain variation in clasper size between D. mauritiana and D. simulans. Since one of these QTL overlapped with the QTL for clasper bristle number identified previously (True ), we focused the fine-resolution mapping of clasper differences on chromosome 3 because this could potentially allow us to map both traits. We found that all regions containing at least one gene responsible for clasper size variation between the species also affect the number of bristles on this structure. This is unlikely to be due to close genetic linkage because at least one of these regions, C1, is very small (Table 2). Instead, it is possible that if genes act to regulate clasper size, nucleotide changes therein that make this structure larger may also result in the development of more bristles because their development is regulated at least in part through lateral inhibition (Heitzler and Simpson 1991).

Our QTL and introgression mapping data also indicate that the effect of the D. mauritiana X chromosome allele on clasper bristle number can be detected only if the D. mauritiana alleles are also present at our mapped region on 3L, which contributes to both clasper size and bristle number. This suggests that clasper bristle number is actually constrained by the size of the clasper.

The sum of effects of regions C1, C2, and C3 is strikingly larger for clasper bristle number (71%) than for clasper size (36%, Table 2). Assuming these estimates are accurate, this would mean that we may have uncovered most loci responsible for interspecific differences in clasper bristle number; while for clasper size, there are other QTL (including the one we detected on the 2nd chromosome). However, a more likely explanation is that the loci interact epistatically for clasper bristle number but not for clasper size (Mackay 2013). This is supported by our data: interactions between regions are additive for clasper size, but for clasper bristle number, we detected a negative epistatic interaction between regions C1 (and perhaps also C2) and C3. Furthermore, previous studies have found pervasive epistasis for other bristle traits (Long ; Gurganus ; Dilda and Mackay 2002). In addition, all clasper bristle number (and clasper size) effects are in the same direction (parental D. mauritiana), suggesting that, in contrast to anal plate bristle number (True ), clasper bristle number (and clasper size) may have evolved under directional selection in D. mauritiana (Zeng ).

The genetic basis of clasper divergence is independent of posterior lobe divergence

It is thought that, in general, the more pleiotropic a mutation is, the less likely it is to contribute to long-term adaptation (Stern and Orgogozo 2008; Wagner and Zhang 2011; Nunes ). Since it was found previously that the QTL for clasper bristle number on 3L overlapped with QTL for both posterior lobe and anal plate traits (True ; Zeng ), we tested whether the 3rd chromosome introgression lines affecting clasper morphology had an effect on these structures. However, we found that the regions responsible for clasper size and clasper bristle number variation have no detectable effect on the size of the anal plates and posterior lobes. This means that the association between clasper bristle number and posterior lobe and anal plate traits found previously is probably due to the low resolution of the QTL, rather than a common genetic basis, and therefore that independent genes underlie the diversification of these traits.

Regions underlying posterior lobe size divergence also affect the size of the anal plates

We found, however, that most regions underlying posterior lobe size differences also affected the size of the anal plates. This is a surprising result because the effect of the introgressions on anal plate size is not consistent with the direction of the difference between D. simulans and D. mauritiana. This means that in all cases, the D. mauritiana alleles cause a reduction in the size of both traits, indicating that these effects could be caused by the same genes. This also suggests that if the fixation of these alleles in D. mauritiana was driven by selection, it was likely due to their effect on the size of the posterior lobes rather than on the anal plates. Given the large estimated effects of the regions we have identified (Table 2), this also means that D. mauritiana may harbor additional major effect genes located elsewhere in the genome that compensate for the effects of our 3L region to give rise to the overall larger anal plates of this species. Despite exhibiting smaller anal plates, none of our introgression lines had a significantly different number of anal plate bristles from D. simulans w501, suggesting that the genetic regulation of anal plate development is different from that underlying clasper development.

Genetic architecture underlying genitalia evolution

“Haldane’s sieve” (Haldane 1924; Turner 1981) has been used to describe the lower fixation probabilities of recessive alleles and predicts that derived beneficial alleles are more likely to be dominant over ancestral ones (Turner 1981; Charlesworth 1992; Noor 1999). More recently however, it has been shown that when adaptation occurs from standing genetic variation or from alleles at mutation–selection balance, allelic dominance has little effect on the probability of fixation of the advantageous allele (Orr and Betancourt 2001; Hermisson and Pennings 2005). It appears that Haldane’s sieve only holds when adaptation occurs from new mutations, as the beneficial alleles will more quickly increase in frequency and more likely reach fixation if the mutation is dominant than if it is recessive (Orr and Betancourt 2001). We found that D. mauritiana alleles are dominant to D. simulans alleles in all regions affecting clasper morphology. Since the claspers of D. simulans and D. sechellia have very similar morphology (True ), large, “hairy” claspers are probably derived in D. mauritiana. Therefore, these results suggest that in this species, sexual selection on clasper traits may be acting not only through standing variation but also through new mutations (Barrett and Schluter 2008).

In contrast, in most regions that underlie posterior lobe and anal plate size differences, D. mauritiana alleles are recessive with respect to those of D. simulans. While we cannot infer the ancestral phenotype for posterior lobes (because the morphology of this structure has evolved in all three species of the D. simulans clade), it is likely that the anal plate size of D. mauritiana is derived (because D. simulans and D. sechellia have very similar anal plates), and therefore the causative genes in these regions must have evolved in D. mauritiana. Therefore, it is unlikely that the smaller posterior lobes evolved through new mutations in these regions in the D. mauritiana lineage (because their frequency in the population would increase very slowly, and they would have probably been lost due to genetic drift).

The mapping resolution that we have achieved is much higher relative to any previous studies of these traits (Liu ; Laurie ; True ; Macdonald and Goldstein 1999; Zeng ). However, even the smallest of these regions still contains 74 genes (Table 2). Therefore we performed RNAi experiments in D. melanogaster against positional candidates with expression differences either between sexes (Chatterjee ) or between species (Masly ) during genital development to identify promising candidate genes for future study in D. mauritiana and D. simulans. Furthermore, this strategy also allows us to exclude positional candidates that do not affect genital development when knocked down with the caveats that the RNAi has high penetrance and the positional candidate gene is also expressed in genital development in D. melanogaster.

Of 56 positional candidate genes with expression differences tested, we found that 20 are required for the normal developmental either of the claspers or the posterior lobes. However, only 6 of these genes had knockdown effects that were consistent both with the effects of the introgressed regions where they are located and expression differences either between sexes (Chatterjee ) or between species (Masly ) found in previous studies (Figure 4 and Table S8).

In region P1, there are two putative candidate genes, and , responsible for the difference in posterior lobe morphology between D. mauritiana and D. simulans. , a male sex lethal gene that affects the size and shape of posterior lobes (located in region P1), is expressed more highly in D. mauritiana genital discs relative to those of D. sechellia (a species that is morphologically more similar to D. simulans than to D. mauritiana) and that when knocked down in D. melanogaster results in elongated and overall larger posterior lobes (Figure 4, B and F and Table S8). RNAi against also affects posterior lobes (Figure 4B and Table S8) consistent with the difference in expression between species (Masly ).

Two genes in regions P2 ( and ) and P3 ( and ) were also consistent between the effect of the relevant introgressed regions (smaller posterior lobes) and RNAi effects as well as the expression differences between D. mauritiana and D. sechellia (i.e., smaller posterior lobes when the gene has lower expression in D. mauritiana than in D. sechellia, and larger posterior lobes when the gene has higher expression in D. mauritiana) (Figure 4B).

These six genes are compelling candidates for further studies to identify the evolved genes underlying diversification of genitalia morphology between D. mauritiana and D. simulans. However, we have not tested all of the genes in our mapped regions and it remains possible that genes for which no expression difference was detected previously (Chatterjee ; Masly ) are involved. Furthermore, it is also possible that at least some of the genes underlying genitalia differences between D. mauritiana and D. simulans have evolved through coding changes rather than through divergence in their expression.

Conclusions

Taken together our results demonstrate the complex genetic basis for fast evolving genitalia structures. Furthermore, we have identified novel roles in genital development for a number of genes, a few of which have effects consistent with the introgressed regions where they are located, and that are thus promising candidates for the divergence of these traits between D. simulans and D. mauritiana.