Genomic responses to the socio-sexual environment in male Drosophila melanogaster exposed to conspecific rivals.

Socio-sexual environments have profound effects on fitness. Local sex ratios can alter the threat of sexual competition, to which males respond via plasticity in reproductive behaviors and ejaculate composition. In Drosophila melanogaster, males detect the presence of conspecific, same-sex mating rivals prior to mating using multiple, redundant sensory cues. Males that respond to rivals gain significant fitness benefits by altering mating duration and ejaculate composition. Here we investigated the underlying genome-wide changes involved. We used RNA-seq to analyze male transcriptomic responses 2, 26, and 50 h after exposure to rivals, a time period that was previously identified as encompassing the major facets of male responses to rivals. The results showed a strong early activation of multiple sensory genes in the head-thorax (HT), prior to the expression of any phenotypic differences. This gene expression response was reduced by 26 h, at the time of maximum phenotypic change, and shut off by 50 h. In the abdomen (A), fewer genes changed in expression and gene expression responses appeared to increase over time. The results also suggested that different sets of functionally equivalent genes might be activated in different replicates. This could represent a mechanism by which robustness is conferred upon highly plastic traits. Overall, our study reveals that mRNA-seq can identify subtle genomic signatures characteristic of flexible behavioral phenotypes.

INTRODUCTION

The ability of individuals to modify their phenotypes in response to the prevailing environment is a crucial determinant of fitness (West-Eberhard 2003). Phenotypic plasticity is widespread, from wholesale developmental switches, such as temperature-dependent sex determination in reptiles (Janzen and Phillips 2006), to flexible, reversible behaviors, such as predator-dependent shoaling in fish (Tien et al. 2004). Flexible plasticity is expected to be of particular importance in highly variable contexts (Komers 1997) because it offers the potential to track the prevailing environment (Bretman et al. 2012).

Highly flexible plasticity is frequently observed in male mating strategies. Across species, males exhibit diverse types of plasticity in response to predicted levels of male–male competition (Parker et al. 1996, 1997). For example, when male–male competition is strong and ejaculates are in limited supply (Hihara 1981; Linklater et al. 2007), males that exhibit plastic ejaculate allocation should be favored (Gage 1991; Wedell et al. 2002). Male–male contests in the form of sperm competition can be separated into “intensity” (the average number of males competing for a given set of eggs) and “risk” (the probability that females in the population have mated or will mate again) (Parker 1990a,b, 1993; Parker et al. 1996, 1997; Wedell et al. 2002; Engqvist and Reinhold 2005). Males can alter their investment in each reproductive bout in response to the likelihood of male–male competition as signaled by a female's mating status (mated versus virgin), by the presence of potential rivals or by assessing the condition of those rivals (Parker 1982; Parker et al. 1996, 1997; Wedell et al. 2002; Bretman et al. 2011a). Altering the number of rival males present has the potential to alter both intensity and risk, and both can also vary with respect to the average levels in a population or species, or to their current levels in a mating bout (Engqvist and Reinhold 2005). Hence, the optimal investment patterns according to chronic or acute exposure to rivals can vary (Parker et al. 1996, 1997).

The responses of males to the presence of same sex conspecifics can include differential investment in reproductive tissue during development (Kasumovic and Brooks 2011), strategic ejaculation of sperm (Wedell et al. 2002), and behavioral plasticity (Bretman et al. 2011a). In D. melanogaster, males increase mating duration (Bretman et al. 2009) and modulate ejaculate composition following exposure to a conspecific potential rival male (Wigby et al. 2009). These responses increase a male's fitness through siring increased numbers of offspring (Bretman et al. 2009). This benefit is associated with an increase in the transfer of seminal fluid components (Wigby et al. 2009) and sperm (Garbaczewska et al. 2013; Moatt et al. 2014). However, there are costs for males that exhibit sustained responses to rivals, evident as reduced lifespan and reduced late-life mating capacity (Bretman et al. 2013).

The duration of mating following exposure of males to conspecific rivals is associated with the length of exposure, rather than to the density of males encountered (Bretman et al. 2010). It takes 18–24 h of exposure to a rival to achieve a significant increase in mating duration (Bretman et al. 2010), with a similar temporal decay in the response following rival removal (Bretman et al. 2012; Rouse and Bretman 2016). The magnitude of the response to rivals continues to increase for at least 5 d of exposure (Bretman et al. 2010). Hence a window of 0–2 d is ideal in order to capture the initiation and establishment of this response. It is this time interval upon which we focussed in this study.

D. melanogaster males use multiple, redundant cues to detect conspecific rivals (Bretman et al. 2011b; Kim et al. 2012). A combination of any two cues from smell, touch, and sound (Bretman et al. 2011b) is required for males to respond. The alternative sensory modalities used suggest that a single initial stimulus (exposure to rivals) may result in the same outcome (increased fitness) via the induction of combinations of alternative sensory pathways leading to physiological responses in behavior (mating duration) and in ejaculate composition. The multifaceted responses of males to their rivals could be mediated by a number of different mechanisms, including via transcriptomic, epigenetic, or proteomic changes. Such changes could result in alterations to cell–cell communication via hormones, neurotransmitters, or other signaling molecules.

In this study, we tested the hypothesis that the responses of D. melanogaster males to their conspecific rivals are, at least in part, mediated by a rich and complex signature of phenotypic plasticity in the male transcriptome. The prediction of observable transcriptional changes in response to the social and sexual environment is supported by previous research. For example, short-term (<20 min) exposure of males to females or other males can alter the transcription of spermatogenesis and odor perception genes (Carney 2007; Ellis and Carney 2011). Males exposed to other males for up to 72 h also show differential expression (DE) in two of three seminal fluid protein genes (but not in four testis genes) tested (Fedorka et al. 2011). However, no study has yet determined whole-transcriptome responses that underlie the adaptive plastic responses of males to their rivals, nor tested for divergent responses in different body parts. The latter is crucial to resolve signatures of tissue-specific responses and avoid “swamping” of subtle signals (Chintapalli et al. 2007). Our design to profile gene expression responses to rivals for 0–2 d offered the opportunity for comparisons with these previous studies (Carney 2007; Ellis and Carney 2011).

A full understanding of how highly flexible plasticity evolves requires a detailed and genome-wide investigation of the underlying molecular and cellular mechanisms (Schlichting and Smith 2002; Hobert 2003; Aubin-Horth and Renn 2009; Ellers and Stuefer 2010). Ultimately, this is essential to answer fundamental questions such as how rapid, transient and reversible behavioral plasticity can be achieved, how the genome might achieve robust and reliable responses against a background of variable behavior, and what proportion of the transcriptome is involved. Such findings will also show how phenotypic plasticity itself is likely to evolve in the face of rapidly changing environments. To date, we lack descriptions of how phenotypic plasticity is mediated via changes in the transcriptome, epigenome, and proteome (Gilbert 2005; Aubin-Horth and Renn 2009; Ellers and Stuefer 2010; Zhou et al. 2012). In part, this may be because behavior is inherently complex, often with an underlying polygenic basis (Bell and Aubin-Horth 2010). Behavioral plasticity can also be highly flexible and reversible—hence gene expression changes underlying behavioral phenotypes may be subtle, transient, and potentially close to current detection thresholds (Carney 2007).

In this study, we addressed the omissions noted above by measuring the genome-wide transcriptional responses of male D. melanogaster exposed only to conspecific rival males from 2 h to 2 d. We profiled the gene expression changes in two different body parts, head–thorax (HT) and abdomen (A), using mRNA-seq in order to enhance the sensitivity of the analyses and to enable us to detect unknown transcripts and expression variants (Marioni et al. 2008; Sultan et al. 2008).

RESULTS

We conducted RNA sequencing on three replicated pools of males exposed to a conspecific rival or maintained singly, and we sampled males following 2, 26, or 50 h of exposure. The samples were taken prior to mating, in the absence of females. The aim was to capture a broad sweep across the phenotypic expression of responses to conspecific rivals. Based on previous research by us and others, we chose the 2 h time point to capture initial gene expression responses prior to the expression of any phenotypic responses, 26 h as the point at which we know males have passed a threshold of significant phenotypic responses to rivals, and 50 h to capture any longer-term changes at a point where the phenotypic responses appear to have ceased (Carney 2007; Bretman et al. 2009, 2010; Ellis and Carney 2011; Fedorka et al. 2011). An additional, important justification for these three time points was the need to avoid strong interactions with circadian effects on gene expression (by sampling cohorts of flies for gene expression analysis at the same absolute time of day, i.e., after 2 h, then 24 h later at 26 h of exposure, and finally a further 24 h later at 50 h of exposure). We divided the flies into head/thorax (HT) and abdomen (A) to identify body part-specific facets of the males’ responses to rivals—specifically to capture the response of the nervous system in the HT and of the reproductive system in the A.

Bioinformatics

The RNA-seq data exhibited good values on a range of quality checks and provided a sensitive quantitative resolution of DE. First, for all positions across the reads, the fastq quality scores were high and showed little variability. Hence there was no evidence for any sequencing errors. Nucleotide distributions (Supplemental Fig. S1) were consistent across all samples. We observed variation in read numbers (35 M–129 M reads per sample), resulting in wide variation in complexity (Supplemental Table S1), defined as the ratio of nonredundant (unique) to redundant (all) reads. The proportion of redundant (R) genome matching reads was ∼71% and for nonredundant (NR) reads ∼50% (Supplemental Fig. S2). The majority of reads mapped to CDS and 5′ and 3′ UTR, indicating high-quality reads in the data set and only minor contamination by other RNA species (Supplemental Fig. S2). For the analysis of expression level variation we retained exon–exon boundary reads and excluded 3′ and 5′ UTRs (∼10% of reads mapped to 3′ UTRs and ∼5% to 5′ UTRs). Reads incident to ncRNAs, such as t/rRNAs, were also excluded. The mRNA sequencing was nondirectional. There was a consistent positive strand bias in the sequencing data, but no systematic variation across samples, body parts, or time periods. Jaccard similarity comparisons (Levandowsky and Winter 1971) highlighted the effect of the variation in read numbers and complexity across samples. We increased replicate and sample comparability by applying a subsampling normalization (Efron and Tibshirani 1993; Gierliński et al. 2015; Schurch et al. 2016) without replacement (I Mohorianu, A Bretman, DT Smith, EK Fowler, T Dalmay, T Chapman, in prep.) at a fixed total (50 M, Supplemental Table S2A). A random subsample was selected for analysis, following a check to ensure that subsamples were representative of the original (example in Supplemental Table S2B). The normalization was highly effective (Supplemental Figs. S3, S4). The remaining minor variation was adjusted using quantile normalization (Supplemental Fig. S4; Bolstad et al. 2003). The efficiency of the normalization was tested using point-to-point Pearson correlation coefficients, calculated for each gene, on the expression profile across the length of the gene. This revealed tight correlations in expression level between the original and subsampled data (Supplemental Fig. S5). However, some replicates were significantly different in terms of per-gene complexity and distributions of point-to-point Pearson correlations, and this prompted the removal of outlier replicates and the identification of the final data set for analysis of two out of three of the original replicates for each treatment (Supplemental Table S2A).

Identification of genes showing DE

We applied a two-step procedure for the DE call using the two levels in the experiment: body parts (HT, A) and treatments (presence or absence of rivals). Body part DE distribution revealed a much higher number of DE transcripts than for ±rivals. We used an offset fold change (OFC), with an empirically determined offset of 20 for the DE call, to filter out low abundance genes. The distributions of replicate-to-replicate differences were tightly centered around zero, as expected (Supplemental Figs. S6, S7). However, for each body part, it was characteristically the few, low abundance genes exhibiting DE that contributed most to replicate-to-replicate differences (Supplemental Fig. S7). These were “leaky” genes, e.g., genes showing DE in a body part in which they are hardly, if at all, expressed (Chintapalli et al. 2007). This is likely to have been due to inaccuracy in the separation of, or movement of abundant mRNAs between, body parts. Hence, we removed such genes showing DE between replicates in the “wrong” body part (Supplemental Fig. S6B).

A summary of the number of DE genes revealed that far more genes showed DE in the HT than in the A (Fig. 1; Supplemental Table S3). There was a marked early up-regulation in the HT in response to rivals, which was shut down by 50 h. The number of genes showing DE across both replicates in the A was lower, but in contrast to the HT increased over time. The functions of the genes showing DE are summarized below.

FIGURE 1.

Summary of the number of differentially expressed genes. The number of genes showing consistent differential expression (>|log2 1.5| ∼ 0.5 OFC) across replicates in the (A) HT and (B) A body parts, following 2, 26, and 50 h of exposure of males to no rivals versus rivals. Up-regulated (+ rivals > no rivals) or down-regulated (+rivals < no rivals). The number of up (down)-regulated genes in the HT at 2, 26, and 50 h, respectively, was 325 (10), 8 (59), 13 (19) and for the A was 4 (1), 4 (8), and 4 (12).

HeadThorax

Exposure to rivals caused a rapid, early up-regulation of 325 genes, many of which were associated with the detection of sensory stimuli (e.g., Gustatory receptor 10a) odorant or pheromone binding (e.g., odorant receptor 1a, 2a, 10a, 13a, 47a, 63a, Odorant-binding protein 83g), hormone metabolism (e.g., shadow, Ecdysone-inducible L1, ETHR), proteolysis (e.g., CG13748, CG6763, CG8563, stall, CG10073, Matrix metalloproteinase 2, CG1632), male courtship behavior (e.g., hector), GPCR (e.g., Serotonin receptor 1A, CG31714, CCK-like receptor at 17D3, rickets, CG32547, Octopamine-Tyramine receptor, SIFamide receptor, Leucokinin receptor, CG7918, CG33696, mangetout, Diuretic hormone 44 receptor 1, Adenosine receptor, CG33639), and cuticle/chitin metabolism (e.g., Lcp65Ac, Acp65Aa, Cuticular protein 64Aa, 64Ad, 76Bd, CG8192, obstructor-B, obstructor-H, CG34220, Ecdysone-dependent gene 78E, CG13837, Cht7, Ccp84Ag, CG13676, CG33263, CG13806)—summarized in Figure 1, Supplemental Table S3, and Supplemental Material 1. The 10 down-regulated genes at 2 h were predominantly pheromone and chemosensory genes (e.g., Olfactory-specific E, Os-C, Chemosensory protein B 42a and 93a, Chemosensory protein A 87a and 86a). By 26 h, the number of up-regulated genes dropped to eight (including a peptidase [E(spl) region transcript m1] and a chitin gene [Chitin deacetylase-like 9]), and down-regulated genes had increased to 59 and included a strong response in hormone metabolism (e.g., Ecdysone-dependent gene 78E), olfaction/chemosensation (e.g., Odorant receptor 46a, Chemosensory protein A 7a, 56a, and 87a, Chemosensory protein B 42a, 53a, 53b, and 93a, Odorant-binding protein 19b and 57e) as well as chitin metabolism (e.g., Cuticular protein 49Af, CG42729), proteolysis (e.g., CG18478, CG15369, CG3117, CG11459, CG30091, CG30098, CG31827, CG33458, CG43335, CG43336), and immunity (e.g., Turandot A, X, and C). By 50 h the number of genes showing up-regulation remained low at 13 (including proteolysis genes such as CG18557, CG33459). The number of down-regulated genes had dropped to 19 and included odorant binding/pheromone genes (such as Os-C, antennal protein 5, 10, Odorant-binding protein 59a, 56b).

Abdomen

The A showed a distinctly different profile of responses, with many fewer genes showing DE (Fig. 1; Supplemental Table S3; Supplemental Material 1). Few genes were consistently up-regulated across time (four genes at each time point), and the number of down-regulated genes increased over time (1, 8, and then 12, at 2, 26, and 50 h, respectively). At 2 h, only a few genes showed consistent DE and this included significant up-regulation in immune genes (Diptericin, Diptericin B, and Attacin-C). By 26 h, the eight down-regulated and four up-regulated genes were mostly of unknown function but included down-regulation in Niemann-Pick type C-1b. We noted that two seminal fluid proteins (Sfps) up-regulated in one replicate at 2 h were down-regulated in the same replicate by 26 h (Sfp79B, Sfp70A4). At 50 h, there was significant down-regulation in seminal fluid protein Sfp79B and ecdysone-related 71Eb. Hence, Sfp79B was up and then consistently down-regulated across 2–50 h. Notable was the observation that there were many changes in Sfp gene expression, but this was mostly inconsistent across replicates (Supplemental Table S3), with many Sfps missing the statistical threshold applied. For example, there was inconsistent up-regulation of Acp26Aa, Acp70A, Acp98AB, Sfps 79B, 51E, 70A4, and 96F at 2 h followed by inconsistent down-regulation at 26 h for Sfps 70, 79, Acp76A, and Acp36DE. By 50 h there was inconsistent down-regulation in Acp24A4, Sfps24Ba 24Bb, 35C, 60F, 24Bc, 33A4, Peb II.

Analysis of functional enrichment among genes showing DE

To complement the above analysis, we first investigated the number of genes showing DE that fell into the major functional categories represented (Table 1; Supplemental Tables S4, S5; see also Supplemental Material 1). To test whether such categories showed significant enrichment within our data, we then compared, for each keyword, the observed frequency density distributions (rivals–no rivals expression) of all genes within each keyword class to a control random distribution of the same number of genes (example in Fig. 2 for GPCR genes in the HT and immune genes in the A; Supplemental Tables S4, S5). This analysis showed that olfactory genes as a whole exhibited dynamic changes in expression over time, initially showing simultaneous increased and decreased gene expression, relative to the control random distribution, at 2 h (Supplemental Table S5) and consistently decreased expression in response to rivals thereafter. The A pattern for olfactory genes was similar. HT pheromone genes showed consistent, significant down-regulation in response to rivals. Sensory genes within the “trichogen” term showed evidence for subtle early up-regulation in the HT, followed by down-regulation at 26 h and then both increased and decreased gene expression, relative to the control random distribution, at 50 h. Chitin HT genes showed a significant increase in expression in response to rivals at 2 h and significantly increased and decreased gene expression, relative to the control, at both 26 and 50 h. The same was true for cuticle genes. Consistent with increased signaling by sensory systems, there was a strong and significant up-regulation of GPCR genes in response to rivals in the HT at 2 h. Most of the response in seminal fluid (Sfp) genes was inconsistent across replicates due to temporal divergence—a large increase in expression at 2 h in response to rivals in the A, followed by significantly decreased expression in one replicate at 26 h and in the other at 50 h (Supplemental Table S6). Many proteolysis genes responded to the presence of rivals. In the HT they showed a slight increase, then decrease in expression in response to rivals over time. There was a similar, though less marked trend, in the A. Immunity genes showed a similar pattern to the proteolysis class.

TABLE 1.

Summary of the patterns of DE in the analysis of functionally enriched gene categories in the HeadThorax (HT)

FIGURE 2.

Differentially expressed genes in two example functionally enriched gene categories. Examples of the frequency density of expression levels for rivals–no rivals for two functional keyword gene groups: (A) G-protein coupled receptor “GPCR” for the HT (94 genes), and (B) “Immune” for the abdomen (Abd) (43 genes). A shift in the positive direction indicates higher expression in the presence of rivals, a shift to the left lower expression. In pink are all the genes corresponding to the keyword category and in green a random control distribution of the same number of genes drawn from the total set of 15,513 genes. For each keyword, the frequency density of gene expression for both replicates is shown for the 2, 26, and 50 h time periods.

Taken together, the functional categories analysis validated the DE call of individual genes and provided evidence for consistent and significant gene expression changes in response to rivals, with the early induction of sensory genes and up-regulation of structural ejaculate components. Immune response genes in both body parts also showed an early increase then decrease in expression over time (summary Fig. 3). An interesting additional finding was that consistent changes within keyword functional categories seemed sometimes to be achieved using different underlying sets of functionally similar genes within each replicate. This suggested that responses to rivals might sometimes be achieved via the activation of alternative, potentially functionally equivalent, underlying pathways. An example was the response of immune genes in the A, where there was consistent up-regulation and then down-regulation of immune genes in both replicates over time though the underlying gene identities were different in each replicate. We conducted an initial exploration of alternative responses by counting the number of times that genes showed inconsistent DE among replicates within keyword categories (Supplemental Table S6). This revealed potential alternative responses in both body parts. At 2 h in the HT, the highest alternative responses were observed in proteolysis and the lowest in GPCR, immune, and olfactory genes. Interestingly, by 26 h alternative responses in most keyword classes had increased and this continued through to 50 h. The number of genes available for the analysis in the A was low. Among categories with the most genes, we observed immune and proteolysis genes as showing alternative responses across time.

FIGURE 3.

Summary model for the responses of males to rivals. Exposure to rivals initiated a rapid early increase at 2 h in the expression of sensory genes, which was shut down by 26 h and then off by 50 h. The pattern in the A was contrasting, with far fewer genes involved and the number of genes changing in expression increased, rather than decreased, over time. In the A there was an early up-regulation in immune genes and by 50 h, changes in seminal fluid protein and hormone genes. Many changes in seminal fluid protein gene expression were inconsistent, and each replicate appeared to switch on different sets of ejaculate component genes.

Variation within mRNA pools

We checked the assumption of low variability within mRNA pools of 40 individuals from which each sample of RNA was extracted. Our analysis showed very low variation within pools overall and only a few alleles were differentially expressed (Supplemental Fig. S8). In the A these included genes with roles in male reproduction, immunity, cuticular/chitin, and odorant/taste proteins. In the HT the differentially expressed alleles included male reproduction, immunity, cuticle/chitin, and mucin genes. As an example, we analyzed the transcript variants for the sex peptide (Acp70A). Two major coding sequence polymorphisms were evident—G to T (Alanine to Serine) in exon 1 (polymorphism c55 of Chow et al. 2010) and T to A (Asparagine to Lysine) in exon 2 (c132 of Chow et al. 2010). The former polymorphism is segregating in the Dahomey wild type used in these experiments (DT Smith, R Evans-Gowing, T Chapman, unpubl.). Expression level variation is reported for these polymorphisms (Chow et al. 2010) associated with several sex peptide-mediated traits.

Validation by qRT-PCR

Variation in gene expression was validated using low-throughput qRT-PCR. We selected three reference genes (that showed low variation across both body parts and time points in our mRNA-seq data). There was good agreement between the overall expression levels and pattern of differences between rivals and no rivals treatments in the RNA-seq versus qRT-PCR data (Fig. 4; Supplemental Fig. S9; Supplemental Table S7). Many different functional categories that formed the DE response were represented among the validated genes, such as multiple seminal fluid, cuticular, odorant binding, and immune genes. To measure concordance, we counted the number of times DE in mRNA-seq data was also observed in the qRT-PCR, and the number of times a lack of DE was observed. Using this analysis, 78% of loci across both replicates were validated (Supplemental Table S7).

FIGURE 4.

qRT-PCR validation of mRNA-seq gene expression levels. Variation in expression level in the qRT-PCR and the corresponding mRNA-seq sequencing expression levels are shown for four example loci (no rivals, light gray; rivals, dark gray); Fbgn0015586 (Acp76A), Fbgn0036110 (Cpr67Fb), Fbgn0000276 (CecA), and Fbgn0028396 (TotA). Shown are the normalized expression levels (RNA-seq) or normalized ΔCt expression levels (qRT-PCR) for A (Acp76A, Cpr67Fb, CecA) and HT (TotA) body parts. Sample labels: 02, 26, or 50 h of exposure, rep 1 or 2. Error bars show ± 10% expression level.

DISCUSSION

The major result was the identification and characterization of the complex genomic responses of D. melanogaster males to their same sex conspecific rivals, prior to mating and in the absence of females. Clear evidence for DE was observed in genes involved in assessing the likely level of sperm competition (sensory genes) and in producing a response to that environment (ejaculate component genes). There was a robust early activation of sensory pathways in the head thorax (HT) and an inconsistent up-regulation of ejaculate component genes in the abdomen (A). Our results have implications both for our understanding of the evolution of phenotypic plasticity and more broadly for the analysis of variation within transcriptomic data.

DE in different classes of genes was evident in the genome-wide responses to rivals. The response to rivals was initiated early by up-regulation of a large number of sensory perception genes in the HT, which then decayed over time. There was a much smaller number of genes showing DE in the A and this appeared to increase over time. The primary response of males to rivals in the HT was a rapid change to the expression of genes involved in sensory perception and chemical communication, particularly odorant signaling. The consistent changes to trichogen cell, GPCR and chitin metabolism are consistent with a change in sensory perception. Trichogen sensory hair cells form an integral part of the sense organs of the peripheral nervous system. They are one of four types of cells that derive from sensory organ precursors and which form mechanosensory, chemosensory, or auditory bristles (Keil 1997) in the peripheral nervous system. This is consistent with the known redundant role of touch/taste, odorants, and sound in the responses of males to rivals (Bretman et al. 2011b). We found signatures of this redundancy in sensory modalities in the divergent types of trichogen cell genes that were activated. Similarly, GPCRs can respond to diverse sensory ligands, including many odorants, light sensitive compounds, pheromones, and neurotransmitters (Brody and Cravchik 2000; Hewes and Taghert 2001). In line with other researchers, we assume that the levels of GPCR mRNAs give an indication of the activation of GPCR proteins (e.g., Churcher et al. 2015; Thulasitha et al. 2015; Quandt et al. 2016), though it is also possible that GPCR genes could be transcriptionally activated by a pathway not involving the same GPCR. There was a strong and consistent profile of early up-regulation followed by later down-regulation and then dampening of expression, again consistent with the activation of alternative sensory modalities. As with the trichogen cell, there was evidence for the activation of GPCR signaling in genes within this class representing all the major sensory modalities. There was also a change in the expression of photoreception and in GPCR genes in the HT at 2 h in response to rivals. This is concordant with the minor role of vision in the response of males to rivals and showed that this facet of the response occurred early. The role of the many genes involved in chitin metabolism in the response to rivals is less clear.

The more limited gene expression response in the A was DE in immune genes and in seminal fluid and protease genes that synthesize and process components of the ejaculate. This is consistent with the qualitative and quantitative changes to ejaculate delivery observed upon exposure of males to rivals (Wigby et al. 2009; Sirot et al. 2011). For two seminal fluid protein genes (Acp26Aa, Acp62F) reported to be down-regulated over time in the presence of rivals (Fedorka et al. 2011), we also found lower gene expression with rivals. We found no consistent pattern of DE in Acp70A in response to rivals, in accordance with previous results (Fedorka et al. 2011). Ejaculate protein genes (Findlay et al. 2008) were found among differentially expressed genes in the A samples, and a few also in the HT, suggesting that some seminal fluid protein genes have functions in and outside the male reproductive tract. The absence of any major signal from sperm genes was interesting as the number of sperm transferred to females can respond to the level of sperm competition (Price et al. 2012; Garbaczewska et al. 2013). However, as reported in Fedorka et al. (2011), we found no evidence for DE in four sperm genes following male exposure, over a similar time course. The time scale for the production of sperm from start to finish is much longer (i.e., many days [Fabian and Brill 2012]) than for the replenishment of seminal fluid components. Therefore, we suggest that males hold a large store of mature sperm that they can apportion strategically, but which takes some time to replenish. However, seminal fluid proteins are more easily exhausted (Linklater et al. 2007), which may select for mechanisms that initiate rapid activation of seminal fluid protein genes when rivals are present.

DE in immune genes was observed in both HT and A samples, which is consistent with findings in humans suggesting that social contact can alter gene expression (Tung and Gilad 2013) and in particular in immune-related genes (Cole 2013). The extent to which immune genes are involved in the assessment of sperm competition threat or in associated changes in behavior or ejaculate investment are not yet known. However, the data are consistent with the emerging idea that immune genes may also respond to different social environments as part of a generalized stress response or to modify the likelihood of disease transmission (Perkins et al. 2009).

In addition to validation via qRT-PCR, external, independent validation of the expression data was also provided by the overlap of genes showing DE with those identified in several previous studies (Carney 2007; Findlay et al. 2008; Ellis and Carney 2011; Fedorka et al. 2011; Supplemental Material 2; Supplemental Table S8). There were also some interesting differences. For example, genes reported to change expression in males upon short-term initiation of courtship (Carney 2007) showed low overlap with those involved in responding to conspecific rivals, as tested here. There was some overlap between socially responsive genes down-regulated following 20-min exposure to males or females in Ellis and Carney (2011) and the genes that were up-regulated at 2 h in our study. This could suggest that early sex-specific responses (Ellis and Carney 2011) in gene expression may be recruited to rival responses. This prompts future investigations of the specificity of genes involved in mediating social interactions.

The data suggested the hypothesis that alternative responses could contribute to the overall response to rivals. Support for this came from the findings that different replicates sometimes responded to the same input stimulus via different pathways to achieve a functionally equivalent output (within replicate variation was also very low). We suggest that alternative responses might represent a biological feature of these data. However, we recognize that further direct tests of this idea are needed, e.g., by using inbred lines and/or transcriptomes from individual males. A prediction is that the transcript profiles of individual males with limited sensory inputs available for detecting rivals would be narrower than for the unmanipulated males tested here. The existence and benefits of redundant pathways and networks is increasingly realized in the study of neurobiology (Greenspan 2012) and might confer robustness and increase the stability of important, fitness-related outputs for individuals responding to complex and variable environments (Bretman et al. 2011b). Alternative pathways underlying the same output are consistent with key elements of the biological response to rivals. Specifically, males can detect rivals using any two of three sensory modalities from smell, touch, and hearing (song). However, as far as is currently known, which of the modalities is used does not matter because the outcomes are functionally equivalent (Bretman et al. 2011b). Different replicates of males may have used different sensory pathways to detect rivals. Alternative Sfp gene responses could be facilitated by the functional redundancy observed in seminal fluid composition (Ravi Ram and Wolfner 2007). Functional redundancy in seminal fluid components is also observed between species. For example, throughout the animal kingdom, seminal fluid typically contains multiple proteases, protease inhibitors, lipases, lectins, and CRISPs (Ravi Ram and Wolfner 2007). Many of the genes that encode these components evolve extremely rapidly and hence Sfp genes often lack orthologs even among very close relatives (Sirot et al. 2014).

MATERIALS AND METHODS

Sample preparation

We tested the responses of males exposed to each other for 2, 26, or 50 h prior to mating and in the absence of females. We chose this design for investigating the effects of exposure to rivals because, as well as being easier to test in an unconfounded manner, the responses of males to each other prior to mating are also larger, in terms of subsequent effects on mating duration, than those observed when females are also present in the mating arena (Bretman et al. 2009). Fly rearing and all experiments were conducted in a 25°C humidified room, with a 12:12 h light:dark cycle. Flies were maintained in glass vials (75 × 25 mm) containing 7 mL standard sugar–yeast (SY) medium (Bass et al. 2007). Wild-type flies were from a large laboratory population originally collected in the 1970s in Dahomey (Benin) and were from the strain used previously in related studies (Bretman et al. 2009, 2010, 2011a,b, 2013). Larvae were raised at a standard density of 100 per vial, on SY medium supplemented with live yeast liquid. At eclosion, males and females were separated using ice anesthesia. Males were collected and stored 10 per vial and 24 h later assigned randomly as rival or focal males, as in our previous studies. The initial housing of males together for 24 h during the collection period has no effect on the subsequent ability of focal males to adapt to the social environment (Bretman et al. 2013). Focal males were then placed one male per vial for 5 d and randomly assigned to the different treatments (±rival, at 2, 26, or 50 h time periods; n = 50 males per treatment per time point). Rival males were given an identifying wing clip under light CO2 anesthesia, a procedure that has no effect on the adaptive response of the focal males (Bretman et al. 2009). Following the 5 d of being housed alone, one rival male was aspirated into each of the focal male +rival treatments at 9 am (rivals were all introduced at the same time of day to control for circadian influences on gene expression). At the same time, males in the corresponding no rival treatments were subjected to the same movements, in order to standardize handling. At 2, 26, and 50 h after introduction of rival males, 50 focal males from each of the ±rivals treatments were aspirated into individual microcentrifuge tubes and immediately frozen in liquid nitrogen. Samples were then stored at −80°C until RNA extraction. This procedure was replicated three times using independent cohorts, to produce three biological replicates.

RNA extraction

Total RNA (Ambion mirVana miRNA Isolation Kit, Life Technologies, AM1561) was extracted from pools of 40 males per treatment, per time point. Whole flies were first separated into HT and A body parts over dry ice. RNA was assessed for quantity and quality using a NanoDrop 8000 Spectrophotometer (Thermo Fisher Scientific) and by running samples on a 1% agarose gel. All yields were above 100 ng/µL in 100 µL in the RNA storage solution (Life Technologies, AM7000) at absorbance λ = 260. RNA was stored at −80°C until use. Nondirectional, single-end RNA-seq was conducted using the Illumina HiSeq2000 platform with 50-nt read length (Baseclear provider). Three samples were run per lane to yield an expected >40 M reads per sample [2 experimental treatments (±rival) × 2 body parts (HT, A) × 3 time periods (2, 26, and 50 h) × 3 biological replicates = 36 samples].

Bioinformatics analysis

The data set was analyzed using standard quality check (QC) procedures (DeLuca et al. 2012; Wang et al. 2012) and new approaches (I Mohorianu, A Bretman, DT Smith, EK Fowler, T Dalmay, T Chapman, in prep.). We assessed (i) read number and complexity (i.e., ratio of nonredundant to redundant reads), (ii) nucleotide composition and strand bias, (iii) genome matching and annotation matching reads (e.g., mRNAs, t/rRNAs, miRNAs, UTRs, introns, intergenic regions), and (iv) correlation and Jaccard similarity indices on the original versus normalized data. Genome matching was performed using PatMaN (Prüfer et al. 2008). To normalize gene abundances, we used subsampling without replacement (Efron and Tibshirani 1993) to a fixed total of 50 M reads (Supplemental Table S2A) followed by quantile normalization (Bolstad et al. 2003) to remove any remaining fine-scale variation. We incorporated the experimental design into the method used for detecting differentially expressed transcripts. First we evaluated the difference in expression between body parts and then the genes expressed only in the target body part (i.e., DE between HT samples for genes expressed in HT; similarly for A samples). The separation of the HT and A body parts during RNA extraction was not 100% precise and, given the high sensitivity of RNA-seq, a “shadow” of a signature of DE in genes that are much more abundantly expressed (and which normally function) in the “other” body part was sometimes obtained (Supplemental Fig. S7). Our method screened out such “leaky” genes.

Functional analysis of transcripts showing DE

We first compiled the set of individual genes showing DE according to a statistically supported threshold and interrogated this set of genes for common functional categories, before testing for significant enrichment/depletion.

Numbers and functions of DE genes

We used a statistical cut off for DE following an analysis of the distribution of DE for each time point and body part for both replicate-to-replicate and rivals/no rivals comparisons. To converge on an appropriate and experiment-wide threshold, we asked what log2(OFC) would give a P value of 0.05 for each time point and treatment. The results showed that thresholds were lower for the A samples, with the average threshold for P < 0.05 being log2(OFC) = 0.49, 0.45 (for rep–rep and ±rivals comparisons, respectively). Adding in the HT data gave an overall log2(OFC) = 0.64, 0.66 for experiment-wide P < 0.05). On the basis of this we chose log2(OFC) = 0.5 (actual OFC of 1.5 and above) as the minimum statistically supported threshold for the detection to DE. This threshold corresponds to the fold change of ∼1.4 and above that can be independently validated using low-throughput methods such as qRT-PCR (Morey et al. 2006). Using this threshold we compiled the set of DE transcripts showing >1.5 OFC for each body part and time period.

Frequency distributions of gene expression for functional categories

In addition to examining DE in individual genes, we tested for changes in the expression of classes of functionally similar gene groups. Genes showing evidence for DE were interrogated for common functional categories using manually curated GO Biological Processes terms (The Reference Genome Group of the Gene Ontology Consortium 2009) from the Biological Processes GO vocabulary at the second to third level of GO hierarchy (The Reference Genome Group of the Gene Ontology Consortium 2009). This was to prevent differences in resolving power introduced by the same term located at different hierarchical levels. Manual investigation of the keywords was useful to ensure the keywords captured a tight, single functional class for which there was a good level of experimental evidence for functional assignment (equivalent to a focused, standard GO analysis [Supplemental Material 1]). To test for significant enrichment of the functional categories, the selected keywords were used to determine, for all the genes in that category, the frequency density distribution of expression for the ±rivals treatments (rivals—no rivals expression). Categories of genes not responding to rivals should show a frequency distribution centered on zero. In contrast, responding categories of genes should show frequency density distributions with a significantly wider spread than the control, or significantly shifted to the right or left. Control frequency density distributions in each case were a randomly drawn set of the same number of genes of similar length, generated using a random uniform distribution with total abundance (sum of abundance across all samples) >100. Changes in location were best captured using the t-test (expected and observed data both quasinormally distributed) and changes in shape by the Kullback–Leibler (KL) divergence (with P and Q being the frequency distributions for the keyword and random control genes, respectively). A KL value close to zero indicated similar distributions, and positive or negative values dissimilar distributions.

RNA-seq detection limits

We first investigated whether there were limits of detection for the qRT-PCR validation. In pilot experiments, sample aliquots were sent to a provider (qStandard) to validate four genes with a range of high- and low-expression values from the RNA-seq data (Che87a, to, TotA, and PebII). The geNorm Kit of six loci (PrimerDesign) was used to identify reference genes with stable expression across time and tissues. Three genes (αTub84B, eIF-1A, and CG13220), whose expression was consistent across all samples (as evidenced by the gene stability measure M, all <1.5; [Vandesompele et al. 2002]) were chosen as reference genes. There was low agreement between qRT-PCR and RNA-seq for genes with very low normalized abundance (Supplemental Fig. S9). Hence we chose genes of interest (GOIs) for the main tests with abundances >100. GOIs with very high normalized abundances (>12,000) were also avoided to remove the potential for saturation.

Reference genes

For the main qRT-PCR validation tests, we selected three reference genes. We assessed variation in expression for reference genes commonly used for qRT-PCR validations (Ling and Salvaterra 2011) in our RNA-seq expression data and from the Fly Atlas microarray expression data (Chintapalli et al. 2007). From this, a list of candidate reference genes of genes not differentially expressed was produced. For the selected reference genes (α-Tub84b, eif-1a, and CG13220), we next confirmed that the presence plot for each gene matched its corresponding FlyBase model (St Pierre et al. 2014). We then selected primers using FlyPrimerBank (Hu et al. 2013) and relevant publications (Ling and Salvaterra 2011) and optimized them for qRT-PCR.

Genes of interest (GOI)

We chose 20 GOIs for validation based on (i) transcript abundance, (ii) the degree of DE combined with functional annotations that intersected with the main biological responses of interest, (iii) clarity and concordance of the gene model from FlyBase with the RNA-seq presence plots, and (iv) the ability to design unique primers (Supplemental Table S9). The aim was to assess whether low-throughput validation of the same biological samples gave the same patterns of gene expression and hence, to validate the patterns being described by these data.

qRT-PCR

For primer optimization we used HT or A RNA from wild-type flies. We removed DNA using TURBO DNase (Life Technologies) on 2 µg total RNA. DNase was deactivated using the resin in the TURBO DNA-free Kit (Life Technologies), and all samples were tested for DNA contamination with no-reverse transcription controls using the three reference genes as targets. The QuantiTect Reverse Transcription Kit (Qiagen) was used to reverse transcribe 1 µg total RNA to cDNA, which was then stored at −20°C. qRT-PCRs were run using a StepOnePlus machine (Life Technologies) using iTaq Universal SYBR Green Supermix with triplicate technical replicates, using 2 ng RNA as template in 20 µL reactions in MicroAmp Fast Optical 96-Well Reaction Plate with Barcode, 0.1 mL (Life Technologies). The qRT-PCR conditions were 95°C for 30 sec, followed by 40 cycles of 95°C for 15 sec and 60°C for 1 min and data acquisition. Following the qRT-PCR, we ran a melt curve analysis (on default settings). All samples showed a single peak on the melt curve. We included a no template control for each GOI on each plate. All primers were used at a final concentration of 500 nM. For primer optimization (Eurofins MWG Operon), we produced standard curves using at least five 1:5 dilutions of RNA starting at 50 ng cDNA. Primer sequences were taken from FlyPrimerBank and citations (Ling and Salvaterra 2011; Hu et al. 2013). Gene names, primer sequences, amplicon length, efficiency, R of primers are given in Supplemental Table S9. We tested for concordance between the qRT-PCR and RNA-seq data by counting the number of FC “agreements,” i.e., the number of times that genes showing DE in the RNA-seq were also found in the qRT-PCR data, and vice versa. GOI were validated if this measure was ≥0.6.

DATA DEPOSITION

The data described in this study are publicly available on the Gene Expression Omnibus (GEO) (Barrett et al. 2013) under accession number GSE55839 (GSM1346870-GSM1347005). All Perl and R scripts will be made available on GitHub.

SUPPLEMENTAL MATERIAL

Supplemental material is available for this article.