A shift to shorter cuticular hydrocarbons accompanies sexual isolation among Drosophila americana group populations

Abstract Because sensory signals often evolve rapidly, they could be instrumental in the emergence of reproductive isolation between species. However, pinpointing their specific contribution to isolating barriers, and the mechanisms underlying their divergence, remains challenging. Here, we demonstrate sexual isolation due to divergence in chemical signals between two populations of Drosophila americana (SC and NE) and one population of D. novamexicana, and dissect its underlying phenotypic and genetic mechanisms. Mating trials revealed strong sexual isolation between Drosophila novamexicana males and SC Drosophila americana females, as well as more moderate bi‐directional isolation between D. americana populations. Mating behavior data indicate SC D. americana males have the highest courtship efficiency and, unlike males of the other populations, are accepted by females of all species. Quantification of cuticular hydrocarbon (CHC) profiles—chemosensory signals that are used for species recognition and mate finding in Drosophila—shows that the SC D. americana population differs from the other populations primarily on the basis of compound carbon chain‐length. Moreover, manipulation of male CHC composition via heterospecific perfuming—specifically perfuming D. novamexicana males with SC D. americana males—abolishes their sexual isolation from these D. americana females. Of a set of candidates, a single gene—elongase CG17821—had patterns of gene expression consistent with a role in CHC differences between species. Sequence comparisons indicate D. novamexicana and our Nebraska (NE) D. americana population share a derived CG17821 truncation mutation that could also contribute to their shared “short” CHC phenotype. Together, these data suggest an evolutionary model for the origin and spread of this allele and its consequences for CHC divergence and sexual isolation in this group.

identify potential mates. Of these, many insect chemosensory signals primarily consist of cuticular hydrocarbons (CHCs)-a broad group of carbon-chain compounds important for many essential functions including sex pheromone signaling, but also environmental adaptation (Blomquist andBagneres 2010, Chung andCarrol 2015;Yew and Chung 2017). Unlike many auditory and visual cues, CHCs are produced and received by both males and females, making them potentially important for both male and female mate choice or species recognition. Indeed, studies on Drosophila chemical communication, particularly in the melanogaster subgroup-have detected clear evidence for sexual isolation based on male choice of female CHCs Oyama 1995, Billeter et al. 2009), while work in other species has demonstrated female choice of species-specific male CHC profiles Mas and Jallon 2005;Curtis et al. 2013;Dyer et al. 2014). Elements of the biochemical pathways of CHC production and the underlying genes are also relatively well understood in D. melanogaster (Pardy et al 2018). Accordingly, single genes contributing to species-specific CHC differences have been implicated as causes of reproductive isolation in several cases, including between D. melanogaster and D. simulans (desatF: Legendre et al. 2008) and D. serrata and D. birchii (mFAS: Chung et al. 2014). This framework provides an excellent resource for identifying candidate genes for sensory sexual signal variation in other Drosophila systems.
Here we use the Drosophila americana group to investigate how species variation in CHC profiles contributes to reproductive isolation, and to identify the underlying chemical and genetic changes that may comprise this variation. This group includes two closely related (MRCA ∼0.5 mya, Morales-Hojas et al. 2011) species-D. novamexicana and D. americana-that occupy distinct geographic and environmental habitats (Davis and Moyle 2019). D. americana is broadly distributed in the United States from the east coast to the Rocky Mountains, and exhibits significant phenotypic and genetic variation among populations throughout (e.g., Caletka and McAllister 2004;Davis and Moyle 2019), while D. novamexicana is localized to the arid southwestern US. Both species have been noted as being associated withand exclusively collected near-willow species of the genus Salix, though the exact nature of this association is not known (Blight and Romano 1953;McAllister 2002). Classical mating studies (Spieth 1951) have documented population-specific variation in reproductive isolation between members of this group (see Table S1 and below). Prior analysis has also shown qualitative differences in CHC composition between males (Bartelt at al. 1986), however the contribution of hydrocarbons to premating isolation has not been assessed. With three populations from this group, one D. novamexicana population and two D. americana populations-one southern (South Carolina, hereafter SC) and one western (Nebraska, hereafter NE)-here we use a com-bination of mating and behavioral studies, chemical analysis, and manipulation, and gene expression and sequence variation analyses to assess the role of CHCs in reproductive isolation. Together, our data suggest that evolution of a male sexual signal-an overall shift in the relative abundance of longer versus shorter cuticular hydrocarbons, due to novel mutation in an elongase gene-has produced complete premating isolation between derived males and females from species that retain the ancestral trait and preference, as proposed in classical models (Kaneshiro 1976(Kaneshiro , 1980 of the evolution of asymmetric sexual isolation.

SC D. am ericana FEMALES DISCRIMINATE AGAINST HETEROSPECIFIC MALES IN MATING TRIALS
We found clear evidence of moderate to strong sexual isolation between SC D. americana and the other two populations (Table 1). Mating rate (average proportion of females mated, in 4×4 mating trials) ranged from very high in intrapopulation crosses and some interpopulation crosses, to <50% in interpopulation combinations of female SC D. americana with each of the other species. Notably, D. novamexicana males were never successful in mating with SC D. americana females, indicating strong sexual isolation in this cross direction between these populations. In contrast, mating rate in the reciprocal cross was 70%. The other pairing with mating rates at or below 50% involved both reciprocal directions between D. americana populations. Across all cross combinations in the 4×4 mating trials, we found that female identity (Kruskal-Wallis test: χ 2 (2) = 8.99, P = 0.012), but not male identity (Kruskal-Wallis test: χ 2 (2) = 0.80, P = 0.67), significantly affected success. Post-hoc comparisons confirmed no pairwise differences in male mating rate between any populations (Table S2). In contrast, the copulation rate of SC D. americana females was significantly lower on average than both NE D. americana (P-adj = 0.025) and D. novamexicana (P-adj = 0.042), whereas NE D. americana and D. novamexicana females did not differ (P-adj = 0.15). Together, these results indicate that SC D. americana females discriminate against heteropopulation males more strongly than do females of the other two populations, with the greatest discrimination against D. novamexicana males. This discrimination produces strongly asymmetric sexual isolation between SC D. americana and D. novamexicana, while premating isolation between the two D. americana populations is bi-directional and more moderate. Importantly, this strong asymmetric isolation reiterates results documented by Spieth (1951,  (previously referred to as D. a. texana) (0-12.5% mating frequency; Table S1), while D. novamexicana females were more receptive in the reciprocal cross (mating frequency of 46.7-86.5%; Table S1). The consistency of these results 70 years apart suggests isolation observed here is not a bi-product of long-term laboratory culture, but instead reflects differences present in the natural populations from which these lines were collected.

EFFICIENCY AND MATING SUCCESS
We video-recorded single-pair matings for each cross combination to evaluate population differences in courtship strategies and whether these differed in interpopulation pairings. Three male courtship behaviors were quantified from these 1×1 trialsdisplay rate, tapping rate, and licking rate. All three behaviors showed similar patterns of among-population variation, when each was evaluated for the effect of male population, female population, or their interaction, using two-way ANOVAs. Males significantly differed with respect to display-rate (F(2, 36) = 3.76, P = 0.03) and tap-rate (F(2, 36) = 4.76, P = 0.017), but only marginally for lick-rate (F(2, 36) = 2.67, P = 0.083). For all three behaviors, this difference is due to SC D. americana males exhibiting higher rates (Tukey HSD post hoc tests; Figure 1; Table  S3). In contrast, we detected no female identity effects on male behavioral rate (display-rate: F(2, 42) = 1.83, P = 0.18; tap-rate F(2, 36) = 0.58, P = 0. 56; lick-rate: F(2, 36) = 0.43, P = 0.65), or any female × male interactions (display-rate: F(2, 36) = 0.72, P = 0.58; tap-rate F(2, 36) = 0.45, P = 0.77; lick-rate: F(2, 36) = 0.65, P = 0.63). Because behavior rates are dependent on copula-tion latency, and not simply the total number of behavioral events, these rates approximate the efficiency with which each courtship behavior results in a mating. Taken together, our results indicate that SC D. americana males have greater courtship efficiency (in terms of display and tap rates), and that males do not significantly vary these behavioral rates depending on female population identity. The single pair mating assays also broadly reiterated the patterns of moderate to strong sexual isolation we observed for SC D. americana in 4×4 trials. As with our 4×4 mating assay, D. novamexicana males never mated with SC D. americana females when paired individually, while mating rate between the two D. americana populations was moderately reduced in both directions. A logistic regression assessing how mating rate varied based on male identity, female identity, and their interaction showed a significant interaction effect only (males: χ 2 (2) = 2.97, P = 0.28; females: χ 2 (2) = 2.52, P = 0.23; males * females: χ 2 (4) = 12.67, P = 0.013)-a difference from 4×4-mating trials where female population identity was the primary predictor of mating rate. This variation between the assays appears to be driven by female NE D. americana accepting fewer D. novamexicana males in 1×1 trials (20%) compared to 4×4 trials (100% , Table 1), suggesting that male density within mating trials might affect copulation success in this particular species combination. Finally, copulation latency differed marginally based on both male identity (ANOVA: F(2) = 3.14, P = 0.055) and the male by female interaction (F(4) = 2.48, P = 0.062), but did not differ based on female species identity (F(2) = 0.68, P = 0.5147). Post-hoc tests suggest that the marginal male population effect is likely due to lower copulation latency in SC D. americana males relative to D. novamexicana males (Table S4).

COMPOSITION
We found that both populations and sexes within populations showed different and distinctive CHC profiles, with the largest differences detected between SC D. americana and the other two lines ( Figure 2). Across all samples (n = 5 replicate samples for each identity) we detected 8 alkene compounds (2 of which are sex-specific) and 4 methyl-branched alkane compounds present in at least one population. A principal component analysis of profiles from unmanipulated flies (unmanipulated-PCA or "U-PCA") summarizing the primary axes of CHC profile composition found that 95.3% of compound variation across all samples was explained by the first three principal components (U-PC1 to 3). Of these, CHC composition varied significantly for both population and sex for U-PC1 (pop: F(2, 429.29), P < 0.0001; sex: F(1, 46.97), P < 0.0001), and U-PC2 (pop: F(2, 10.54), P = 0.00043; sex: F(1, 181.69), P < 0.0001), but had no sex or population effect for U-PC3 (pop: F(2, 0.31), P = 0.737; sex: F(1, 0.57), P = 0.458; Figure 2). Notably, SC D.
americana of both sexes differed from the other two populations along the U-PC1 axis (Tukey HSD; Nov-SC: P-adjust < 0.0001; NE-SC: P-adjust < 0.0001; Nov-NE: P-adjust = 0.28), which also explains most of the CHC variation between samples (69.9%). U-PC1 was positively loaded for most of the shorter carbon chain length compounds in both alkene and methylbranched alkane classes of compounds, and negative values were strongly loaded for the compounds with longer carbon chain length (Table S5); consequently, this axis can be interpreted as a composite of average compound length across the CHC profile as a whole ( Figure 2). Accordingly, both sexes of SC D. americana had a higher abundance of longer-chain alkenes and methyl-branched alkanes than did either NE D. americana or D. novamexicana.
In contrast to U-PC1, U-PC2 appeared to primarily differentiate sexes within populations, but also NE D. americana males from males of the other two populations. This axis was most heavily loaded for two compounds: C21:1 and Me-C28, followed by C27:1 and Me-C26, with much smaller loadings for all other compounds (Table S5). C21:1 is a male-specific compound that is not detected in females; Me-C28, Me-26, and C27:1 abundance was also consistently different between sexes, with males always having more of these CHCs than females within each population ( Figure 2). Therefore, negative U-PC2 values can be interpreted as primarily representing more male-specific profiles, while positive U-PC2 values correspond to more female-like profiles. With respect to the detected species difference in U-PC2, posthoc tests reveal that this was driven by NE D. americana (Tukey HSD; Nov-SC: P-adjust = 0.99; NE-SC: P-adjust = 0.0062; Nov-NE: P-adjust = 0.0088). Furthermore, we find that abundance of the individual male-specific compound C21:1 differs between NE D. americana and males of the other two populations, but not between D. novamexicana and SC D. americana (Tukey HSD; Nov-SC: P-adjust = 0.98; NE-SC: P-adjust = 0.016; Nov-NE: P-adjust = 0.022). (Results from tests of individual compound differences can be found in Table S6).

INTER-POPULATION MALES
We found that patterns of sexual isolation could be modified by specifically changing the CHC profiles of males of different populations. To do this, we used perfuming assays to manipulate the pheromone profile of males by co-housing them with either intra-or inter-population males, and then evaluated how this manipulation influenced mating rates among populations. For these analyses, we focused on two pairings in two separate, analogous experiments-D. novamexicana and SC D. americana ("Nov-SC pair"), and the D. americana population pair ("NE-SC pair")because SC D. americana shows the strongest sexual isolation from the two other populations, and the largest differences in CHC composition. Within each perfuming experiment, we generated four combinations of target and donor male identities: each population perfumed by same-population males (control) and each population perfumed by hetero-population males. Each class of perfumed males was evaluated for mating rate specifically with SC D. americana females, because these females were the most discriminating against heteropopulation males in our previous mating assays ( Table 2). The observed mating rate of males perfumed with same-population males (the control manipulation) was similar to that observed in unperfumed mating trials: D. novamexicana and NE D. americana males perfumed with their own population performed poorly with SC D. americana females, relative to SC D. americana perfumed with their own males, which had a 100% mating rate (Table 2). In strong contrast, hetero-population perfumed males showed altered mating rates relative to control perfumed males. For Nov-SC pairings, D. novamexicana males perfumed with SC D. americana males successfully mated with SC D. americana females 94% of the time, a significantly higher mating rate than D. novamexicana males perfumed with their own males (25%; Mann-Whitney U-test: U = 0, P = 0.025). (Note that copulation success between conspecific perfumed D. novamexicana males and SC D. americana remains low (25%) but differs from the complete mating isolation observed in unperfumed trials, possibly because of changes in behavior due to male-male co-housing during the perfuming manipulation.) The reciprocal treatment also showed a significant effect of perfume source: male SC D. americana perfumed with D. novamexicana males displayed copulation success of 19%, significantly lower success compared to conspecific-perfumed SC D. americana (100%, Mann-Whitney U-test: U = 0, P = 0.018). For NE-SC pairings, SC D. americana male mating rate was also significantly reduced when perfumed with NE D. americana males (50%), compared to same-population perfumed SC D. americana males (100%, Mann-Whitney U-test: U = 0, P = 0.018). In contrast, mating rate of NE D. americana males perfumed with SC D. americana males increased slightly (56%), but did not significantly differ from the control (same-population-perfumed) males (44%, Mann-Whitney U-test: U = 5.5, P = 0.54).
Overall, these results indicate that perfuming males with heteropopulation CHCs influences the frequency of successful copulation with SC D. americana females. SC D. americana males perfumed with profiles of either D. novamexicana or NE D. americana had significantly reduced mating rates with their own females. Conversely, perfuming D. novamexicana males with SC D. americana males significantly increased their mating rate with D. americana females, almost completely reversing the pattern of sexual isolation observed for unmanipulated males. This large shift in mating rate for heteropopulation-perfumed D. novamexicana indicates that differences in male CHC profiles play a critical role in sexual isolation in the Nov-SC pair. In contrast with the three other classes of heteropopulation-perfumed males, perfuming NE D. americana with SC D. americana males did not significantly change their mating rate. It possible that this manipulation might have been less effective at altering the CHC profile of NE D. americana males (see next section), and/or that mating isolation in the NE-SC pair might depend on more complex factors than a simple shift in CHC composition alone (see Discussion).

THE DONOR MALE PROFILE
Using CHCs extracted from an additional set of perfumed males (see Methods), we confirmed that our heteropopulation perfuming manipulation produced quantitative changes in composite male CHC profiles in both the Nov-SC and NE-SC pairs, by shifting male CHCs closer to the donor male profile. PCAs were performed separately for each pairing (N-PCA for the Nov-SC pair, and A-PCA for the NE-SC pair, hereafter), as were analyses of differences among classes of perfumed males. In both cases, of the first three PCs of CHC composition, PC1 (i.e., N-PC1 or A-PC1) varied by both donor and target male identity (Table S7), indicating that perfuming significantly shifted CHC profiles along the primary axis of variation in both perfuming experiments. Nonetheless, the specific donor-target manipulation that was most successful in this regard differed between the Nov-SC and NE-SC pairs. In the Nov-SC pair, D. novamexicana males perfumed with SC D. americana males significantly differed in CHC composition from control (D. novamexicana) perfumed samples (N-PC1) (t(5.9) = −2.8, P = 0.031); however, heteropopulation-perfumed SC D. americana did not differ from same-population SC D. americana-perfumed samples along the same axis (t(3.53) = 1.91, P = 0.14). In contrast, in the NE-SC pair, SC D. americana perfumed with NE D. americana showed a significant shift in A-PC1 (t(4.85) = 4.23, P = 0.0088), but there was no significant shift for SC-perfumed NE D. americana compared to same-population-perfumed controls (t(4.68) = 0.055, P = 0.96). The difference between these pairs might be attributed to high variation seen among samples in some groups (in particular same-population perfumed SC D. americana). Regardless, it is clear from these results that this perfuming assay can produce significant, detectable differences in CHC profiles after heteropopulation perfuming, most notably in D. novamexicana males.

CONTRIBUTES TO CHC VARIATION BETWEEN SPECIES
From among a set of 23 candidate genes whose orthologs have functions related to cuticular hydrocarbon variation, both transcriptome and sequence comparisons implicated one specific elongase locus as potentially causal in CHC composition differences between our species. Using whole-transcriptome RNA-seq data from the same three populations (previously generated in Davis and Moyle 2020), we found that three of our candidate genes in males and two candidates in females were in the upper 10th percentile of genes most differentially expressed between populations (i.e. they showed a stronger species effect than 90% of all 11,301 expressed genes in the transcriptome dataset) (Table 3). Of these candidate genes, only CG17821 showed elevated gene expression specifically in SC D. americana compared to the other two populations; this pattern was observed in both sexes but is more pronounced in females ( Figure 3).
Closer inspection of our CG17821 sequences revealed that alleles in NE D. americana and D. novamexicana share a thymine insertion mutation that causes a premature stop codon four amino acids upstream of the end of the gene, compared to the annotated gene model in outgroup D. virilis and the allele in our SC D. americana stock ( Figure 3); these four terminal amino acids are not present in RNA transcripts for either NE D. americana and D. novamexicana. Among our three lines, this suggests an association between the truncated protein product of CG17821 and the "short" CHC phenotype we observe. This association is further supported by additional CHC and sequence data from another D. americana population, originally collected in Anderson, IN (species stock center line 15010-0951.00, henceforth IN D. americana). Lamb et al. 2020 characterized CHC divergence between females from this IN D. americana stock and another D. novamexicana line (species stock center 15010-1031.04). The profile for the D. novamexicana line had no long compounds (>C30), a low abundance of C29 compounds, and presence of multiple short alkenes (<C27), consistent with the "short" CHC profile we observe for our D. novamexicana line. Likewise, the IN D. americana profile they observed is broadly consistent with the "long" profile we observed in SC D. americana, as both show presence of compounds longer than 30 carbons, a high abundance of C29 length alkenes, and absence of compounds shorter than 27 carbons. In addition, using a genome assembly from the same IN D. americana stock (Kim et al. 2020), we found that this population also has the non-truncated allele of CG17821. These data are consistent with our inference that truncation of this allele could contribute to the difference between "long" and "short" CHC phenotypes among populations in this group.

EVOLUTIONARY HISTORY OF ELONGASE CG17821
Several lines of evidence indicate that CG17821 has a complex evolutionary history within the D. americana subgroup that likely includes post-speciation introgression. The presence of the nontruncated allele in outgroup D. virilis ( Figure 3) indicates that the truncated allele is derived from a thymine insertion mutation event that took place within the D. americana subgroup. Moreover, D. novamexicana and our western D. americana (NE) line share an identical derived allele at CG17821, indicating an evolutionary history at this locus that disagrees with expected phylogenetic relationships among these three taxa (i.e., where the two D. americana populations are expected to be most closely related). Gene trees for each of our 23 candidate genes confirmed that CG17821 has a topology that is discordant with the expected species relationships, placing NE D. americana as sister to D. novamexicana rather than grouping it with SC D. americana (Figure 3). Four other CHC candidate genes also show this discordant topology, most notably CG18609 that is located immediately downstream (within 1kB) of CG17821, according to the annotated genome of D. virilis ( Figure 3). For four out of five of these loci-including both CG17821 and CG18609-using the allele from the IN D. americana stock instead produces a genealogy where the D. americana populations group together, as expected from the species tree ( Figure 3).
The observation of phylogenetically discordant sites can be due to several factors, notably incomplete lineage sorting (ILS) or post-speciation introgression (see also Discussion). Two lines of evidence strongly support introgression. First, evidence from SNP variants, across the whole transcriptome and specifically at CG17821, is more consistent with introgression between D. novamexicana and our western D. americana (NE) populationas determined by Patterson's D-statistic (Durand et al. 2011). In particular, we found that genome-wide, of 34,441 total SNPs detected in 11,301 loci within the transcriptome dataset (Davis and Moyle 2020), 14,902 supported D. americana populations as sister (in accordance with the species tree), 13,372 SNPs grouped D. novamexicana and NE D. americana, and 6167 grouped D. novamexicana with SC D. americana. The significant excess of shared variants between D. novamexicana and our NE D. americana line (D = 0.369, P < 0.0001), is consistent with a history of introgression between the progenitors of these two populations. Second, gene tree topologies at and around CG17821 also indicate this specific genomic region shares recent ancestry between D. novamexicana and NE D. americana due to introgression. CG17821, the adjacent downstream CHC candidate CG18609, and the next nearest gene (List, a neurotransmitter located ∼14 kB downstream) all display shared derived SNPs between D. novamexicana and the NE D. americana line, and discordant gene tree topologies that group these two lines as sister taxa. In contrast, topologies for the next two nearest up-or down-stream genes in this region reflect the species tree, indicating that recent shared ancestry between D. novamexicana and the NE D. americana extends across a genomic window between 18 and 68 kB long around the specific region containing CG17821 ( Figure 3C). Estimates of linkage disequilibrium (LD) from D. melanogaster autosomes indicate that most LD decays within Bold genes in these columns indicate genes in the top 10% of differentially expressed genes within sex. 200 base pairs, and r 2 (the correlation between SNPs) decreases to <0.1 within a 1 kB window (Franssen et al. 2015); therefore, the size of the putatively-introgressed window observed here is well outside the range expected under ILS, as discordant ancestry due to sorting from ancestral variation is expected in small blocks (Hudson and Coyne 2002;Gao et al. 2015). Note that inversions could also be responsible for capturing ancestrally segregating variation in larger genomic regions than expected from LD measures, and populations within this group are known to differ in the presence/absence of inversions, including at the specific genomic region containing these genes (Reis et al. 2018). However, the inversion at this location is only present in southern D. americana (i.e., within the "D. a. texana" chromosomal form) to the exclusion of NE D. americana, D. novamexicana, and D. virilis, and therefore could not be responsible for the patterns of discordant variation we observe for these genes here. Coupled with genetic and phenotypic evidence above, these results indicate introgression of these genes between D. novamexicana and western D. americana-and not ILS-is most likely responsible for their shared allele at CG17821 and therefore potentially for their similarity in "short" CHC phenotypes. We also evaluated whether there was evidence of elevated protein evolution at CG17821, CG18609, or any other of our candidate genes, based on estimates of the ratio of nonsynonymous to synonymous substitutions (d N /d S ) at these loci, compared to a genome-wide average rate estimated from our transcriptomewide dataset (see Methods). For 4397 genes transcriptome-wide, the median and mean d N /d S were 0.065 and 0.0991, respectively (± 0.133 s.d.) ( Figure S3). Of our candidate genes, only three had an estimated d N /d S greater than 1 standard deviation above this mean, including the two candidates within the putatively introgressed region-CG17821 (d N /d S = 0.284) and CG18609 (d N /d S = 0.249)-as well as CG6660 (d N /d S = 0.38714)(full results in Table S8). Notably, these rates of protein change fall within the top 8% of all genes analyzed ( Figure S3, all gene data reported in Supporting Information). While these estimates are below the threshold for unambiguous evidence of positive selection (d N /d S > 1), only nine genes in the whole dataset met this criterion. Of the three candidate loci that exhibit elevated protein evolution, CG6660 does not show patterns of genealogical, allelic, or gene expression variation that match our observed CHC phenotypic variation ( Figure 3A and B, and above). For both CG17821 and CG18609, we also estimated d N /d S using just D. virilis and SC D. americana sequence comparisons, to determine if elevated protein evolution in these genes is solely driven by the branch leading to the D. novamexicana-NE D. americana allele, or if this pattern is more general across the clade. We found our estimate of d N /d S was still modestly elevated between D. virilis and SC D. americana (d N /d S = 0.284) even when this branch was removedindicating CG17821 has experienced consistently elevated pro-tein evolution across the whole clade. In contrast, CG18609 protein evolution is estimated to be lower (d N /d S = 0.147) when just considering divergence between D. virilis and SC D. americana, suggesting that most of the accelerated protein evolution in this gene occurred after the split of the D. novamexicana/NE D. americana lineage. The timing of elevated protein evolution in CG18609 therefore also generally coincides with CHC profile differences observed here and, like CG17821, CG18609 is also an elongase. However, CG18609 does not show patterns of gene expression that match the "long" versus "short" CHC phenotypes we observe (instead, its expression is significantly reduced in NE D. americana; Figure 3A) so while it's possible that this tandem elongase gene also plays a role in CHC phenotype divergence, it is unlikely to be responsible for the primary axis of CHC variation described here.
Finally, we note that our d N /d S analysis also revealed two odorant-binding proteins (OBPs)-to be among the fastest evolving loci in our transcriptome-wide dataset: the orthologs of Obp99d (second fastest in our dataset; d N /d S = 1.589), and Obp56g (40th fastest; d N /d S = 0.689; see Supporting Information). OBPs are thought to facilitate olfactory processing by chaperoning odorants to the olfactory receptor neurons or terminating neuron activity by clearing odorants form the surrounding (Sun et al. 2018). OBP genes have been previously shown to evolve quite rapidly in multiple insect groups (Foret and Maleszka 2006), although rates in these specific OBPs are lower across six species in the Drosophila melanogaster group (d N /d S > 0.24; Vieira et al. 2007) than we detect here. Notably, both these OBPs are known to expressed in chemosensory organs in D. melanogaster-including adult labellum (Obp56g; Galindo and Smith 2001), antenna and maxillary palps (Obp99d; Hekmat-Scafe et al. 2002), and wing sensilla (both loci; He at al. 2019)consistent with roles in pheromone perception. Moreover, one of these loci-Obp56g-is a known seminal fluid protein in the D. melanogaster group (Findlay et al 2008), that also changes in expression within females specifically in response to mating (Mc-Graw et al. 2004). Given these roles, the elevated rates we observe here suggest both genes might contribute to variation in behavioral responses to pheromone and other sexual stimuli among our species (see also discussion).

Discussion
Identifying genetic and evolutionary mechanisms involved in the earliest steps of reproductive isolation between species is essential for understanding the factors that drive speciation. Evolutionary divergence in sexual signals may be an especially potent contributor to this process (Schemske 2000;Coyne and Orr 2004;Ritchie 2007, Schluter 2009, van Doorn et al. 2009). However, demonstrating the connection between sensory signal divergence and emerging reproductive isolation can be challenging, as it requires identification and demonstration of the direct role of specific signals mediating sexual isolation between species and knowledge of the specific mechanistic changes that have given rise to lineage differences in this signal. Here we have demonstrated that sexual isolation between laboratory populations in the D. americana group is based on female choice of male chemical signals, and identified both the specific phenotypic shift between species in pheromone chemistry as well as a genetic variant likely contributing to this phenotypic change and the mating isolation that results from it. Together these data support a clear role for sensory signal divergence in the evolution of premating isolating barriers between populations in this closely related group, and provide insight into how relatively simple genetic and phenotypic mechanisms can cause strong isolation even at early stages of evolutionary divergence.

RESPONSIBLE FOR ISOLATION BETWEEN D. novam exicana AND SC D. am ericana
Differences in copulation success between populations are the product of variation in male choice (via differences in courtship intensity and copulation attempts, depending upon female identity), female choice (via differences in acceptance rates of males, depending upon their identity), or a combination of these factors. Disentangling these alternatives is critical for identifying the specific trait and corresponding preference variation responsible for the emergence of prezygotic species barriers. Here we showed that female choice of a male sensory signal results in strong sexual isolation among D. americana group populations, specifically SC D. americana female preference for CHC profiles of their own males and rejection of males with alternative profiles. The role of CHCs is exceptionally clear in the case of premating isolation between male D. novamexicana and female SC D. americana: strong sexual isolation in this cross can be almost entirely reversed by perfuming D. novamexicana males with SC D. americana male CHCs. In contrast, we find little evidence for differential male preference of females, even though female CHC profiles are also divergent between populations. Instead, our behavioral data indicate that males did not alter their courtship behavior in response to female species identity. Interestingly, these observations also suggest evidence for female choice of male courtship behaviors, whereby SC D. americana male courtship consistently induced females to copulate sooner than courtship behaviors of other population males. Moreover, our finding of strong asymmetric mating isolation between male D. novamexicana and specific (southeastern) populations of D. americana recapitulates patterns of isolation described in this group more than 70 years ago (Spieth 1951; Table S1), indicating that these mat-ing isolation patterns reflect natural variation among populations from which these lines were collected.
These findings fit within a body of studies that have identified either male or female choice of sensory signals as critical for sexual isolation among Drosophila species (patterns reviewed Yukilevich and Patterson 2019). In Drosophila melanogaster, many studies have found that male choice of specific female CHC compounds play a role in isolation between closely related heterospecifics (Coyne et al. 1995, Billeter et al. 2009) as well as between intraspecific populations (Wu et al. 1995;Hollocher et al. 1997;Yukilevich and True 2008). Such patterns of CHC-mediated male mate-discrimination have also been associated with allelic variation in CHC elongases. For example, D. sechellia reproductive isolation from closely related D. simulans results from male mate-discrimination based on female CHC profiles (Shahandeh et al. 2018) and the CHC elongase eloF has been demonstrated to inhibit interspecific mating in the same species (Combs et al. 2018). Female choice of male CHC variation has received comparably less mechanistic attention, however has been demonstrated as a prezygotic barrier in several Drosophila groups, most notably between D. santomea and D. yakuba Mas and Jallon 2005), between the mycophagous D. subquinaria and D. recens (Curtis et al. 2013, Dyer et al. 2014, and between D. mojavensis males reared on different cactus substrates (Havens and Etges 2013). Our results, therefore, complement a growing body of evidence that shows female choice can be an important factor dictating patterns of sexual isolation among closely related Drosophila populations and species.

CHC DIVERGENCE AND THE ROLE OF ELONGASES
Dissecting the finer details of phenotypic divergence in sensory signals can help pinpoint underlying mechanisms and associated genes responsible, and more clearly demonstrate how these signal changes contribute to isolation between populations. Here we found that CHC divergence between our lines occurs primarily on the basis of compound length, with D. novamexicana and NE D. americana profiles both having similar enrichment of shorter compounds ("short" phenotype) compared to an enrichment of longer compounds for SC D. americana ("long" phenotype). These differences in abundance of shorter-versus longer-chain compounds were observed across both sexes, and across both alkenes and methyl-branched alkanes (Figure 2, Table S6). In contrast, we find no evidence for variation in other features such as double bond or methyl branch location or number. These features themselves suggest that the striking difference between "shorter" and "longer" profile phenotypes could be due to variation in fattyacid elongase activity, which globally influences the carbon chain length of CHC precursors, and therefore can have consistent downstream effects on both alkenes and methyl-branched alkanes, because both are modified after the elongation step (Pardy et al. 2018).
Strikingly, analyses of gene expression and sequence variation revealed that CG17821 and CG18609-two putative fattyacid elongases (Szafer-Glusman et al. 2008;Gaudet et al. 2011)-could functionally contribute to our observed variation between longer and shorter CHC profiles. Sequences at both loci are identical between NE D. americana and D. novamexicana, differ from SC and IN D. americana lines, and show modestly elevated protein evolution on the branch leading to D. novamexicana/NE D. americana, all of which broadly coincide with our observed differences between "short" versus "long" CHC phenotypes. The variation we observed at CG17821 is particularly interesting, because we observe a thymine insertion mutation that results in premature truncation of CG17821 specifically in D. novamexicana and NE D. americana ( Figure 3) and because our gene expression data also indicate reduced expression of the CG17821 allele specifically in these lines. In the latter case, while we might expect gene expression differences to be tissue-specific (that is, observed specifically in the oenocytes: the CHC-producing organs), our data indicate that the populationspecific expression signal at this locus is strong enough to detect from whole-body transcriptome data. Both observed sequence and expression changes at CG17821 could individually produce a greater abundance of short CHC products either by lowering elongase protein levels or reducing enzyme activity. Therefore, although our current data cannot differentiate which change may have occurred first (or whether they are pleiotropic), either could produce the pattern of phenotypic difference observed between populations. Overall, these data suggest allelic variation in CG17821 primarily underlies the major axis of CHC divergence observed between SC D. americana and the other populations, consistent with a hypothesis of a simple underlying basis for chain length variation. Moreover, because this axis of CHC divergence appears primarily responsible for sexual isolation between D. novamexicana and SC D. americana, this points to a large role for simple allelic change at this genomic location in the emergence of a strong isolating barrier between these two populations.

SEXUAL ISOLATION IN THIS GROUP
Together, our data demonstrate sexual isolation between SC D. americana and D. novamexicana is due to CHC divergence in compound length and suggest that variation between truncated and non-truncated alleles of the elongase CG17821 contribute to this phenotypic variation. Moreover, both genome-wide SNP variation, as well as localized variation specifically around this locus ( Figure 3C), indicate that this phenotypic variation involves a history of introgression. These data point to a model hypothesis for the evolutionary history of transitions involved in the change in CHC profiles between species and, potentially, in the emergence and expression of sexual isolation that depends upon this phenotype.
First, the distribution of both "long" versus "short" CHC phenotypes (shown here and in Lamb et al. 2020), and allelic variation CG17821, indicate that the "short" phenotype and the CG17821 truncated allele are derived states that arose within the D. americana group. This shift most likely occurred in western lineages that gave rise to contemporary D. novamexicana. The evolutionary forces responsible for the persistence and spread of this phenotype are not yet known. The relationship between variation in environmental factors, insect stress physiology, and features of CHC length and branching is known to be complex. Prior work has associated aspects of CHC divergence with abiotic variation such as latitude (Frentiu and Chenoweth 2010;Rajpurohit et al. 2017, but see Gibbs et al. 2003), and physiological traits such as desiccation resistance (reviewed Chung and Carroll 2015), that suggest a role for natural selection in shaping CHC composition, but the specific sources(s) of selection can be challenging to pinpoint. In the D. americana group, species habitats are differentiated primarily on the basis of water availability and these two species, as well as populations within them, differ in key physiological traits such as desiccation resistance (Davis and Moyle 2019). However, sexual selection might also contribute to shaping evolution at CHC loci-as evidenced by the importance of CHC variation for mating success outlined here and elsewhere. We also find evidence that two odorant binding proteins (OBPs)-Obp99d and Obp56g-are rapidly evolving across this group. OBPs are known to be important for mediating olfactory behavioral responses during sexual interactions (Laughlin et al. 2008;Leal 2013, Sun et al. 2018) as well as host-plant preference (Matsuo et al. 2007, Comeault et al. 2017. Therefore these loci could be evolving due to natural selection, sexual selection, or both, including in response to changes in pheromone profiles described here, or to other factors such as this group's close but little investigated habitat association with willow (Salix sp.) trees (Blight and Romano 1953;McAllister 2002, personal observations). Interestingly, our analysis indicates that the elongase CG17821 has experienced modestly elevated protein evolution across the whole D. virilis sub-clade; this suggested history of sustained selection indicates this locus might have played important roles in long-term CHC-mediated adaptive divergence across this group. In comparison, the downstream elongase CG18609 has evidence of elevated protein change primarily on the branch leading to D. novamexicana/NE D. americana, possibly suggesting that this acceleration occurred after impactful changes to CG17821.
Based on our data associating phenotypic divergence with sexual isolation, the appearance of this new CHC phenotype would have reduced sexual compatibility between derived "short" males and females with strong preferences for the "long" ancestral CHC profile. Persistence of this phenotype would have required a broadening of female preference to accommodate males with the derived ("short" CHC) pheromone phenotype (e.g., as has been observed, for example, in male Ostrinia moths; Roelofs et al. 2002). Our data support this expectation, as D. novamexicana females are more accepting of both (putatively derived) D. novamexicana and (ancestral) SC D. americana male CHC phenotypes, while SC D. americana females discriminate against derived "short" CHC phenotypes. Interestingly, this model for the evolution of asymmetric sexual isolation is broadly consistent with Kaneshiro's (1976Kaneshiro's ( , 1980 model for peripatric speciation. Kaneshiro observed that females from derived populations frequently have broad preferences for both derived and ancestral male phenotypes; he proposed that this was due to relaxed selection on narrow female preferences in genetically bottlenecked island populations, where founder effects have led to the loss of elements of male courtship. Although Kaneshiro's model explicitly invokes genetic drift in male trait evolution, our observations indicate his model for the origin of sexual isolation asymmetry could extend more generally to any case where evolutionary change affects a trait important for male sexual signaling. In the case described here, evolutionary change in male CHC profiles (possibly due to selection acting on a CHC elongase gene(s)), accompanied by an apparent broadening of female preferences for these profiles in derived populations, has resulted in the emergence of strong premating asymmetry specifically between females with ancestral trait preferences and males with derived trait values-akin to the model outlined by Kaneshiro. Intriguingly, our data also support the subsequent movement of this CHC phenotype from D. novamexicana into western D. americana lineages. One possible explanation for shared variation in the CG17821-CG18609 locus and CHC phenotypes is that this arose from segregating ancestral variation present in both D. americana and D. novamexicana (that is, is due to ILS). Prior evidence (Caletka and McAllister 2004) as well as data here indicate that some observed site discordance between populations of D. americana and D. novamexicana is consistent with ILS. However, our analysis strongly supports the additional occurrence of introgression between D. novamexicana and western D. americana, including specifically of this trait from D. novamexicana into western D. americana populations. Both their significant excess of shared genome-wide variation (as indicated by the Dstatistic), and a genomic region of at least 18 kB of shared recent ancestry surrounding CG17821/CG18609 that accompanies a similar shift to the "short" CHC phenotype in western D. americana, support this inference. Note also that while D. novamexi-cana and D. americana are reported to be allopatric, limited and sporadic field collections of D. novamexicana mean there is an incomplete understanding of the density and extent of this species' historical range. As a result, this species could have been in closer geographical contact with western D. americana populations during the period since their initial split (∼500,000 years ago), which could help explain evidence for introgression after speciation inferred here. Recent work (Sramkowski et al. 2020) describing rare D. novamexicana-like pigmentation alleles in geographically disparate D. americana populations, similarly suggests evidence of more recent gene exchange between these species.
Given the effect of the "short" CHC phenotype on sexual isolation between D. novamexicana and the SC D. americana population, this introgression is expected to have consequences for reproductive isolation among D. americana lineages. Interestingly, our data indicate that, even though NE D. americana shows the general shift to shorter chain CHC length associated with the introgressed region, its current patterns of sexual isolation differ from those seen in D. novamexicana. One explanation might be that this introgression-mediated shift in CHCs occurred on a novel D. americana genomic background, with potential consequences for the expression of introgressed CHCaffecting loci, and therefore for patterns and strengths of CHCmediated isolation. A concrete example of these background effects can be seen for CG18609, where NE D. americana shares the D. novamexicana allele but nonetheless exhibits significantly reduced expression of this locus compared to the other two lineages ( Figure 3A). A differential history of allelic exchange (with D. novamexicana) across the range of D. americana, plus variable genomic background effects on the expression of CHC loci, could contribute to the more complex mating relationships observed here, and elsewhere, among D. americana populations. For example, Spieth's mating analyses (Table S1) identified up to 10-fold variation in mating success among nine disparate populations of D. americana. The possibility that this variation is influenced by differential gene exchange with D. novamexicana (including of loci with major effects on a signaling phenotype) is testable with work examining mating success between geographically diverse D. americana populations-particularly between eastern and western populations-and its covariation with CHC phenotypic and genotypic variation.
Regardless of the collateral consequences for isolation among D. americana populations, our data clearly support the role of divergent cuticular hydrocarbon profiles-specifically a general shift in carbon chain length-in sexual isolation between our D. novamexicana and SC D. americana populations. They also implicate a potentially causal role for the gene CG17821 in determining CHC phenotype via a global change in CHC elongation activity early in the generation of these compounds. These data provide a strong example of how a recently derived allele in a single gene with large phenotypic effects on a sexual signal could underpin asymmetric sexual isolation between closely related species. Moreover, they suggest multiple (behavioral, biochemical, and molecular) lines of evidence that chemosensory processes are evolving rapidly and dynamically across this group.

EXPERIMENTAL FLY STOCKS
Three stocks were obtained from the University of California San Diego Drosophila Species Stock Center (DSSC): a Drosophila novamexicana stock from San Antonio, NM (15010-1031.08); and two D. americana stocks, one from Chadron, NE (15010-0951.06, NE D. americana throughout); and one from Jamestown, SC, (15010-1041.29, SC D. americana throughout). All stocks were originally collected between 1946 and 1953. D. americana has sometimes been divided into two subspecies according to presence (D. americana americana) or absence (D. a. texana) of a chromosomal fusion of the X-and 4-chromosomes that shows a distinct latitudinal cline (McAllister 2002), however because sub-specific differences apart from this fusion have not been consistently supported, our two lines are treated as populations from within a single heterogeneous species here. The D. novamexicana and NE D. americana populations used here are the same as those collected and used by Spieth 1951. All fly stocks were reared on standard cornmeal media prepared by the Bloomington Drosophila Stock Center (BDSC) at Indiana University, and were kept at room temperature (∼22°C). Every assay in this study used virgins isolated within 8 hours of eclosion and aged for 7 days prior to the start of experiments, similar to the 8 days used by Spieth (1951).

4×4 UNPERFUMED MATING ASSAY
We performed trials in which four virgin males and four virgin females were paired and observed for mating behavior, following the design used by Spieth (1951) that allows for behavioral interactions that might not otherwise be observed in similar singlepair assays. Within each trial, all males are from a single population, as are all females, so are no-choice with respect to the genotype of a mating partner; crosses are varied by pairing males and females of alternative lines. For each trial, four males and four females were transferred to a single vial without anesthetization and observed for 3 h. The number and duration of each copulation event were recorded for each trial. This assay was repeated for a total of five replicates for each possible population combination, in reciprocal (i.e., 9 cross types, each of N = 5). Each trial used 7-day old virgins and testing was started within 30 minutes of lights on in the morning.

COURTSHIP BEHAVIOR ASSAY
To quantify and evaluate differences in courtship behaviors between cross types, flies were observed in single pair (no-choice) mating assays. We used a modified FlyPi setup-that combines a Raspberry Pi (Raspberry Pi Foundation, Cambridge, UK), pi camera, and 3D printed parts (Chagas et al. 2017)-to record courtship behaviors. Assays were performed in a modified cell culture plate consisting of six 3 cm-diameter culture wells, each with a small amount of cornmeal media in the bottom, that allowed six total crosses to be recorded simultaneously. For each assay, individual virgin male and female flies of a given cross were aspirated without anesthetization to a cell culture well; after six total crosses were set up, the plate was videotaped for a contiguous 3-h period. The six cross combinations assessed in any particular video trial were randomized to account for variance that might otherwise be explained by date. As in the 4×4 mating assay, we performed five replicates of all possible population combinations in reciprocal.
Behavioral features were analyzed and scored manually by the same individual to avoid subjective variation among researchers. Three courtship behaviors-male display events, male tapping events, and male licking events-in addition to copulation were scored in each 1×1 trial using the following criteria. A male "display" was counted when a male performed a combination of back and forth movements and occasional wing flicks while maintaining sustained orientation in front of and facing the female. "Tapping" events were defined when a male used his tarsus to touch the abdomen of the female when oriented behind her. "Licking" events were defined when the male's mouthparts contacted the female genital arch. For male display, tapping, and licking events, individual events were scored separately only if a different behavior (including sitting still/walking away) was observed between instances of the defined behavior. This criterion minimized overcounting of discrete behavioral events, especially those with difficult to view aspects such as number of times a male extruded mouthparts during a contiguous licking event. In addition to these behaviors, copulation success and latency to copulation were also recorded for each trial. Note that females also engage in tapping and other behaviors, but because these appear to be less consistent and are more difficult to observe and score, they were not addressed in this study.

HYDROCARBONS
Cuticular hydrocarbons were extracted from pooled samples by placing five 7-day old virgin flies of a single sex and species identity in a 1.8 mL glass vial (Wheaton 224740 E-C Clear Glass Sample Vials) with 120 μL of hexane (Sigma Aldrich, St Louis, MO, USA) spiked with 10 μg/mL of hexacosane (Sigma Aldrich). After 20 min, 100 μL of the solution was removed to a sterilized 1.8 mL glass vial (Wheaton 224740 E-C Clear Glass Sample Vials) and allowed to evaporate overnight under a fume hood. Extracts were stored at −20°C until analysis. Five replicate samples consisting of five flies per sample (25 flies total) were prepared for each sample type, and all replicates were extracted on the same day.
Gas chromatography-mass spectrometry (GC/MS) analysis was performed on a 7820A GC system equipped with a 5975 Mass Selective Detector (Agilent Technologies, Inc., Santa Clara, CA, USA) and a HP-5ms column ((5%-Phenyl)methylpolysiloxane, 30 m length, 250 μm ID, 0.25 μm film thickness; Agilent Technologies, Inc.). Electron ionization (EI) energy was set at 70 eV. One microliter of the sample was injected in splitless mode and analyzed with helium flow at 1 mL/ min. Two different temperature gradients were used depending on the sample type. For CHC analysis of unmanipulated males and females of each species, the following parameters were used: column was set at 50°C for 0 min, increased to 210°C at a rate of 35°C/min, then increased to 280°C at a rate of 3°C/min. The MS was set to detect from m/z 33 to 500. For analysis of samples from perfuming trials, the parameters were modified to increase resolution and sensitivity for less abundant compounds: the column was set at 40°C and held for 3 min, increased to 200°C at a rate of 35°C/min, then increased to 280°C at a rate of 3°C/min and held for 15 min. Chromatograms and spectra were analyzed using MSD ChemStation (Agilent Technologies, Inc.). CHCs were identified on the basis of retention time and electron ionization fragmentation pattern. Compounds are identified in this study as: CXX:Y for alkenes or Me-CXX for methyl-branched alkanes, where XX indicates the length of the carbon chain, and Y indicates number of double bonds, e.g., C21:1 is a 21-carbon alkene with a single double bond.
The abundance of each compound was quantified by normalizing the area under each CHC peak to the area of the hexacosane signal using homebuilt peak selection software (personal correspondence, Dr. Scott Pletcher, Univ. of Michigan).

PERFUMING MANIPULATION AND MATING ASSAY
"Perfuming" involves co-housing target flies in vials filled predominantly with flies of a desired donor identity, so that the CHC profile of the target flies is altered via physical transfer of CHCs from donor flies (Coyne et al. 1994;Dyer et al. 2014;Serrato-Capuchina et al. 2020). Perfuming was performed by placing two "target" males with 15 "donor" males within a single vial. All males were 1-day old virgins when perfuming vials were established, and the wings of donor males were removed (under anesthetic) to distinguish them from target males. For the hetero-population perfuming treatments, target males were co-housed with donor males of a different population; samepopulation (control) treatments paired donor and target males of the same identity. Following 7 days of perfuming, four male flies from the same perfuming conditions were transferred (without anesthetization) to vials containing 4 female SC D. americana. Within each experiment (Nov-SC, or NE-SC), we ran one trial of each perfuming condition (2 heterospecific-perfumed types, and 2 conspecific-perfumed types) in parallel on the same day (four 4×4 trials in total).
For CHC analysis, males were perfumed as described for the perfumed mating assay, for both Nov-SC and NE-SC experimental pairs, with the exception that two to three target males (rather than just 2) were co-housed with 15-18 donor males, enabling two parallel perfuming vials to generate five target male flies per identity, for each individual CHC extraction. This was then replicated four times for each identity to reach N = 4 biological replicates for this analysis.

CANDIDATE GENE SELECTION
Our candidate list was generated by searching Flybase (Flybase.org) for annotated Drosophila melanogaster genes with one of the protein-coding domains (as identified by InterProt) that have known functions in CHC synthesis (Pardy et al. 2018): "fatty acid desaturase," "fatty acid desaturase domain," "Cyp4g," "ELO family," or "fatty-acid-synthase"-resulting in an initial list of 34 genes. To this, we added the pigmentation genes ebony and tan as they have been shown to alter CHC variation among both D. melanogaster (Massey et al. 2019) and D. americana group species (Lamb et al. 2020). With these 36 genes, we identified orthologs in Drosophila virilis (via Flybase using OrthoDB version 9.1, Zdobnov et al. 2017) and then evaluated whether each of these loci had transcript expression in our three populations. To do so, we used the BLASTn function of BLAST+ version 2.6.0 (Camacho et al. 2009) to search for matches between the D. virilis orthologs and previously published wholebody transcriptome data generated from the same three populations using RNA-seq (Davis and Moyle 2020). Our analyses used only gene expression data from the control (ambient) conditions for both males and females from this study, and excluded any desiccation stress treatment data. Of the 36 initial genes (listed in Table S9), 30 were found to have unambiguous 1-to-1 orthologs in D. virilis and, of these, 23 had transcripts present within the Davis and Moyle (2020) dataset. This final set of 23 genes was used to evaluate gene expression and sequence variation among our three lines. For downstream analyses, we also identified alleles of these 23 loci in a genome assembly of D. americana Anderson (15010-0951.00, also referred to as A01) generated by Kim et al. 2020, using the BLASTn function of BLAST+ version 2.6.0 (Camacho et al. 2009) with 1/−1 match/mismatch scoring parameters, and retaining the top BLAST hits for each locus.

STATISTICAL ANALYSES
All statistical analyses and figure construction in this study was performed with R version 3.4.3. For the unperfumed 4×4 mating dataset, we used Kruskal-Wallis tests to compare copulation success within 3 h of observation with male or female species identity (one test for each sex). In addition, we used Wilcoxon rank sum tests to perform post hoc pairwise comparisons for each sex to evaluate which species differed from one another in female acceptance of males or in male courtship success of females. We also calculated the reproductive isolation index RI 1 for prezygotic barriers as defined by Sobel and Chen (2014). This index is appropriate for comparisons between no-choice tests, and described by the equation RI 1 = 1 -(heterospecific success)/(conspecific success).
In 1×1 mating trials, mating rate was evaluated using a logistic regression with a chi-square test. For copulation latency, a two-way ANOVA was performed with latency in minutes as the dependent variable, to test for the effects of male identity, female identity, and their interaction. This was followed by a post hoc Tukey HSD in which we assessed which male identities (populations) were specifically different for copulation latency. For each of the three courtship behaviors (displays, tappings, and lickings), we converted the count data for each trial to a rate per unit time within a trial, by dividing the recorded counts for each trait by the latency to copulation (in minutes) or, if no copulation occurred, the maximum amount of observed time (180 minutes). Describing male behaviors in terms of rates accounts for differences among males in their courtship efficiency; for example, it allows us to differentiate males that performed fewer courtship behaviors because they were rapidly accepted by a female, from males that displayed lower courtship intensity across the total 3 h monitored period. For each of the three behavior rates (displayrate, tap-rate, lick-rate) as dependent variables, we performed an ANOVA with male identity, female identity, and their interaction as independent variables.
The perfumed 4×4 mating experiment was analyzed using planned contrasts that compared the mating success of males that had the same target identity but different (con-versus heterospecific) donor perfumes, when paired with SC D. americana females. This enabled us to specifically assess the effects of donor perfume variation on mating success of a given target male species. Each pairwise test was performed using a non-parametric Mann-Whitney U test.
Our primary analyses of differences in CHC composition were done after summarizing major axes of variation in CHC phenotypes using a set of Principal Component Analyses (PCA). Separate PCAs were performed on each CHC experiment within this study: one PCA for the initial (unperfumed) parental populations (denoted U-PCA), and one PCA for each of the perfumed experimental pairs-Nov-SC (N-PCA) and NE-SC (A-PCA; Figure S2). Within each dataset, a one-way ANOVA was used on each of the first 3 PCs to examine effects of sex and species (for the unperfumed dataset) or male perfume identity (in each perfuming study) on CHC composition. Additionally, for the perfumed datasets, T-tests were used to compare differences in PC values between target males that were paired with different (con-and hetero-specific) donor males to assess the effects of our perfuming manipulation on major axes of CHC variation. Factor loadings for N-PC an A-PC analyses, as well as individual compound differences in perfuming pairs can be found in the Supporting Information (Tables S10, S11, S12, and Figures S1,  S2).
For all gene expression analyses, expression is quantified in transcripts per million (TPM), and therefore is normalized within each sample. Here, we analyzed datasets separated by sex, as Davis and Moyle (2020) showed that the majority of genes have differential expression based on sex, whereas our primary interest here is in differences between species that might be implicated in female mating choices and male CHC profile variation. With the dataset for each sex, we ran one-way ANOVAs on TPM of every expressed locus (11301 genes total)-including our 23 target candidate genes-to determine loci for which gene expression varied by population. Each gene was ranked according to their resulting F-value (Table 3; uncorrected P-values are given in Table S13). This allowed us to evaluate which candidate genes had more pronounced expression differences between populations, compared to all other genes in the dataset. Candidate genes with greater differential expression between populations than the majority of the transcriptome could suggest these genes impact observed phenotypic differences, even if the number of tests performed make finding significance at alpha = 0.05 difficult.
We also estimated rates of nonsynonymous to synonymous substitutions (d N /d S ) for 4397 genes in this transcriptome-wide dataset, to quantify molecular evolution transcriptome-wide in this group and to evaluate evidence for positive selection specifically in candidate loci. We used a pipeline for obtaining genomewide estimates of d N /d S modified from Wu et al. 2018. First, as in Davis and Moyle 2020, transcript sequences from each population were aligned to the D. virilis reference genome (Flybase version 1.7, www.flybase.org) to identify loci with expression in the populations studied here. Then, for each gene with transcripts present in all populations and coding sequence (CDS) annotated in the D. virilis reference, a consensus fasta was generated for each line (population) for the longest splice variant of a given gene. These consensus fasta sequences were then aligned to each other using PRANK (Löytynoja and Goldman 2005) with codons enforced and ten bootstrap replicates, allowing us to obtain orthologous gene sequence alignments among the three populations used here and the D. virilis outgroup. We then calculated d N /d S from these aligned sequences in PAML version 4.9 (Yang 2007) using model M0 in CodeML-a maximum likelihood model for codon substitution. As PAML uses a tree-based model for computing d N /d S that is sensitive to use of the correct gene tree for a given gene (Mendes et al., 2016) consensus gene trees were constructed individually for each gene using RAxML v8.3 with the GTRGAMMA model with 100 bootstraps (Stamatakis 2014). Lastly, for both CG17821 and CG18609-genes with shared derived alleles for NE D. americana and D. novamexicana-we also computed the d N /d S value using only sequences for D. virilis and SC D. americana to determine if observed estimates of protein evolution in these genes is limited to the derived pair or if this pattern is consistent across the clade. Full results of d N /d S for each gene are reported in the supplementary file, with annotation of any known function of orthologs (according to Flybase) for the 100 loci with the highest estimated rates of protein evolution.
To evaluate possible introgression of candidate gene alleles between populations in this group, we constructed additional gene trees for each candidate gene using orthologs from the Anderson IN D. americana genome (from Kim et al. 2020) in place of our western (NE) D. americana sequences using RAxML version 8.3 (Stamatakis 2014), with alignments created using PRANK. Additional gene trees were generated in the same manner for neighboring genes within 50 kbp around CG17821/CG18609, to estimate size of a window of shared ancestry.
To assess evidence for genome-wide introgression, we used our transcriptome data in conjunction with the D. virilis genome to generate a genome-wide set of SNPs for our three focal populations and the D. virilis outgroup. To do so, we used bcftools (Li 2011) to generate a VCF (variant calling format) file and call SNP sites between the 4 taxa after filtering for low depth and masking heterozygous sites from individual taxa. With this, we calculated Patterson's D-statistic (Durand et al. 2011) using the Dtrios program from Dsuite (Malinsky et al. 2020) to calculate the three taxa D-statistic as well as an overall P-value using standard errors generated from the default 20 jackknife blocks. We used (((SC americana, NE americana), novamexicana), virilis)(figure 3B top tree) as the expected species tree topology for this analysis.

Supporting Information
Additional supporting information may be found online in the Supporting Information section at the end of the article. Table S1: Mating rate data from Table 3 of Spieth 1951.  Table S2: Posthoc comparisons for differences in male courtship success (based on proportion mated) between male populations regardless of female identity. Table S3: Posthoc comparisons of pairwise species differences in male behavior rates. Tukey HSD post-hoc mean differences and P-values between male identities for each behavior. Table S4: Posthoc comparisons of pairwise species differences for copulation latency based on male population identity. Tukey HSD post-hoc mean differences and P-values between male identities reported. Table S5: Loadings for principal component (U-PC) axes for each individual compound detected from GC/MS analysis of unmanipulated (unperfumed) samples. Table S6: Analyses of individual cuticular hydrocarbon compounds in unmanipulated (unperfumed) samples, from two-way ANOVAs testing for population or sex differences. Table S7: Analyses of the first three PCs of CHC compound variation from each perfuming pair, from one-way ANOVAs examining the effects of donor and target male identity on each PC. Table S8: Clade-wide d N /d S values for 23 candidate genes associated with CHC function. For genes with elevated clade-wide d N /d S , we evaluated d N /d S between D. virilis and SC D. americana only to determine if accelerated evolution is limited to the most derived pairs. Table S9: List of candidate genes and presence/absence of transcript expression in transcriptomes of each population used here (as determined by BLASTn hits to D. virilis ortholog). Table S10: Nov-SC perfume sample loadings for principal component (N-PC) axes for each individual compound detected from GC/MS analysis. Table S11: NE-SC perfume sample loadings for principal component (A-PC) axes for each individual compound detected from GC/MS analysis. Table S12: Analyses of individual cuticular hydrocarbon compounds in perfumed samples, from t-tests comparing target males with con-and heterospecific donor perfumes. For example, D. novamexicana C21:1 tests the difference between C21:1 values in D. novamexicana target males perfumed with D. novamexicana or SC D. americana donor males. Compound abundance values were log-transformed for analysis. Table S13: Species differences in quantitative gene expression for 23 candidate genes associated with CHC function, analyzed separately by sex, with F-value and un-adjusted P-value of ANOVA test. Figure S1: Relative log-scale abundance of compounds for perfumed males from NE-SC (left) and Nov-SC (right) experimental pairs for each major compound type (line = mean, box = central quartile, whiskers = S.E, n = 4 for all bars). Figure S2: PC1 and PC2 axes showing individual samples for each male perfuming identity for the Nov-SC pair (panel A) and the NE-SC pair (panel B). Figure S3: Distribution of clade-wide d N /d S values from 4397 genes. The median d N /d S value was 0.065, while the mean and standard deviation were 0.0991 ± 0.133. Supporting Information