Mating Changes Sexually Dimorphic Gene Expression in the Seed Beetle Callosobruchus maculatus.

Sexually dimorphic phenotypes arise largely from sex-specific gene expression, which has mainly been characterized in sexually naïve adults. However, we expect sexual dimorphism in transcription to be dynamic and dependent on factors such as reproductive status. Mating induces many behavioral and physiological changes distinct to each sex and is therefore expected to activate regulatory changes in many sex-biased genes. Here, we first characterized sexual dimorphism in gene expression in Callosobruchus maculatus seed beetles. We then examined how females and males respond to mating and how it affects sex-biased expression, both in sex-limited (abdomen) and sex-shared (head and thorax) tissues. Mating responses were largely sex-specific and, as expected, females showed more genes responding compared with males (∼2,000 vs. ∼300 genes in the abdomen, ∼500 vs. ∼400 in the head and thorax, respectively). Of the sex-biased genes present in virgins, 16% (1,041 genes) in the abdomen and 17% (243 genes) in the head and thorax altered their relative expression between the sexes as a result of mating. Sex-bias status changed in 2% of the genes in the abdomen and 4% in the head and thorax following mating. Mating responses involved de-feminization of females and, to a lesser extent, de-masculinization of males relative to their virgin state: mating decreased rather than increased dimorphic expression of sex-biased genes. The fact that regulatory changes of both types of sex-biased genes occurred in both sexes suggests that male- and female-specific selection is not restricted to male- and female-biased genes, respectively, as is sometimes assumed.

Introduction

Differences in morphology, physiology and behavior between the sexes are nearly ubiquitous in sexually reproducing organisms and offer some of the most striking examples of intraspecific variation in the natural world. Evolution of sexual dimorphism requires both distinct selection in the sexes and a sex-specific genetic architecture. Because the sexes largely share the genome, the majority of sexually dimorphic phenotypic variation arises and persists via differential regulation of gene expression in females and males (Ellegren and Parsch 2007; Williams and Carroll 2009; Parsch and Ellegren 2013). This has now been demonstrated in several taxa (Arbeitman et al. 2002; Hahn and Lanzaro 2005; Yang et al. 2006a; Lebo et al. 2009; Prince et al. 2010; Naurin et al. 2011). Sex determination hierarchy (Williams and Carroll 2009), sex-linkage (Meisel et al. 2012) and miRNA’s (Fagegaltier et al. 2014) participate in the generation and maintenance of sex-biased expression. Genes expressed at a higher level in males (i.e., male-biased genes) are thought to evolve primarily under male-specific selection and those with higher expression in females (i.e., female-biased genes) primarily under female-specific selection (Ellegren and Parsch 2007). In line with this, both comparative genomic (Harrison et al. 2015) and experimental evolution (Hollis et al. 2014; Immonen et al. 2014) studies have shown that sex-biased genes evolve regulatory differences in response to variation in the mating system. Studies across taxa have also demonstrated that male-biased genes generally tend to evolve faster at both regulatory and sequence levels than female- and un-biased genes (Zhang et al. 2004, 2007; Grath et al. 2009; Harrison et al. 2015), presumably due to stronger sexual selection in males.

In addition to the sex of the organism, changes to the internal (De Gregorio et al. 2001; Zhou et al. 2014) and external abiotic environments (Smith et al. 2013; Hoang et al. 2015) as well as social conditions (Carney 2007; Cummings et al. 2008; Immonen and Ritchie 2012) are major causes of variation in the gene expression patterns. A static genome can track environmental changes by altering the regulation of gene expression, and understanding the transcriptomic basis of phenotypic plasticity can thus provide important insights into how a single genotype can generate different phenotypes (Zhou et al. 2012; Kijimoto et al. 2014; Perry and Mank 2014). Thus far, studies have explored to what extent sexual dimorphism changes across ontogeny (Arbeitman et al. 2002; Perry et al. 2014) and under different developmental conditions (Wyman et al. 2010; Ledon-Rettig and Moczek 2016). The likely possibility that the pattern of sex-biased expression changes during adult life, due to changes in reproductive state, has to our knowledge not been previously investigated.

Here, we test for plasticity in gene expression between the sexes in response to mating in Callosobruchus maculatus seed beetles, a model species in sexual selection and conflict studies (Zuk et al. 2014). Mating triggers major physiological and behavioral changes that differ between the sexes, and is therefore likely to affect the expression of sex-biased genes. Molecular social interactions that occur during and after mating are predicted to evolve under sex-specific selection due to sperm competition, cryptic female choice and interlocus sexual conflict (Rice and Gavrilets 2014). For example, males transfer a cocktail of seminal fluid proteins and peptides (Sfps) to females, many of which exercise physiological effects in females (Wolfner 2009). In C. maculatus, 98 Sfps have thus far been identified (Bayram et al. 2016). Many of these proteins differ in their abundance across populations and are associated with variation in male sperm competitiveness and ability to influence female fecundity and lifespan (Goenaga et al. 2015). Studies in D. melanogaster have shown that females undergo dramatic changes not only in egg production following mating (Chapman et al. 2001), but also in sleeping and locomotion (Isaac et al. 2010), eating (Carvalho et al. 2006), re-mating (LaFlamme et al. 2012), immunity (Peng et al. 2005; Innocenti and Morrow 2009) and ageing (Chapman et al. 1995). Whilst such a comprehensive analysis of female mating responses in seed beetles is yet to be undertaken, SFPs injected into females alter egg laying rate, female remating rate and male competitive fertilization success (Yamane et al. 2015). These phenotypes likely result from postmating changes in expression of hundreds to thousands of genes in females in both reproductive and somatic tissues (Lawniczak and Begun 2004; McGraw et al. 2004, 2008; Mack et al. 2006; Innocenti and Morrow 2009; Dalton et al. 2010; Gioti et al. 2012).

Compared with females, male postmating physiological and molecular changes are poorly understood. In male C. maculatus, high mating frequency results in a rapid decline in the ejaculate volume (Rönn et al. 2008), which should require a subsequent renewal of the ejaculate components (i.e., lipids, proteins and polysaccharides). In line with this, males show a postmating increase in metabolic rate (Immonen et al. 2016b). A handful of genes are known to change expression following mating in male D. melanogaster (Ellis and Carney 2010) and Ceratitis capitata (Gomulski et al. 2012). Studying the transcriptomic basis of mating responses in each sex allows identifying candidate genes and genetic pathways underlying reproductive phenotypes. Thus far, however, our knowledge of such mating responses is largely limited to females of D. melanogaster (McGraw et al. 2004, 2008; Mack et al. 2006; Innocenti and Morrow 2009; Dalton et al. 2010; Zhou et al. 2014), Anopheles mosquitoes (Kocher et al. 2008; Rogers et al. 2008; Alfonso-Parra et al. 2016) and Ceratitis capitata fruit flies (Gomulski et al. 2012).

In the present study, we tested how gene expression dimorphism changes as a result of mating, using RNA-seq. First, we characterized sexual dimorphism in transcript abundance in virgin beetles, comparing both reproductive and somatic tissues. We also explored whether any genes show sex differences in isoform-specific expression due to alternative splicing. This is the first study to characterize sex-biased expression in C. maculatus. Second, we investigated how mating changes gene expression, specifically testing how each sex may alter the patterns of sexually dimorphic expression. Under the premise that female reproduction involves primarily female-biased genes and male reproduction male-biased genes, we predicted that females more commonly upregulate female- and males upregulate male-biased genes in response to mating. This would subsequently result in an increase in sexual dimorphism relative to the virgin state.

Materials and Methods

Sample Preparation

Callosobruchus maculatus is a facultatively aphagous pest of legume plants, occurring in the subtropical and tropical regions throughout the world. It has a polygynandrous mating system. Females lay eggs on seeds and the larval development and pupation are completed within the seeds. Adults are ready to mate as soon as they emerge from the seed. A near isogenic line (from “South-India” stock population, generated by 5 generations of inbreeding) was used in this experiment. The beetles were reared on mung beans (Vigna radiata) in laboratory climate cabinets at 29°C, 60% RH and a 12-L:12-D light cycle. The parental generation laid eggs over a 24-h period in noncrowded conditions with a surplus of host beans. Beans were isolated prior to adult hatching and virgin males and females were subsequently collected and stored in isolation. In the mating treatment group, 1-day-old males and females were briefly paired and were separated after mating was terminated. Same-aged mated and virgin beetles were kept individually for 24 h with beans, after which they were snap-frozen with liquid nitrogen and stored in −80°C. This time point was chosen as it corresponds to the time when females lay most of their eggs (Fox 1993; Immonen et al. 2016a) and males still undergo increased postmating metabolism, presumably due to ejaculate renewal (Immonen et al. 2016b). We separated the head and thorax (the somatic tissues) from the abdomen (the reproductive tissues) under ice, pooled tissues from six individuals per sample and extracted total RNA using RNeasy Mini Kit (Qiagen), following the manufacturer’s protocol. DNase digestion was applied using DNase I (RNase-Free DNase set by Qiagen). The RNA quality and quantity was assessed and affirmed using NanoDrop, Qubit and Bioanalyzer. Three replicate samples were prepared per sex, treatment and tissue type, resulting in a total of 24 samples.

cDNA Library Generation and Illumina Sequencing

The RNA-seq libraries were prepared using the Illumina TruSeq stranded mRNA sample preparation kit according to the manufacturer’s guidelines (Illumina 2013). Poly-A containing RNA was purified from total RNA using poly-T oligo attached magnetic beads, after which mRNA was fragmented and reverse transcribed to first strand cDNA using random primers. The cDNA fragments were ligated to adapters and purified cDNA libraries enriched with PCR. All sequencing was performed using Illumina HiSeq 2500 sequencing technology producing 100-bp length paired-end reads. All libraries were sequenced in two lanes. The two technical replicates of each sample were subsequently pooled in the assembly to maximize the coverage. The library generation and sequencing were performed by SNP&SEQ Technology Platform at Uppsala University.

Transcriptome Assembly and Annotation

Due to the lack of closely related reference genome, the transcriptome was de novo assembled, using all libraries as well as additional samples from different developmental stages. All samples were deeply sequenced to cover high and low expressed genes. In total, more than 492 million read pairs were sequenced and used to generate a reference transcriptome with Trinity assembler (Grabherr et al. 2011; Haas et al. 2013). The de novo assembly has been described in detail in Sayadi . We separately mapped back the reads from each sample to the assembled transcriptome in order to obtain the read counts as the measure of expression using RSEM method. Trinity transcript (gene isoform) clusters with the longest isoform were selected as representatives of genes in the subsequent gene-level analysis (hereafter referred to as genes for convenience). The assembled transcriptome generated 218,192 transcripts which correspond to 145,883 putative genes. De novo transcriptome assemblies generate far more transcripts than what is a biologically realistic number of genes for several reasons (Sayadi ). We thus considered only the predicted genes that contained an ORF (a total of 29,812 genes). Further, we required a gene to be expressed by at least 2 counts per million reads (cpm) in at least three libraries. This reduced the dataset into 12,412 expressed genes in the abdomen libraries and 11,143 in the head and thorax libraries for the gene-level analysis for sex-biased expression. The female libraries used for the mating response analyses contained 12,333 expressed genes and the male libraries 12,855 expressed genes.

Statistical Analyses

The analyses were performed using edgeR v.3.10.5 (Robinson et al. 2010) and limma v.3.24.15 (Smyth 2005) packages within Bioconductor (Gentleman et al. 2004), using R v.3.2.2 (RDevelopmentCoreTeam 2011). We followed the pipeline offered in edgeR (Robinson et al. 2010): after normalization (by computing scaling factors using the trimmed mean of M values, or TMM, and adjusting the library sizes accordingly to account for putative differences in the RNA composition between libraries) a generalized linear model with likelihood ratio tests was used, where a negative binomial model is fitted with a Cox-Reid profile-adjusted likelihood ratio method to estimate tag-wise dispersion. Putative differential exon usage between the sexes (i.e., alternative splicing) was tested by fitting a negative binomial generalized log-linear model at the exon level and subsequently comparing the log fold change of each exon to that of the entire gene, as implemented in edgeR. In all analyses, we used a statistical significance threshold of 5% False Discovery Rate (FDR) unless stated otherwise (Benjamini and Hochberg 1995). Overrepresentation of Gene Ontology terms (Biological Process and Molecular Function) was tested with the GOstats package v.2.34.0, using a conditional hypergeometric test with a P value cutoff < 0.05 (Falcon and Gentleman 2007). The gene universe was defined as all the expressed genes in a given condition.

Our inferential tests include: (1) a test of differential expression between virgin males and females within each tissue type, to identify sex-biased genes. To call sex-bias, we required at least a 2-fold difference in expression between the sexes (i.e., log2FC ≥ 1 or −1) (Montgomery and Mank 2016), in addition to the statistical significance threshold. Considering only genes with large expression difference between the sexes limits the cases where sex-bias in expression may result from differences in the composition of the shared tissues rather than from regulatory differences per se (Montgomery and Mank 2016). Expression differences due to different reproductive organs naturally remain. We present the data on expression differences between females and males with different fold-change cut-offs in the supplementary figure S1, Supplementary Material online. (2) We examined patterns of sex-biased expression in genes with a single versus multiple expressed transcripts to test whether any genes show sex differences in isoform-specific expression in virgins, due to alternative splicing, separately for the head and thorax and abdomen. (3) To achieve our main goal of identifying mating induced changes in sexually dimorphic expression, we fitted a two-way ANOVA-type model (in edgeR) with a specific contrast to test for the interaction effect between sex and mating status, separately for each tissue class. Here, interactions signify effects of mating status on sexual dimorphism. We then compared gene expression changes in response to mating within each sex (for each tissue class). These results were combined to assess how each sex contributes to changes in sex-bias (i.e., a switch from male- to female-bias, or vice versa, or to un-biased expression). For identifying genes that change sexually dimorphic expression due to mating, we required that they show a significant mating-by-sex interaction and respond significantly to mating in at least one of the sexes. The additional criterion of the main effect of mating allows investigating which sex is causing the putative sex-bias plasticity with high statistical confidence. Moreover, for classifying those plastic genes that switch to opposite sex-bias due to mating we also required at least a 2-fold expression difference (with 5% FDR) between the sexes in both virgin and mated states, for the reasons stated earlier. Proportion tests were used for assessing whether sex-biased genes were overrepresented among the mating responsive genes in either sex, as well as to test for deviations in the up/downregulation patterns of sex-biased mating responsive genes from the expected (the expected patterns were calculated as the proportion of all up/downregulated genes).

We refer to a pattern of expression where mated females change to more resemble males as de-feminization (i.e., when female-biased genes decrease and/or male-biased genes increase expression relative to virgin females). We refer to de-masculinization of males as a pattern where mated males increase the expression of female-biased and/or decrease male-biased genes relative to their virgin state.

Results

Sexual Dimorphism of the C. maculatus Virgin Transcriptome

First, we analysed male and female transcriptomes separately in the abdomen and head/thorax in virgin beetles, in order to characterize how many and what types of genes differ in expression between the sexes in somatic and reproductive tissues of C. maculatus. We identified a total of 6,645 differentially expressed genes in the abdomen (54% of the 12,412 genes expressed in the abdomen, supplementary fig. S1a, Supplementary Material online) and 1,398 genes in the head and thorax (13% of the 11,143 expressed genes, supplementary fig. S1b, Supplementary Material online), with at least 2-fold expression difference between virgin females and males (FDR 5%). Even when using a highly stringent statistical criteria of 0.1% FDR (i.e., Padj≥ 0.001), 51% (6,343 genes) and 9% (1,018 genes) of the transcriptome was sex-biased in the abdomen and head/thorax tissues, respectively. In the abdomen, there were more male- than female-biased genes (≥ 2-fold change: 54% vs. 46% of the sex-biased genes, respectively). This was due to a higher number of highly (≥ 10-fold change) male-biased genes (10%) in the transcriptome compared with female-biased (4%) (see supplementary fig. S1a, Supplementary Material online). The opposite was true in the head and thorax, where 65% of the sex-biased genes showed higher expression in the females, and this was the case even for the highly sex-biased genes (≥10-fold change: 2% female- vs. 0.5% male-biased, supplementary fig. S1b, Supplementary Material online). There was a large overlap of sex-biased genes between the reproductive and somatic tissues (1,040 genes, fig. 1). This pattern was, however, greater for the female-biased genes, of which 23% were shared between the tissues in contrast to 7% of the male-biased genes (supplementary fig. S2, Supplementary Material online). The overlap of both types of sex-biased genes between the tissues was mainly due to the head and thorax: 74% of the sex-biased genes expressed in the head and thorax were also sex-biased in the abdomen, while only 16% of the abdomen sex-biased genes were shared between the tissue classes (fig. 1). The 50 most sex-biased genes of both tissue classes are presented in tables 1–4.

F

(A) A Venn diagram of the numbers of sex-biased genes in virgins (i.e., Padj <0.05 and at least 2-fold sex-bias) in the abdomen and the head + thorax. (B) Numbers of genes showing plasticity in sex-biased expression due to mating (i.e., sex-by-mating interaction and mating response at least in one sex with Padj < 0.05) in each tissue class. Overlap of the genes that respond to mating in each sex and show sex-bias plasticity in (C) the abdomen and (D) the head and thorax tissues.

Female-biased genes of the abdomen were enriched with biological processes related to metabolism, cellular response to stress, response to hormones and regulation of transcription, for example, supplementary table S1a, Supplementary Material online. Similar terms were also enriched among the female-biased genes of the head and thorax (supplementary table S1c, Supplementary Material online). Male-biased genes of the abdomen showed overrepresentation of biological processes related to oxidation–reduction, carbohydrate metabolism and G-protein coupled receptor signaling pathway, for example, supplementary table S2a, Supplementary Material online. Male-biased genes in the head and thorax were enriched, for example, with genes involved in visual perception, detection of chemical stimulus and neurotransmitter transport (supplementary table S2c, Supplementary Material online). We note that the GO analyses are a useful exploratory tool, but should be interpreted as guidance rather than as conclusive evidence of functional differences (Pavlidis et al. 2012). The lack of closely related and well-annotated genomes can mean that some of the potentially most interesting genes under sex-specific selection may lack annotation due to faster divergence. In addition to the Gene Ontology annotations, we also explored whether the recently identified Sfp encoding genes of C. maculatus (Bayram et al. 2016) were overrepresented among the male-biased genes of the abdomen. These 98 genes were more common than expected among the male-biased genes when considering all the genes with the Padj < 0.05 (1.0% vs. 0.78%; X2 = 111, P = 0.0078, supplementary fig. S3, Supplementary Material online), but they were not significantly overrepresented among those with >2-fold sex difference. Although 15 of these 98 genes were indeed strongly male-biased, most were only moderately male-biased in expression and some even showed female-bias (supplementary fig. S3, Supplementary Material online).

In addition to sex differences in the amount of expression, males and females can also differ in the patterns of gene splicing. In order to test for such putative differences in alternative splicing, we examined changes in the isoform-specific expression differences between males and females. We considered transcripts that show at least a 2-fold expression difference as sex-biased, with 5% FDR. In the abdomen, 20,683 transcripts were retained in the analysis after pre-filtering (transcripts with > 2 cpm in at least three libraries were retained), corresponding to 11,839 genes. The majority of the genes had only a single expressed transcript in our dataset (7,258 genes) and 54% of these were sex-biased (3,903 genes). Among the genes with multiple expressed transcripts (4,581), the average number of isoforms per gene was 2 with a maximum of 17. 68% of these genes (3,128) had at least one sex-biased isoform. This is a significantly greater proportion than for the single-transcript genes (X2 = 245.05, df = 1, P value < 2.2e-16). We also compared the numbers of sex-biased transcripts from multi- versus single-transcript genes: 61% (6,080) of the sex-biased transcripts were from the former, while 3,903 transcripts (39%) were from the single-transcript genes (X2 = 949.48, df = 1, P value < 2.2e-16). These patterns are consistent with the notion that the evolution of sex-biased expression is less constrained in multi-exon compared with single-exon genes. This interesting observation should, however, be verified with a well-annotated genome. We found that 28% of all the multi-transcript genes (1,371, 12% of the whole transcriptome) showed significant sex differences in the isoform-specific expression in the abdomen, in line with alternative splicing. However, the proportion of these genes was lower compared with the multi-transcript genes that showed a uniform pattern of sex-bias across the isoforms (1,837 genes, 40%; X2 = 144, df = 1, P value < 2.2e-16). The putative alternatively spliced genes of the abdomen were enriched with biological processes related to metabolic processes and oxidation–reduction (supplementary table S3a, Supplementary Material online). Among the Sfp genes of C. maculatus, 27 out of 98 showed signs of alternative splicing. For 14 genes, the transcripts showed a mixture of male-biased and un-biased expression (supplementary table S4a, Supplementary Material online), two showed different degrees of male-bias (supplementary table S4a, Supplementary Material online), two a mixture of male-, female- and un-biased expression (supplementary table S4b, Supplementary Material online), two a mixture of male- and female-bias (supplementary table S4b, Supplementary Material online) while the remaining 8 showed a mixture of female- and un-biased expression (supplementary table S4c, Supplementary Material online).

In the head and thorax, 18,278 expressed transcripts (from 10,551 genes) were retained after filtering (transcripts with > 2 cpm in at least three libraries were retained), of which the majority had only one expressed transcript per gene (6,586). Overall, 3,965 genes had more than one expressed transcript (on an average there were 2 expressed transcripts per gene, with a maximum of 16). Similar to the abdomen, multi-transcript genes of the head and thorax showed sex-biased expression (in at least one of the transcripts) more often than the single-transcript genes (16% vs. 11%, respectively: X2= 61.914, df = 1, P value = 3.588e-15). A total of 316 (6.6% of the multi-transcript genes and 3% of the whole transcriptome) genes showed sex differences in isoform-specific expression in the head and thorax. Similar to the abdomen, this is a lower proportion than of those multi-transcript genes that showed a uniform sex-bias across the isoforms (390 genes, 10%: X2= 27.8, df = 1, P value = 1.312e-07). The alternatively spliced genes in the head and thorax were enriched with biological processes related to carbohydrate metabolism (supplementary table S3c, Supplementary Material online). The alternatively spliced genes were largely distinct in each tissue class (supplementary fig. S4, Supplementary Material online).

Plasticity in Gene Expression Due to Mating

Our main goal was to study whether and how the sex-biased genes identified in virgins respond to mating in each sex. We did this by first testing for significant two-way interactions between mating treatment and sex: 1,305 genes showed a significant interaction in the abdomen and 309 in the head and thorax. Of these, 1,041 genes (15.7% of the virgin sex-biased genes) in the abdomen and 243 (17.4%) in the head and thorax also showed a significant main effect of mating in at least one of the sexes and were thus considered further as genes showing plasticity in the degree of sexual dimorphism due to mating. The abdomen and the head and thorax showed largely different responses to mating: only 59 genes showed a similar interaction effect in both tissue groups (fig. 1), despite the large shared proportion of the sex-biased genes between them (fig. 1). Females and males both contributed to the sex-bias plasticity but largely of different genes (fig. 1).

Next, we estimated how many of the identified “plastic” sex-biased genes changed status from the virgin state. The majority of the genes showed subtle changes in the degree of sex-bias, but in 2% of the genes in the abdomen and 4% in the head and thorax the sex-bias status were altered (fig. 2table 5). Most of such cases changed from showing sex-biased expression to showing no bias (79 genes in the abdomen and 48 genes in the head and thorax, table 5). However we also identified 9 and 6 genes that changed from female- to male-biased expression, while 8 and 5 genes changed from male- to female-biased expression in the abdomen and head + thorax, respectively (i.e., showed significant sex-bias by 2-fold as both virgin and mated; table 5 and fig. 3). Only a minority of these genes (5 out of 28) switched sex-bias due a simultaneous but opposite expression change in both sexes. In 6 genes, the switch to male-bias occurred due to upregulation by mated males whereas in 7 genes there was a strong downregulation in the mated females without any change in the male expression. Switches to female-bias occurred mainly due to increased expression only in females (6 genes) but 2 genes switched to female-bias by strong downregulation only in males (fig. 3). See supplementary table S5, Supplementary Material online, for the data and annotations for these genes that switch sex-bias.

F

Sex-biased expression in virgin and mated beetles for the plastic genes in (A) the abdomen (1,041 genes) and (B) the head and thorax (243 genes). Grey dashed lines indicate a 2-fold difference in expression between the sexes (i.e., logFC = 1 or − 1. Positive values indicate female- and negative male-biased expression). Off diagonal segments show the genes that have changed sex-bias status due to mating (upper left segment: from male- to female-bias, lower right segment: from female- to male-bias).

Table 5

Numbers of Genes Showing a Significant Switch in Sex-Bias Status Due To Mating

Tissue	F-b to M-b	F-b to Nb	Nb to F-b	M-b to F-b	M-b to Nb	Nb to M-b
Abdomen	9	49	10	8	30	19
H+T	6	13	11	5	35	10

Note.—The columns denote changes in sex-bias from virgin to mated state.

F-b, female-bias; M-b, male-bias; Nb, no sex-bias.

F

Expression response to mating in males and females for the genes that switch sex-biased expression in the two tissue groups (i.e., sex-biased ≥ 2-fold in both virgin and mated beetles, with FDR 5%). Two left-most columns show the degree of sex-bias in the virgin and mated beetles (positive log2FC indicates female-biased expression) and the two right-most columns the degree of mating response in males and females for the same genes (positive log2FC indicates higher expression in the mated beetles). Absence of a record (white) indicates no significant difference between mated and virgin individuals.

To further test how each sex responds to mating and how this might alter sexually dimorphic expression, we focused on examining the transcriptomic changes of sex-biased genes (as identified in virgins) separately within each sex in each tissue class. Nearly 6 times more genes responded to mating in the female abdomen compared with males: 2,286 genes significantly changed expression in the females and 382 genes in the males (Padj < 0.05). Female mating response was enriched with virgin sex-biased genes: 59% (1,351) were sex-biased, a significantly greater proportion than in the whole abdomen transcriptome (X2df= 24.11, P < 0.0001). Female-biased genes were more common than expected relative to their frequency among the virgin sex-biased genes (771 genes, 57% of the sex-biased mating responsive genes: X2df= 61.41, P < 0.0001). However, opposite to our prediction, a great majority of these female-biased genes were downregulated in mated relative to virgin females (707 genes, 92%: X2df= 452.21, P < 0.0001), whereas male-biased genes responded mainly by upregulation (541 genes, 93%: X2df= 509.11, P < 0.0001) (fig. 4), resulting in de-feminization of the female abdomen relative to the virgin state. The genes that increased expression due to mating were enriched for example with biological processes involved in oxidation–reduction, carbohydrate metabolism, lipid metabolism and proteolysis (supplementary table S6a, Supplementary Material online), while the downregulated genes showed for example overrepresentation of protein modification process, regulation of signal transduction and DNA repair (supplementary table S6b, Supplementary Material online).

F

Mating response in the female (A and B) and male (C and D) abdomen and the head + thorax, respectively, as a function of average expression (logCPM). Sex-biased genes in virgins (Padj < 0.05 and at least 2-fold sex-bias) are highlighted. Positive values indicate a higher expression in the mated state relative to virgin.

In contrast to female abdomen, in the male abdomen the virgin sex-biased genes were significantly underrepresented among the mating responsive genes (118 genes, 31%: X2df = 76.3, P < 0.0001; fig. 4). Male-biased genes were also not more common than expected (68 genes, 57%: X2df = 0.785, P = 0.3756). However, those male-biased genes that did respond were more often downregulated than expected (41 genes, 60%: X2df= 25.71, P < 0.0001), which is again opposite to our prediction. Female-biased genes were, however, not more commonly upregulated than expected in mated males’ abdomen (31 genes, 62%: X2df= 0.921, P = 0.33). Functional enrichment analysis showed that males also increased expression of oxidation–reduction genes in the abdomen (supplementary table S7a, Supplementary Material online), for example. We note that these were different redox genes than those that responded in females (results not shown).

In the head and thorax, 558 genes responded to mating in females and 421 in males (Padj < 0.05). Sex-biased genes were more common than expected among the female mating responsive genes (87 genes, 16%: X2df = 4.41, P = 0.0350). Although female-biased genes were not overrepresented relative to their expected proportion (61 genes, 70%: X2df = 0.99 1, P = 0.3192), they were much more commonly downregulated due to mating (48 genes, 79%: X2df= 25.0 1, P < 0.0001). This is a similar pattern to the abdomen female-biased genes. Male-biased genes did not show a deviant pattern from expected (14 genes, 54%: X2df= 0.0041, P = 0.9495) (fig. 4). Females’ upregulated, for example, signal transduction genes and, in contrast to the abdomen, downregulated metabolic and oxidation–reduction genes in the head and thorax (supplementary table S7c and d, Supplementary Material online).

Similar to the male abdomen, sex-biased genes were not particularly common among the mating responsive genes in the male head and thorax (50 genes, 12%: X2df= 0.17 1, P = 0.6837), and male-biased genes were also not overrepresented (16 genes, 32%: X2df = 0.089 1, P = 0.7653). However, the male-biased genes that did change expression responded exclusively by downregulation (16 genes, 100% X2df= 8.5 1, P = 0.00356) and a great majority of female-biased genes by upregulation (32 genes, 94%: X2df= 53.0 1, P < 0.0001) (fig. 4), suggesting a de-masculinization of the male somatic tissues due to mating. See supplementary table S7b and c, Supplementary Material online, for the enriched Gene Ontology terms.

The full result outputs for the sex-biased genes of virgins, genes that show plasticity in sex-bias, as well as the mating responsive genes of females and males are presented in supplementary table S8, Supplementary Material online.

Discussion

Characterizing sex-specific patterns of plasticity in adult gene expression represents a gap in our understanding the functional operation and evolutionary dynamics of sex-biased genes. Here, we addressed this by investigating how stable sex-bias is across different reproductive states and how sex-biased genes are regulated within each sex. We examined sex-biased expression in reproductive and somatic tissues of C. maculatus in order to investigate how genes with sex-biased expression respond to changes in reproductive status. Mating clearly induces major sex-specific behavioral and physiological changes. We predicted that females would primarily increase the expression of female- and males of male-biased genes, under the premise that the degree of relative expression between the sexes reflects the functional importance of the focal gene product for that sex. This is generally assumed in molecular evolutionary studies which predict that male-biased genes evolve under male-specific and female-biased genes under female-specific selection (Ellegren and Parsch 2007). Our findings yield a number of novel insights. We show that mating affects the expression of hundreds of genes in a largely sex-specific manner (fig. 1). This results in a change in the magnitude of sexually dimorphic expression in both reproductive and somatic tissues (fig. 2). However, the status of sex-bias does not change due to mating for the majority of the virgin sex-biased genes. Both sexes respond to mating by changing the expression of both male- and female-biased genes but do so in largely opposite directions (fig. 4). Contrary to our prediction, mating decreases dimorphic expression of sex-biased genes rather than increases it. Below, we will present the characteristics of sexually dimorphic gene expression in C. maculatus and discuss in detail how each sex responds mating.

Sexually Dimorphic Transcriptome of C. maculatus

Advances in high transcriptome profiling have afforded new insight into the widespread occurrence on sex-biased gene expression in diverse species including brown alga (Martins et al. 2013), nematodes (Albritton et al. 2014), Drosophila (Parisi et al. 2004; Zhang et al. 2004, 2007; McIntyre et al. 2006; Baker et al. 2007; Assis et al. 2012; Allen et al. 2013), mosquitoes (Hahn and Lanzaro 2005; Baker et al. 2011), birds (Mank et al. 2007), zebrafish (Small et al. 2009), and mice (Yang et al. 2006a). However, despite the fact that beetles are uniquely diverse, and include many pest species, genomic tools in this large insect order are still poor. Characterization of sexual dimorphism in the transcriptome has thus far only been conducted in the red flour beetle Tribolium castaneoum (Prince et al. 2010) and the horned beetle Onthophagus taurus (Kijimoto et al. 2014). In line with previous studies, we found extensive sex differences in gene expression: in the abdomen, 54% of the expressed genes were sex-biased and in the head and thorax 13% showed sex-biased expression (supplementary fig. S1, Supplementary Material online). Extremely male-biased genes (by at least 10-fold) were more common in the abdomen than extremely female-biased genes. This is in accordance with many other studies that find more genes with stronger male-biased expression, particularly in the reproductive tissues (reviewed in Ellegren and Parsch 2007; Zhang et al. 2007; Small et al. 2009; Martins et al. 2013). In contrast, we find that in the somatic tissues of C. maculatus female-biased genes were nearly twice as numerous as male-biased genes (supplementary figs. S1, Supplementary Material online). This could reflect the generally more pleiotropic nature of female-biased genes, which in Drosophila show a wider expression breadth than male-biased genes (Assis et al. 2012; Meisel et al. 2012). In accordance, we found that the overlap between the female-biased genes of the abdomen and the head and thorax was far greater than for the male-biased genes (supplementary fig. S2, Supplementary Material online), suggesting that female-biased genes are less specific to sex-limited tissues also in C. maculatus.

Functional analysis revealed sex-specific variation in metabolism of carbohydrates and nitrogen compounds, and in oxidation–reduction processes. One of the most female-biased genes include the gene Cathepsin B-like cysteine protease (CmCatB) of C. maculatus (logFC = 10.4) (table 1). This gene is expressed in the larval gut where it is involved in the hydrolysis of proteins and counter-defence of plant dietary toxins (Koo et al. 2008). It is nearly 2,000 times higher expression in the abdomen of both virgin and mated females suggests that females have specific dietary requirements compared with males. We studied the transcriptome of C. maculatus under aphagous conditions, and therefore CmCatB is not involved in the digestion of food but rather in the decomposition of stored macromolecules, perhaps in the fat body where it is expressed in the moth Helicoverpa armigera (Yang et al. 2006b). Given the importance of diet and metabolism in sex-specific life-history evolution, the metabolic genes with dimorphic expression make interesting, yet under-appreciated, candidates under sex-specific selection for future studies.

Table 1

The 50 Most Sex-Biased Genes with Annotations. Female-Biased (log2FC >1) Genes in the Abdomen

Gene ID	logFC	FDR	GO	Type	Term	BLAST hit
TR72008|c0_g2	13.1	1E-117
TR70149|c0_g1	12.5	2E-243
TR59852|c0_g2	12.1	0E-00	GO:0005634	CC	Nucleus	Uncharacterized protein LOC663065 isoform X2, Tribolium castaneum
			GO:0016567	BP	Protein ubiquitination
			GO:0004842	MF	Ubiquitin-protein transferase activity
			GO:0016874	MF	Ligase activity
			GO:0046872	MF	Metal ion binding
			GO:0008270	MF	Zinc ion binding
			GO:0006511	BP	Ubiquitin-dependent protein catabolic process
			GO:0007275	BP	Multicellular organismal development
TR19738|c0_g2	11.7	0E-00
TR18182|c9_g3	11.7	9E-207
TR11356|c0_g2	11.6	7E-138
TR74056|c0_g1	11.5	2E-171
TR64757|c0_g1	11.4	1E-42
TR29702|c2_g6	11.3	0E-00
TR44608|c4_g2	11.2	5E-258
TR25775|c0_g2	11.2	2E-36
TR38999|c0_g1	11.1	0E-00	GO:0006468	BP	Protein phosphorylation	Hypothetical protein YQE_12596, partial, Dendroctonus ponderosae
			GO:0004672	MF	Protein kinase activity
			GO:0005524	MF	ATP binding
TR41216|c5_g4	11.0	9E-128
TR44677|c1_g2	11.0	6E-49
TR55432|c1_g2	11.0	0E-00
TR32008|c0_g1	11.0	0E-00	GO:0006685	BP	Sphingomyelin catabolic process	Hypothetical protein D910_08519, Dendroctonus ponderosae
			GO:0016798	MF	Hydrolase activity, acting on glycosyl bonds
			GO:0008152	BP	Metabolic process
			GO:0016787	MF	Hydrolase activity
			GO:0004767	MF	Sphingomyelin phosphodiesterase activity
TR64521|c0_g2	10.9	0E-00
TR66638|c0_g1	10.9	0E-00	GO:0003723	MF	RNA binding	Protein bicaudal C, Tribolium castaneum
			GO:0003676	MF	Nucleic acid binding
TR66019|c0_g2	10.8	9E-09
TR64213|c0_g1	10.7	2E-240
TR23777|c0_g3	10.7	1E-110
TR25403|c5_g2	10.6	3E-84
TR61942|c0_g1	10.5	2E-289
TR54080|c0_g1	10.4	0E-00	GO:0016787	MF	Hydrolase activity	Esterase, partial, Leptinotarsa decemlineata
TR54080|c0_g1	10.5	2E-14	GO:0008152	BP	Metabolic process
TR73901|c6_g2	10.5	2E-279
TR34425|c0_g3	10.6	5E-42
TR33970|c0_g1	10.4	0E-00	GO:0008233	MF	Peptidase activity	Cathepsin B-like cysteine protease, Callosobruchus maculates
			GO:0006508	BP	Proteolysis
			GO:0050790	BP	Regulation of catalytic activity
			GO:0004197	MF	Cysteine-type endopeptidase activity
			GO:0008234	MF	Cysteine-type peptidase activity
TR10926|c0_g1	10.3	0E-00
TR57758|c0_g1	10.5	6E-57
TR10388|c4_g1	10.3	0E-00	GO:0016787	MF	Hydrolase activity	Trypsin-1 isoform X1, Bactrocera cucurbitae
			GO:0008233	MF	Peptidase activity
			GO:0008236	MF	Serine-type peptidase activity
			GO:0004252	MF	Serine-type endopeptidase activity
			GO:0006508	BP	Proteolysis
TR36924|c3_g1	10.4	2E-43
TR27866|c0_g3	9.9	0E-00
TR52542|c5_g1	9.9	0E-00
TR3585|c2_g1	9.9	2E-02
TR7023|c1_g2	10.0	7E-55
TR18182|c4_g1	9.8	5E-160
TR28153|c0_g1	9.8	0E-00
TR62432|c0_g1	9.8	4E-195
TR9715|c0_g2	9.8	3E-04
TR68457|c0_g1	9.4	3E-33	GO:0046983	MF	Protein dimerization activity	Hairy, Leptinotarsa decemlineata
			GO:0006355	BP	Regulation of transcription, DNA-templated
			GO:0003677	MF	DNA binding
TR29717|c0_g3	9.7	0E-00
TR12240|c0_g1	9.6	9E-71
TR68765|c0_g2	9.5	0E-00	GO:0055114	BP	Oxidation–reduction process	Salicyl alcohol oxidase-like protein, Phaedon cochleariae
			GO:0050660	MF	Flavin adenine dinucleotide binding
TR10416|c0_g2	9.4	5E-26
TR58271|c0_g2	9.5	3E-158
TR36192|c2_g2	9.5	3E-269	GO:0006508	BP	Proteolysis	Cathepsin B-like proteinase, Triatoma vitticeps
			GO:0004197	MF	Cysteine-type endopeptidase activity
			GO:0016787	MF	Hydrolase activity
			GO:0008233	MF	Peptidase activity
			GO:0050790	BP	Regulation of catalytic activity
			GO:0008234	MF	Cysteine-type peptidase activity
TR16870|c0_g2	9.5	4E-21
TR45052|c0_g1	9.5	4E-38
TR24355|c0_g1	9.4	5E-244

CC, cellular component; MF, molecular function; BP, biological process.

Previous studies using proteomic analyses identified 98 seminal fluid proteins in C. maculatus (Goenaga et al. 2015; Bayram et al. 2016). The genes encoding these Sfps were overrepresented among the subtly male-biased genes in the abdomen. Although some of Sfp encoding genes were extremely male-biased (>10 logFC), several were also female-biased (supplementary fig. S3, Supplementary Material online). In line with our finding, many Sfp gene transcripts are expressed in the reproductive tracts of both males and females in a tiger mosquito Aedes albopictus (Boes et al. 2014). This shows that despite their important reproductive function in males, Sfp’s also play important roles in females. Indeed, it is possible that the manner by which some of these seminal fluid proteins act in mated females is by providing an extra dose of proteins that are produced endogenously within females (Arnqvist and Rowe 2005).

Male-biased genes in the head and thorax were enriched with genes involved in visual and chemical perception and neurotransmitter transport. These processes likely reflect the importance of detection of these sensory stimuli for males as part of mate searching, courtship and pre-mating competition. For example, two odorant binding genes were among the 50 most male-biased genes and showed over 25- and 10-fold higher expression in the males (odorant receptor 124, pheromone-binding protein-related protein 3, respectively, table 4). We note that the virgin males and females used in this study were naïve and held in isolation, and therefore these results stem from innate sex differences rather than direct responses to social interactions.

Table 4

The 50 Most Sex-Biased Genes with Annotations. Male-Biased (log2FC < −1) Genes in the Head and Thorax

Gene ID	logFC	FDR	GO	Type	Term	BLAST hit
TR45430|c1_g5	−9.7	3E-08
TR43196|c0_g1	−8.6	3E-23	GO:0055085	BP	Transmembrane transport	Putative inorganic phosphate cotransporter, Tribolium castaneum
			GO:0016021	CC	Integral component of membrane
TR45058|c0_g1	−8.4	3E-09
TR60993|c0_g2	−8.2	5E-12
TR12288|c0_g1	−8.1	3E-14
TR45430|c1_g4	−7.7	4E-11
TR68274|c1_g3	−6.0	1E-44
TR7204|c0_g1	−6.0	2E-10
TR7035|c2_g4	−5.7	0.0003	GO:0005576	CC	Extracellular region	Attacin-like immune protein, Diabrotica virgifera virgifera
TR56173|c0_g1	−5.5	1E-11
TR48097|c0_g1	−5.5	0.0005
TR14545|c0_g1	−5.3	3E-23	GO:0008234	MF	Cysteine-type peptidase activity	Putative gut cathepsin L-like cysteine protease, Callosobruchus maculates
			GO:0016787	MF	Hydrolase activity
			GO:0008233	MF	Peptidase activity
			GO:0006508	BP	Proteolysis
TR57561|c0_g2	−5.1	1E-19	GO:0008152	BP	Metabolic process	Esterase, Leptinotarsa decemlineata
			GO:0016787	MF	Hydrolase activity
TR34187|c0_g1	−4.9	1E-19	GO:0006030	BP	Chitin metabolic process	Peritrophic matrix protein 2-C precursor, T. castaneum
			GO:0008061	MF	Chitin binding
			GO:0005576	CC	Extracellular region
TR52512|c0_g4	−4.9	2E-07
TR68287|c1_g1	−4.8	4E-21
TR33586|c1_g1	−4.8	6E-10
TR29249|c0_g1	−4.8	8E-51
TR7634|c0_g1	−4.7	5E-23	GO:0016787	MF	Hydrolase activity	Esterase, L. decemlineata
			GO:0008152	BP	Metabolic process
TR50077|c0_g2	−4.6	3E-17	GO:0005549	MF	Odorant binding	Odorant receptor 124, partial, T. castaneum
			GO:0016020	CC	Membrane
			GO:0007608	BP	Sensory perception of smell
			GO:0004984	MF	Olfactory receptor activity
			GO:0050911	BP	Detection of chemical stimulus involved in sensory perception of smell
TR48099|c0_g1	−4.5	3E-10	GO:0005975	BP	Carbohydrate metabolic process	Hypothetical protein YQE_09695, partial, D. ponderosae
			GO:0016020	CC	Membrane
			GO:0016798	MF	Hydrolase activity, acting on glycosyl bonds
TR28685|c0_g2	−4.4	0.0003	GO:0016740	MF	Transferase activity	Glutathione S-transferase D2, T. castaneum
			GO:0008152	BP	Metabolic process
			GO:0004364	MF	Glutathione transferase activity
TR56870|c0_g5	−4.4	0.0197
TR53421|c0_g1	−4.2	9E-22	GO:0003993	MF	Acid phosphatase activity	Venom acid phosphatase Acph-1-like, T. castaneum
			GO:0006470	BP	Protein dephosphorylation
TR32496|c0_g1	−4.2	0.0032
TR10366|c0_g1	−4.2	1E-24
TR7204|c0_g2	−4.2	1E-09
TR72160|c0_g1	−4.1	7E-25	GO:0004497	MF	Monooxyge se activity	Cytochrome P450 412a2, L. decemlineata
			GO:0005506	MF	Iron ion binding
			GO:0055114	BP	Oxidation–reduction process
			GO:0020037	MF	Heme binding
			GO:0046872	MF	Metal ion binding
			GO:0016491	MF	Oxidoreductase activity
			GO:0016705	MF	Oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen
TR33543|c0_g2	−4.0	0.00043	GO:0004222	MF	Metalloendopeptidase activity	Zinc metalloproteinase s-13, T. castaneum
			GO:0008233	MF	Peptidase activity
			GO:0008237	MF	Metallopeptidase activity
			GO:0016787	MF	Hydrolase activity
			GO:0008270	MF	Zinc ion binding
			GO:0006508	BP	Proteolysis
			GO:0046872	MF	Metal ion binding
TR29759|c1_g1	−4.0	2E-14
TR17569|c0_g1	−4.0	1E-09
TR64329|c5_g1	−4.0	3E-30
TR57575|c0_g1	−4.0	0.0120	GO:0005856	CC	Cytoskeleton	Hpothetical protein YQE_08828, partial, D. ponderosae
			GO:0005525	MF	GTP binding
			GO:0051258	BP	Protein polymerization
			GO:0000166	MF	Nucleotide binding
			GO:0043234	CC	Protein complex
			GO:0005737	CC	Cytoplasm
			GO:0008152	BP	Metabolic process
			GO:0003924	MF	GTPase activity
			GO:0005200	MF	Structural constituent of cytoskeleton
			GO:0007017	BP	Microtubule-based process
			GO:0005874	CC	Microtubule
TR37369|c3_g7	−3.9	1E-49	GO:0055114	BP	Oxidation–reduction process	Acyl-CoA Delta(11) desaturase-like, T. castaneum
			GO:0006629	BP	Lipid metabolic process
			GO:0006633	BP	Fatty acid biosynthetic process
			GO:0016491	MF	Oxidoreductase activity
			GO:0016021	CC	Integral component of membrane
TR7035|c2_g2	−3.9	8E-05	GO:0005576	CC	Extracellular region	Attacin-like immune protein, D. virgifera virgifera
TR1055|c0_g1	−3.9	0.00906	GO:0030145	MF	Manganese ion binding	Hypothetical protein D910_02561, D. ponderosae
			GO:0005737	CC	Cytoplasm
			GO:0008235	MF	Metalloexopeptidase activity
			GO:0004177	MF	Aminopeptidase activity
			GO:0004672	MF	Protein kinase activity
			GO:0006508	BP	Proteolysis
			GO:0019538	BP	Protein metabolic process
			GO:0005622	CC	Intracellular
			GO:0006468	BP	Protein phosphorylation
TR56982|c4_g2	−3.8	2E-20	GO:0016787	MF	Hydrolase activity	Putative glycosyl hydrolase, Chrysomela lapponica
			GO:0008152	BP	Metabolic process
			GO:0005975	BP	Carbohydrate metabolic process
TR34693|c1_g1	−3.8	8E-19	GO:0016787	MF	Hydrolase activity	Glycoside hydrolase family protein 28, C. maculatus
			GO:0008152	BP	Metabolic process
			GO:0071555	BP	Cell wall organization
			GO:0005975	BP	Carbohydrate metabolic process
			GO:0004650	MF	Polygalacturonase activity
			GO:0005576	CC	Extracellular region
TR52603|c0_g2	–3.7	0.02513
TR33725|c0_g1	−3.7	3E-11
TR28685|c0_g1	−3.7	2E-50	GO:0004364	MF	Glutathione transferase activity	Glutathione S-transferase D2, T. castaneum
			GO:0016740	MF	Transferase activity
			GO:0008152	BP	Metabolic process
TR37895|c0_g1	−3.7	8E-19
TR57761|c0_g1	−3.6	3E-23
TR10418|c0_g2	−3.5	3E-43
TR19023|c0_g2	−3.4	1E-27	GO:0003993	MF	Acid phosphatase activity	Hypothetical protein YQE_11758, partial, D. ponderosae
TR19023|c0_g2	−3.4	1E-27	GO:0006470	BP	Protein dephosphorylation
TR12918|c3_g8	−3.4	4E-07
TR10887|c0_g1	−3.4	7E-25
TR15912|c4_g11	−3.4	7E-06
TR38223|c0_g1	−3.3	2E-11	GO:0005549	MF	Odorant binding	Pheromone-binding protein-related protein 3, Zootermopsis nevadensis

CC, cellular component; MF, molecular function; BP, biological process.

In addition to whole genes that are sex-biased in expression, an additional level of sexual dimorphism is found at loci that produce multiple transcripts through alternative pre-mRNA splicing or the use of alternative transcription initiation and termination sites (McIntyre et al. 2006). Sex-specific alternative splicing is a major player in sex determination including both cellular differentiation of somatic tissues (Nagoshi et al. 1988; Penalva and Sanchez 2003) and the gonads (Gan et al. 2010). In addition, pervasive sexual dimorphism in the patterns of alternative splicing in sexually mature adults is also common in Drosophila (McIntyre et al. 2006; Telonis-Scott et al. 2009), but has overall been less well explored. Here, we examined sex-biased expression of 4,581 genes with multiple expressed transcripts in the reproductive tissues. In total, 28% of these genes showed evidence for sex-specific alternative splicing, which makes up 12% of the whole expressed transcriptome. Even in the somatic tissues, 7% of the multi-transcript genes were alternatively spliced between the sexes. These results support the growing evidence of the prevalence of sex-specific alternative splicing in adults (McIntyre et al. 2006; Telonis-Scott et al. 2009; Hartmann et al. 2011). Our data also suggests that sex-biased splicing is much more common in the reproductive and sex-specific tissues of the abdomen (supplementary fig. S4, Supplementary Material online), in line with our gene-level analysis and the patterns seen in Drosophila (Telonis-Scott et al. 2009). Metabolic processes were overrepresented among the alternatively spliced genes in both tissue types.

We also specifically examined the exon-specific expression patterns among the 98 Sfp genes of C. maculatus and found that 27% showed evidence of sex-biased splicing, suggesting that alternative splicing may be a rather common way of regulating distinct phenotypes in each sex. Most of these genes exhibited a combined expression pattern of male-biased and un-biased transcripts (supplementary table S4a, Supplementary Material online). A similar mixture of male- and un-biased expression has previously been reported for a gene bmarg encoding for an arginase most likely secreted into the seminal fluid of Bombyx mori silkworm, where bmarg-r is expressed only in the male reproductive organs, while bmarg-f is expressed uniformly in the fat body and muscle of both sexes (Nagaoka et al. 2011). Alternative splicing of a gene encoding for a prostate enzyme transglutaminase 4, which is secreted into semen, has also been observed in humans (Cho et al. 2010). More detailed analyses of the candidate sex-specifically spliced Sfp genes identified here will be necessary to dissect the regulatory patterns and function of the isoforms in each sex.

Mating Induces Subtle but Numerous Changes in Sex-Biased Expression

Perhaps unsurprisingly, mating affects the expression of more genes in females than in males, and therefore the absolute contribution of females to the plasticity in virgin sex-biased expression is greater than of males (fig. 1). However, this picture is more moderated when reproductive and somatic tissues are considered separately. Male mating response in the abdomen contributes proportionally more to the plasticity in virgin sex-bias than female response: 54% of the male response genes show a change in sex-bias due to mating in contrast to 38% of the female response genes (fig. 1). The opposite is true for the head and thorax, where only 19% of the male response genes and 32% of the female response genes show plasticity in sex-bias (fig. 1). The great majority of the sex-bias plasticity is subtle and does not change the direction of the virgin-state bias. Sex-bias is typically characterized in reproductively naïve adults for categorizing genes for molecular evolutionary analyses (Jiang and Machado 2009). This practice assumes that the sex-bias status of a gene in virgins is representative of its pattern of expression throughout the adult lifespan during which sex-specific selection may operate. A study of sex-biased expression throughout ontogeny has recently demonstrated that many genes show different evolutionary patterns depending on whether their sex-biased expression is conserved or stage-specific (Perry et al. 2014), suggesting the importance of considering the patterns of plasticity also in adults. Our results provide evidence that the sex-bias status is rather stable in the face of altered mating status and therefore offer general assurance for the way sex-bias is typically characterized in adults.

Studies of sex-biased expression typically focus on the molecular evolutionary patterns, whereas analyzing the functional roles of sex-biased genes in males and females has received less attention. Here, we attempted to bridge some of this gap by asking how sex-biased genes respond to mating in each sex. Despite the highly sex-specific nature of postmating changes in physiology and behavior, sex-biased genes were the primary targets of mating only in the female abdomen (57% of the responsive genes), although they were also more common than expected in the female head and thorax (16%). In both types of tissues, females pre-dominantly downregulated female-biased genes, and in the abdomen a great majority of the male-biased genes responded by increased expression (fig. 4). This result shows that females also require the products of male-biased genes in their postmating reproductive physiology. Despite the general finding that male-biased genes are also expressed in females, they are assumed to be subject to only male-specific selection (Ellegren and Parsch 2007). Our results instead suggest that many male-biased genes may also experience female-specific selection which raises the possibility that, at least in some cases, the faster regulatory divergence often documented for male-biased genes (Ellegren and Parsch 2007) may have resulted from sex-specific selection operating in both sexes rather than in males alone. In males, the sex-biased genes responded to mating in a reversed fashion compared with females: male-biased genes were downregulated in both tissue groups, and female-biased genes more commonly upregulated than expected in the head and thorax (fig. 4). Taken together, these results suggest that sex-specific phenotypes involve expression of both types of sex-biased genes in both sexes. This finding echoes a previous study on a wild turkey that shows how the degree of masculinity of males does not only depend on the expression of male- but also on the relative expression of female-biased genes (Pointer et al. 2013).

In line with several previous studies (Mack et al. 2006; Innocenti and Morrow 2009; Gomulski et al. 2012; Alfonso-Parra et al. 2016), we found that mating increases the expression of genes involved in metabolism of carbohydrates and lipids, and reduction of oxidative stress in the female abdomen (supplementary table S6a, Supplementary Material online). These are the same pathways that show higher expression in virgin males. For example, the great majority of the mating responsive redox genes in the female abdomen are male biased and show increased expression due to mating in females (supplementary fig. S5, Supplementary Material online). These patterns likely reflect a shift in metabolic requirements for males and females that depend on mating status. In C. maculatus, virgin males show a higher metabolic rate and are more active than virgin females (Berger et al. 2016), as male pre-mating activities are characterized by behaviors such as active mate searching, harassment of other males, vigorous chasing and mounting of females. In females, mating induces oogenesis and oviposition, including the search for suitable oviposition sites, no doubt increasing metabolic demands in females. Hence, the masculinization seen in females after mating is likely to large extent due to an up-regulation of the metabolic machinery in mated females. Some of the observed expression changes in metabolic and redox genes may also contribute to nourishing and protecting stored sperm from oxidative damage (Baer et al. 2009; Prokupek et al. 2009; Shaw et al. 2014). In contrast to Drosophila (McGraw et al. 2004; Innocenti and Morrow 2009), the female mating response in C. maculatus did not involve antimicrobial responses. This may reflect differences between species, but it is also possible that the timing of expression of immunity genes differs. Some immunity genes may also go undetected because of the difficulty in annotating rapidly evolving genes based on sequence similarity, such as those involved in immunity and reproduction.

Conclusions

Our results show that sexual dimorphism is a pervasive feature of the C. maculatus transcriptome that is not only limited to the reproductive tissues. Mating induces changes in the degree of sex-bias for thousands of genes, demonstrating the plastic nature of adult gene expression. Although a number of genes changed sex-bias classification, this was relatively rare (4% in head and thorax and 2% in abdomen) suggesting that characterizing sex bias from reproductively naïve individuals is generally predictive of their sex-bias status also after mating. On the other hand, our results highlight the difficulty of predicting how genes may be involved in encoding reproductive phenotypes based on their sex-bias status. The fact that regulatory changes in both types of sex-biased genes occurred in the mating response of both sexes, suggest that selection may not be limited to the sex showing overall higher expression of the focal gene. More studies should therefore aim to link sex-biased expression with reproductive phenotypes under sex-specific selection.

Supplementary Material

Supplementary data are available at Genome Biology and Evolution online.

Click here for additional data file.

Table 2

The 50 Most Sex-Biased Genes with Annotations. Female-Biased (log2FC >1) Genes in the Head and Thorax

Gene	logFC	FDR	GO	Type	Term	BLASTHit
TR45052|c0_g1	9.8	7E-27
TR7171|c0_g2	9.7	1.E 96	GO:0016491	MF	Oxidoreductase activity	Cytochrome P450 412a2, Leptinotarsa decemlineata
			GO:0005506	MF	Iron ion binding
			GO:0004497	MF	Monooxyge se activity
			GO:0046872	MF	Metal ion binding
			GO:0020037	MF	Heme binding
			GO:0055114	BP	Oxidation–reduction process
TR37635|c0_g1	9.5	2E-21
TR23129|c0_g2	9.4	4E-27
TR38999|c0_g1	9.4	1E-67	GO:0006468	BP	Protein phosphorylation	Hypothetical protein YQE_12596, Dendroctonus ponderosae
			GO:0005524	MF	ATP binding
			GO:0004672	MF	Protein kinase activity
TR74056|c0_g1	9.1	7E-20
TR20045|c4_g3	8.9	9E-173
TR72359|c0_g1	8.9	4E-05
TR33787|c0_g4	8.9	2E-21	GO:0016874	MF	Ligase activity	Hypothetical protein TcasGA2_TC010120, Tribolium castaneum
			GO:0008152	BP	Metabolic process
TR25835|c0_g2	8.8	1E-14
TR68289|c3_g2	8.7	1E-109
TR28694|c0_g1	8.7	7E-20
TR54960|c0_g3	8.7	4E-14	GO:0003676	MF	nucleic acid binding	Hypothetical protein TcasGA2_TC012080, T. castaneum
			GO:0003723	MF	RNA binding
			GO:0000166	MF	Nucleotide binding
TR58271|c0_g2	8.5	5E-11
TR35084|c0_g1	8.5	5E-18
TR29702|c2_g6	8.4	1E-61
TR18180|c0_g1	8.4	6E-56
TR44554|c0_g1	8.3	2E-76
TR18182|c4_g1	8.3	5E-54
TR12240|c0_g1	8.3	1E-53
TR18182|c9_g4	8.2	7E-49
TR44608|c4_g2	8.2	0.0003
TR28684|c0_g4	8.2	1E-12	GO:0004861	MF	Cyclin-dependent protein serine/threonine kinase inhibition	Protein HEXIM, T. castaneum
			GO:0000122	BP	Negative regulation of transcription from RNA polymerase II promoter
			GO:0017069	MF	snRNA binding
			GO:0071901	BP	Negative regulation of protein serine/threonine kinase
			GO:0005634	CC	Nucleus
			GO:0005737	CC	Cytoplasm
TR57758|c0_g1	8.1	9E-14
TR4301|c0_g2	7.9	5E-12
TR66638|c0_g1	7.9	1E-48	GO:0003676	MF	Nucleic acid binding	Protein bicaudal C, T. castaneum
			GO:0003723	MF	RNA binding
TR61957|c3_g2	7.8	0E+00
TR75409|c0_g1	7.8	2E-37	GO:0046872	MF	Metal ion binding	Uncharacterized protein LOC103314013, T. castaneum
TR62432|c0_g1	7.8	7E-22
TR69060|c0_g3	7.8	3E-14
TR59852|c0_g2	7.8	6E-61	GO:0016567	BP	Protein ubiquitination	Uncharacterized protein LOC663065 isoform X2, T. castaneum
			GO:0005634	CC	Nucleus
			GO:0016874	MF	Ligase activity
			GO:0046872	MF	Metal ion binding
			GO:0008270	MF	Zinc ion binding
			GO:0006511	BP	Ubiquitin-dependent protein catabolic process
			GO:0007275	BP	Multicellular organismal development
			GO:0004842	MF	Ubiquitin-protein transferase activity
TR44673|c0_g9	7.7	1E-16
TR44247|c0_g1	7.7	4E-50
TR55946|c0_g2	7.6	7E-75
TR45396|c0_g1	7.5	9E-24
TR1169|c3_g9	7.5	8E-23
TR54908|c0_g1	7.5	8E-107
TR49895|c0_g1	7.4	1E-11	GO:0003743	MF	Translation initiation factor activity	Eukaryotic translation initiation factor 2 subunit 1, T. castaneum
			GO:0003723	MF	RNA binding
			GO:0006413	BP	Translation l initiation
			GO:0003676	MF	Nucleic acid binding
			GO:0005850	CC	Eukaryotic translation initiation factor 2 complex
TR57251|c0_g1	7.4	7E-55
TR18182|c9_g3	7.4	9E-23
TR71884|c0_g1	7.2	5E-39	GO:0005634	CC	Nucleus	G2/mitotic-specific cyclin-A, T. castaneum
TR44662|c0_g1	7.1	4E-36	GO:0020037	MF	Heme binding	Hypothetical protein TcasGA2_TC000751, T. castaneum
			GO:0004601	MF	Peroxidase activity
			GO:0006979	BP	Response to oxidative stress
			GO:0055114	BP	Oxidation–reduction process
TR36341|c0_g2	7.0	1E-16
TR63881|c0_g1	7.0	3E-20
TR29995|c3_g9	6.9	4E-16
TR64757|c0_g1	6.9	3E-05
TR8233|c0_g2	6.8	3E-33
TR1205|c0_g3	6.7	8E-16	GO:0016021	CC	Integral component of membrane	Hypothetical protein YQE_10313, D. ponderosae
			GO:0022857	MF	Transmembrane transporter activity
			GO:0006810	BP	Transport
			GO:0055085	BP	Transmembrane transport
			GO:0005215	MF	Transporter activity
			GO:0022891	MF	Substrate-specific transmembrane transporter activity
			GO:0016020	CC	Membrane
TR61942|c0_g1	6.7	0.0001
TR68765|c0_g1	6.6	4.E-160	GO:0050660	MF	Flavin adenine dinucleotide binding	Hypothetical protein TcasGA2_TC015713, T. castaneum
			GO:0055114	BP	Oxidation–reduction process

CC, cellular component; MF, molecular function; BP, biological process.

Table 3

The 50 Most Sex-Biased Genes with Annotations. Male-Biased (log2FC < −1) Genes in the Abdomen

Gene ID	logFC	FDR	GO	Type	Term	BLAST hit
TR6663|c0_g1	−16.5	0E+00
TR27296|c0_g1	−14.8	0E+00
TR2855|c0_g1	−14.4	0E+00
TR19004|c0_g5	−14.4	0E+00
TR64284|c2_g6	−14.3	0E+00	GO:0016491	MF	Oxidoreductase activity	Unknown, Dendroctonus ponderosae
			GO:0055114	BP	Oxidation–reduction process
TR6663|c0_g5	−14.3	0E+00
TR71857|c0_g2	−13.9	0E+00	GO:0005737	CC	Cytoplasm	Unknown, Dendroctonus ponderosae
			GO:0004177	MF	Aminopeptidase activity
			GO:0030145	MF	Manganese ion binding
			GO:0008235	MF	Metalloexopeptidase activity
			GO:0019538	BP	Protein metabolic process
			GO:0006508	BP	Proteolysis
TR27305|c0_g1	−13.8	0E+00
TR47060|c0_g1	−13.6	0E+00
TR34460|c0_g3	−13.6	0E+00
TR64463|c0_g1	−13.4	0E+00
TR51551|c0_g1	−13.4	0E+00
TR17032|c0_g1	−13.4	0E+00	GO:0005509	MF	Calcium ion binding	Unknown, Dendroctonus ponderosae
TR43223|c0_g1	−13.3	0E+00
TR49381|c0_g1	−13.3	1E-276
TR29424|c0_g1	−13.2	7E-270	GO:0005737	CC	Cytoplasm	Unknown, Dendroctonus ponderosae
			GO:0005622	CC	Intracellular
			GO:0008235	MF	Metalloexopeptidase activity
			GO:0004177	MF	Aminopeptidase activity
			GO:0030145	MF	Manganese ion binding
			GO:0019538	BP	Protein metabolic process
			GO:0006508	BP	Proteolysis
TR7979|c0_g2	−13.2	7E-251
TR55513|c0_g1	−13.1	0E+00
TR29301|c0_g2	−13.0	0E+00
TR73858|c0_g1	−13.0	0E+00
TR44273|c0_g1	−13.1	3E-150	GO:0005524	MF	ATP binding	Heat shock cognate 71 kDa protein-like isoform X1, Tribolium castaneum
			GO:0000166	MF	Nucleotide binding
TR66585|c0_g1	−13.0	4E-268	GO:0005634	CC	Nucleus	Hypothetical protein YQE_11918, partial, Dendroctonus ponderosae
			GO:0006334	BP	Nucleosome assembly
TR17086|c0_g1	−13.0	0E+00
TR70731|c0_g1	−13.0	0E+00
TR16878|c0_g1	−12.9	0E+00
TR27262|c0_g1	−12.9	0E+00
TR22315|c0_g1	−12.7	0E+00
TR28238|c1_g4	−12.8	2E-282
TR30376|c0_g1	−12.7	0E+00
TR38674|c0_g1	−12.9	6E-234	GO:0005506	MF	Iron ion binding	Unknown, Aedes aegypti
			GO:0016491	MF	Oxidoreductase activity
			GO:0020037	MF	Heme binding
			GO:0046872	MF	Metal ion binding
			GO:0016705	MF	Oxidoreductase activity
			GO:0004497	MF	Monooxygenase activity
			GO:0055114	BP	Oxidation–reduction process
TR12307|c0_g1	−12.9	1E-282
TR41652|c0_g1	−12.7	0E+00
TR18203|c0_g1	−12.6	0E+00
TR39204|c1_g1	−12.7	0E+00
TR74166|c0_g2	−12.6	0E+00	GO:0005737	CC	Cytoplasm	Unknown, Dendroctonus ponderosae
			GO:0016301	MF	Kinase activity
			GO:0016740	MF	Transferase activity
			GO:0016773	MF	Phosphotransferase activity
			GO:0004335	MF	Galactokinase activity
			GO:0005524	MF	ATP binding
			GO:0000166	MF	Nucleotide binding
			GO:0008152	BP	Metabolic process
			GO:0016310	BP	Phosphorylation
			GO:0006012	BP	Galactose metabolic process
			GO:0046835	BP
TR49828|c0_g2	−12.6	0E+00	GO:0016491	MF	Oxidoreductase activity	Hypothetical protein TcasGA2_TC010329, Tribolium castaneum]
			GO:0008152	BP	Metabolic process
			GO:0055114	BP	Oxidation–reduction process
TR59474|c0_g1	−12.7	6E-259
TR3816|c0_g1	−12.7	6E-308	GO:0005874	CC	Microtubule	Hypothetical protein D910_07655, Dendroctonus ponderosae
			GO:0005871	CC	Kinesin complex
			GO:0000166	MF	Nucleotide binding
			GO:0005524	MF	ATP binding
			GO:0008017	MF	Microtubule binding
			GO:0003777	MF	Microtubule motor activity
			GO:0008152	BP	Metabolic process
			GO:0007018	BP
TR14537|c0_g1	−12.5	0E+00
TR69544|c0_g1	−12.6	6E-261	GO:0008270	MF	Zinc ion binding	Hypothetical protein YQE_01746, partial, Dendroctonus ponderosae
			GO:0004181	MF
			GO:0006508	BP	Proteolysis
TR22773|c0_g1	−12.6	2E-239
TR36788|c0_g2	−12.6	0E+00
TR810|c0_g2	−12.6	4E-207
TR57874|c0_g3	−12.5	0E+00	GO:0004252	MF	Serine-type endopeptidase activity	Hypothetical protein D910_08079, Dendroctonus ponderosae
			GO:0006508	BP	Proteolysis
TR47457|c0_g1	−12.5	7E-264
TR15902|c0_g1	−12.4	0E+00
TR15343|c0_g2	−12.5	3E-312
TR31224|c1_g1	−12.4	0.0E+00	GO:0004252	MF	Serine-type endopeptidase activity	Hypothetical protein D910_03668, Dendroctonus ponderosae
			GO:0006508	BP	Proteolysis
TR24113|c2_g5	−12.5	0E+00
TR38284|c0_g1	−12.4	2E-228	GO:0008177	MF	Succinate dehydrogenase activity	Succinate dehydrogenase flavoprotein subunit, mitochondrial like, Tribolium castaneum
			GO:0006099	BP	Tricarboxylica acid cycle
			GO:0016020	CC	Membrane
			GO:0055114	BP	Oxidation–reduction process
			GO:0016491	MF	Oxidoreductase activity
			GO:0022900	BP	Electron transport chain
			GO:0050660	MF	Flavin adenine dinucleotide binding
			GO:0005743	CC	Mitochondrial inner membrane

CC, cellular component; MF, molecular function; BP, biological process.