Divergent evolutionary trajectories shape the postmating transcriptional profiles of conspecifically and heterospecifically mated cactophilic Drosophila females

Abstract

Postmating-prezygotic (PMPZ) reproductive isolation is hypothesized to result from divergent coevolutionary trajectories of sexual selection and/or sexual conflict in isolated populations. However, the genetic basis of PMPZ incompatibilities between species is poorly understood. Here, we use a comparative framework to compare global gene expression in con- and heterospecifically mated Drosophila mojavensis and D. arizonae female reproductive tracts. We find striking divergence between the species in the female postmating transcriptional response to conspecific mating, including differences in differential expression (DE), alternative splicing (AS), and intron retention (IR). As predicted, heterospecific matings produce disrupted transcriptional profiles, but the overall patterns of misregulation are different between the reciprocal crosses. Moreover, we find a positive correlation between postmating transcriptional divergence between species and levels of transcriptional disruption in heterospecific crosses. This result indicates that mating responsive genes that have diverged more in expression also have more disrupted transcriptional profiles in heterospecifically mated females. Overall, our results provide insights into the evolution of PMPZ isolation and lay the foundation for future studies aimed at identifying specific genes involved in PMPZ incompatibilities and the evolutionary forces that have contributed to their divergence in closely related species.

Fig. 1. Evidence of postmating-prezygotic isolation between D. mojavensis and D. arizonae.

a Egg production over seven days of D. arizonae females following a con- or heterospecific mating (n = 10 and 11 females, respectively). b Proportion of eggs oviposited by D. arizonae females in a, that hatch into first instar larva. c Proportion of eggs oviposited by D. arizonae females following a con- or heterospecific mating that were determined to be fertilized using DAPI staining (n = 445 and 338 eggs, respectively). These results indicate an overall decrease in reproductive success of heterospecifically mated D. arizonae. d Levels of sperm motility within seminal receptacles in both D. mojavensis and D. arizonae females mated con- and heterospecifically at 1 and 5 days postmating. P-values of correlations are noted: *P < 0.05, **P < 0.01, ***P < 0.001, ns not significant (P > 0.05). Boxplots represent the median with 25th and 75th percentiles, and whiskers show the 1.5 interquartile range. Bar graphs identify mean and standard error bars.

Fig. 2. Experimental design used for differential expression and alternative splicing as a response to conspecific (blue) and heterospecific (yellow) matings in D. mojavensis and D. arizonae.

RNA-seq libraries were constructed from LRT of virgins, con- and heterospecifically mated females at 45 min and 6 h postmating. All comparisons were performed between virgin females and mated females for each species.

Fig. 3. Conspecific postmating transcriptional response in D. mojavensis and D. arizonae.

a Number of conspecific mating responsive genes with significant patterns of DE, AS and IR following conspecific matings in D. mojavensis and D. arizonae. All comparisons were performed against virgin females at 45 min and 6 h postmating (FDR_α = 0.05). Genes in the DE-AS-IR category, showed significance in two or more of the individual categories (DE, AS and/or IR). Venn diagrams comparing significant b DE genes, c AS genes, and d IR genes from conspecific matings in each species. Overlapping genes represent the conserved postmating response, while species-specific genes indicate divergence in the postmating response. The distinct expression and splicing patterns indicate substantial differences in the postmating transcriptional response between the species.

Fig. 4. Transcriptional correlations between conspecific mating responsive DE genes in D. mojavensis vs D. arizonae.

Scatterplots represent the relationship between relative fold changes (log₂) for conspecific matings in D. mojavensis vs D. arizonae at 45 min and 6 h postmating. Log₂ fold changes are relative to virgin females. All conspecific mating responsive genes in each species are included (FDR_α = 0.05). Pearson’s R² correlation coefficients and linear method trend-lines (with 95% confidence intervals shaded) between the species are indicated. P-values of correlations are noted: ***P < 0.001.

Fig. 5. Gene ontology analysis of DE genes.

Functional analysis for significant DE genes, indicating distinct biological process gene ontology enrichment categories between the species for genes differentially regulated relative to virgins in con- and heterospecific crosses. The fold enrichment of detected genes within each enriched category is indicated.

Fig. 6. Comparison of differentially regulated genes in conspecific and heterospecific crosses in D. mojavensis and D. arizonae.

Overlapping genes represent conspecific mating responsive genes that were properly regulated in heterospecific crosses. Genes unique to the conspecific cross represent conspecific mating responsive genes that were misregulated in heterospecific crosses. Genes unique to the heterospecific cross are those that are not part of the normal conspecific mating response but had disrupted transcriptional profiles in the heterospecific cross (“heterospecific only” genes).

Fig. 7. Transcriptional correlations for DE genes between con- vs heterospecific matings in D. mojavensis and D. arizonae.

Scatterplots represent the relationship of relative fold changes (log₂) between con- vs heterospecific matings at 45 min and 6 h postmating in a D. mojavensis and b D. arizonae. Log₂ fold changes are relative to virgin females. Genes colored in blue are all conspecific mating responsive genes. Genes colored in yellow are those that were DE in the heterospecific cross but not the conspecific cross (“heterospecific only”). Pairwise Pearson’s R² correlation coefficients and linear method trend-lines (with 95% confidence intervals shaded) are shown. The black regression line is fitted through all the data (blue and yellow combined). P-values of correlations are noted: ***P < 0.001.

Fig. 8. Relationship between postmating conspecific expression divergence and the level of heterospecific disruption for conspecific-responsive DE genes.

Scatterplots show the relationship between postmating conspecific divergence (e.g., D. mojavensis conspecific expression divergence = |conspecific D. arizonae−conspecific D. mojavensis | ) and the level of heterospecific disruption (e.g., heterospecific disruption in D. mojavensis = |heterospecific D. mojavensis−conspecific D. mojavensis | ) for a D. mojavensis and b D. arizonae. All DE genes detected in con- or heterospecific crosses at 45 min and 6 h postmating are included (FDR_α = 0.05). Genes colored in blue are all conspecific mating responsive genes. Genes colored in yellow are those that were DE in the heterospecific cross but not the conspecific cross (“heterospecific only”). Pairwise Pearson’s R² correlation coefficients and linear method trend-lines (with 95% confidence intervals shaded) are shown. The black regression line is fitted through all the data (blue and yellow combined). P-values of correlations are noted: **P < 0.01, ***P < 0.001, ns not significant.

Fig. 9. Median pairwise ω (dN/dS) for genes that were misregulated in heterospecific crosses.

Blue boxes represent mating responsive genes in conspecific crosses that were misregulated in heterospecific crosses. Yellow boxes represent genes that were differentially expressed only in heterospecific crosses. The median ω for the genome background (blue dash line) is indicated. Significant comparisons against the genome background using the Dunn method for joint ranking are indicated. P-values are noted: ***P < 0.001. Boxplots represent the median with 25th and 75th percentiles, and whiskers show the 1.5 interquartile range.