Effects of oxidative stress on sex-specific gene expression in the copepod Tigriopus californicus revealed by single individual RNA-seq
Oxidative stress reflects the imbalance of pro-oxidants and antioxidants. Prolonged oxidative stress can induce cellular damage, diseases and aging, and the effects may be sex-specific. Tigriopus californicus has recently been proposed as an alternative model system for sex-specific studies due to the absence of sex chromosomes. In this study, we used comparative transcriptomic analyses to assess sex-specific transcriptional responses to oxidative stress. Male and female individuals were maintained separately in one of three treatments: 1) control conditions with an algae diet, 2) pro-oxidant (H2O2) conditions with an algae diet or 3) decreased antioxidant conditions (reduced carotenoid due to a yeast diet). Single individual RNA-seq was then conducted for twenty-four libraries using Ligation Mediated RNA sequencing (LM-Seq). Variance in gene expression was partitioned into 62.3% between sexes, 26.85% among individuals and 10.85% among treatments. Within each of the three treatments, expression was biased toward females. However, compared to the control treatment, males in both pro-oxidant and decreased antioxidant treatments differentially expressed more genes while females differentially expressed fewer genes but with a greater magnitude of fold change. As the first study of copepods to apply single individual RNA-seq, the findings will contribute to a better understanding of transcriptomic variation among individuals as well as sex-specific response mechanisms to oxidative stress in the absence of sex chromosomes.
Introduction<!>Copepod treatments, maintenance and sampling<!>RNA extraction, library preparation and transcriptome sequencing<!>Reads mapping and differential expression analysis<!>Partitioning of gene expression variation<!>Single individual RNA-seq<!>Gene expression analysis<!>Functional enrichment analysis<!>Candidate genes responding to oxidative stress<!>Discussion<!>Conclusion
<p>A free radical is any chemical species that contains one or more unpaired electrons (Halliwell and Gutteridge, 2015). The diatomic oxygen molecule is a free radical with two unpaired electrons which have parallel spins in separate orbits (Halliwell and Gutteridge, 2015). Thus, if oxygen attempts to accept a pair of electrons to oxidize another molecule, both electrons have to be in the parallel spin as well to fit in the orbits (Halliwell and Gutteridge, 2015). Since a pair of electrons have opposite spins in the same orbit, oxygen prefers to accept the electrons one at a time, and this sequential reduction of oxygen leads to the generation of various types of Reactive Oxygen Species (ROS), such as superoxide and hydrogen peroxide (Halliwell and Gutteridge, 2015). ROS could be derived from exogenous sources such as pollutants or UV radiation, as well as endogenous ROS production, such as through the mitochondrial respiratory chain (Monaghan et al., 2009). At the same time, organisms have a series of enzymatic and non-enzymatic defense systems to eliminate the harmful ROS (Lesser, 2006). The imbalance between the production of ROS and the elimination by antioxidant systems would result in a disturbance of redox status and cause a state of oxidative stress (Lushchak, 2011). Oxidative stress may cause damage to lipids, proteins and DNA, and subsequently lead to cellular damage, diseases and aging (Lesser, 2006). A positive correlation has been observed between oxidative stress resistance and longevity by using chemicals such as paraquat or ethanol to introduce oxidative stress (Finkel and Holbrook, 2000; Vermeulen et al., 2005; Niveditha et al., 2017). Moreover, the effect of oxidative stress may be sex-specific. Under oxidative stress female Drosophila melanogaster exhibit lower ROS, higher levels of antioxidant enzymes, lower mortality and longer lifespan than males (Niveditha et al., 2017). Similarly, female Wistar rats have been shown to have higher antioxidant gene expression, lower oxidative damage on mitochondrial DNA, and longer lifespan than males (Borrás et al., 2003).</p><p>This study focuses on Tigriopus californicus, an intertidal copepod species that makes an ideal laboratory organism due to its short generation time, amenability to laboratory culture and extensive genetic resources. A high-quality reference genome has recently been published (Barreto et al., 2018), providing valuable information for further gene-based functional studies. The species is a particularly interesting alternative model for understanding sex differences because it does not have sex chromosomes. Instead, sex determination is polygenic with a small environmental component (Ar-rushdi, 1962; Voordouw and Anholt, 2002; Alexander et al., 2014; Alexander et al., 2015). This means that T. californicus provides a simplified system for understanding sexual dimorphism, as sex-specific effects cannot be attributed to sex chromosomes, where expression is complicated by dosage compensation, a process organisms utilize to equalize gene expression between sexes with different copy numbers of sex chromosomes (e.g. Disteche, 2012). T. californicus is also free from the complication of genes that have different fitness effects in males and females being moved on or off sex chromosomes, thereby adjusting the level of sexual conflict (Spigler et al., 2011; Hill, 2013; Dean et al., 2014).</p><p>Next-generation RNA sequencing (RNA-seq) is based on high-throughput sequencing techniques and ideal for profiling of large quantities of expression data before and after stress treatments, where regulation-related genes and pathways for stress tolerance could be identified. Since T. californicus is small (~1 mm), previous RNA-seq work pooled a large number of individuals, sometimes at an equal sex ratio, to increase the total amount of RNA for transcriptome sequencing (Barreto et al., 2011; Schoville et al., 2012; Pereira et al., 2016; Lima and Willett, 2017; DeBiasse et al., 2018; Graham and Barreto, 2019). Pooling reveals average expression levels but obscures variation among individuals. This pooling design relies heavily on the biological averaging hypothesis, which refers to the assumption that RNA levels average out when pooled (Kendziorski et al., 2005). However, this assumption may not hold in the presence of pooling bias, which is the value differences measured between pool and the mean of corresponding individuals (Kendziorski et al., 2003; Kendziorski et al., 2005; Mary-Huard et al., 2007; Rajkumar et al., 2015). The overall pooling bias results from imperfect averaging of the individuals in the pool and also the differences caused by log transformation of expression values during the normalization process (Kendziorski et al., 2005; Mary-Huard et al., 2007). Both simulation and real data suggested that pooling samples for gene expression analysis will increase the pooling bias as well as false discovery rates (Mary-Huard et al., 2007; Rajkumar et al., 2015). Moreover, within population variation will be reduced or neglected in the analysis of a pooled sample, decreasing the ability to identify valuable information at the individual subject level (Kendziorski et al., 2005; Rajkumar et al., 2015).</p><p>In addition to the bias from pooling multiple individuals, previous work on gene expression in T. californicus has ignored sex differences. It has been well documented that sex differences exist across taxa in stress tolerance, lifespan and gene expression (Marotti et al., 2010; Archer et al., 2013; Mank, 2017; Niveditha et al., 2017; Tower, 2017). In T. californicus, males are generally less tolerant to stresses (Willett, 2010; Kelly et al., 2012; Foley et al., 2019). Whether these sex differences in stress tolerance are associated with different transcriptomic responses is still not clear, and thus it seems indispensable to sequence transcriptomes of males and females separately. The development of cost-effective Ligation Mediated RNA sequencing (LM-Seq) protocol (Hou et al., 2015) makes it possible to prepare sequencing libraries for single male or female T. californicus individuals and subsequently detect their gene expression patterns under a certain stress.</p><p>Oxidative stress could be induced by various abiotic factors such as environmental and nutritional conditions. Hydrogen peroxide (H2O2) is a pro-oxidant which has been widely used to induce oxidative stress in aquatic invertebrates (e.g. Seo et al., 2006; Rhee et al., 2011). It can easily diffuse through cell membranes and increase ROS levels directly (Abele-Oeschger et al., 1997). In contrast, a nutritional yeast diet indirectly induces oxidative stress by decreasing the input of carotenoids obtained from an algae diet (Hill et al., 2018), thereby weakening antioxidant defenses (Ndhlala et al., 2010). Yeast fed copepods have been shown to be more sensitive to ultra-violet radiation and copper, suggesting the important roles of carotenoids in protecting cells against oxidation (Davenport et al., 2004; Caramujo et al., 2012). This study uses H2O2 as a pro-oxidant to increase the production of ROS and also uses a yeast diet to induce oxidative stress by decreasing the input of antioxidants. RNA-seq will be employed to assess sex-specific transcriptomic responses to these two different types of oxidative stress, and to evaluate possible mechanisms of tolerance.</p><!><p>Copepods were collected from Santa Catalina Island, USA (33°44'66"N, 118°48'50"W), and maintained in the laboratory prior to the start of the experiment. Three culture media were prepared for three treatments (1) control conditions: 0.1g ground Tetramin flakes (Tetra, Germany) and 0.1g powdered Spirulina (Nutrex Hawaii, USA) per 1L three-times-filtered seawater (from Santa Catalina island); (2) pro-oxidant conditions: 0.1g Tetramin and 0.1g Spirulina per 1L filtered seawater with 2 mM H2O2 (30% solution, EMD Millipore, Germany); and (3) decreased antioxidant conditions: 0.2g ground nutritional yeast (Bob's Red Mill, USA) per 1L filtered seawater. Each treatment group was started with 50 gravid females with red or orange eggs because these eggs are about to hatch in a short time (Hawkins, 1962). Cultures were kept in a 25°C incubator with a 12-hour light:dark cycle. Adults were removed from dishes the next day. Juveniles were moved to a new petri dish with the corresponding culture medium every week until males and females were distinguishable. Males and females were separated into different dishes and maintained in the same corresponding culture medium. Under the decreased antioxidant conditions, copepods lost the orange-red coloration due to the lack of catenoids.</p><p>Samples were only collected during the light cycle. Twenty-four adult individuals (4 replicates * (2 treatments + 1 control) * 2 sexes) were sampled 35 days after the start of the experiment. Each individual was rinsed with autoclaved three-times-filtered seawater and then transferred to a piece of filter paper. The egg sacs of individual females were removed using probes. Each individual was then immediately transferred into a round bottom nuclease-free microtube with 30–50 1.0 mm diameter zirconia/silica beads (BioSpec Products, USA). The sample was flash-frozen in liquid nitrogen and stored at −80 °C until RNA extraction.</p><!><p>RNA extraction was conducted using the Direct-zol RNA MicroPrep Plus kit according to the manufacturer's instructions (Zymo Research, USA). Briefly, each sample was lysed in 300 μL TRIzol Reagent (Ambion, USA) and homogenized using a TissueLyser (Qiagen, USA). Two shaking steps were carried out at 25 Hz for 2 minutes to ensure uniform disruption and homogenization of samples. An equal volume of ethanol (95–100%, 300 μL) was added to the lysed sample, and then the mixture was transferred into a Zymo-Spin IC Column from the kit. After in-column DNase I treatment and three RNA wash steps, 15 μL nuclease-free water was used to elute RNA from the column. RNA was stored at −80 °C until RNA sequencing library preparation.</p><p>Twenty-four libraries were constructed based on the LM-Seq protocol, which has been shown to be efficient for library preparation with RNA input as low as 10 ng (Hou et al., 2015). Briefly, mRNA was isolated from the total extracted RNA using the NEBNext Poly(A) mRNA Magnetic Isolation Module (NEB, USA), and then fragmented with heat (85 °C, 7 min). The fragmented mRNA was reverse-transcribed with SmartScribe reverse transcriptase (Clontech, USA) using a random hexamer oligo. After RNA removal by RNase A (Epicentre, USA) and RNase H (NEB), a partial Illumina 5' adaptor was ligated to single-stranded cDNA using T4 RNA ligase 1 (NEB). Ligated cDNA was then amplified by 20 PCR cycles using FailSafe PCR system (Epicentre), during which the full-length Illumina adaptor and index primer were included to generate the final library for Illumina sequencing. Additional dual size-selection process was conducted using AMPure XP beads (Beckman Coulter, USA). The library was quantified with the Qubit dsDNA HS Assay Kit (Invitrogen, USA) and quality-assessed by Bioanalyzer 2100 system (Agilent, USA). All the primer sequences used in the library preparation were listed in Hou et al., 2015. Libraries were pooled and sequenced on one lane of the Illumina HiSeq 4000 platform to obtain paired-end 150 base pairs (bp) reads at Fulgent Genetics (CA, USA).</p><!><p>Raw reads were filtered with the parameters of base quality ≥ 20 and trimmed read length ≥ 36 bp by Trimmomatic v0.32 (Bolger et al., 2014), and trimmed reads were evaluated by FastQC v0.11.4 (Andrews, 2010). Mapping strategies for T. californicus need to contend with the high level of transcriptomic divergence among populations (nuclear divergence 0.7% to 2.8%; Pereira et al., 2016). Because a well assembled and annotated reference genome for the San Diego populations has recently become available (Barreto et al., 2018), we chose to map our Santa Catalina Island population reads to this reference genome and subsequently extract the expression levels of each annotated gene. Because mitochondrial divergence between T. californicus populations is particularly large (9.5% ~ 26.5%) (Barreto et al., 2018), we replaced the San Diego mitochondrial genome with the corresponding mitochondrial genome (Santa Catalina Island) for reads mapping. HISAT2 v2.1.0 (Kim et al., 2015) was used to map all the trimmed reads to the modified reference genome, and StringTie v1.3.4d (Pertea et al., 2015) was utilized for quantification of expression levels of all annotated genes.</p><p>To test our mapping strategy, we used published sequencing reads from multiple populations (Barreto et al., 2014; Lima and Willett, 2017; DeBiasse et al., 2018) and mapped them to the San Diego reference genome (Table S1). Reads from San Diego (SD) samples yielded alignment rates ranging from 41% to 96%. Reads from non-SD populations (Santa Cruz and Bodega Bay) yielded alignments of 55% - 75%. Although the SD and non-SD data came from different studies, the overlap in mapping efficiencies support our strategy of mapping Santa Catalina Island reads to the SD reference genome. Since 20 PCR amplification cycles were used for library preparation, non-T. californicus contamination will be amplified as well as T. californicus sequences. Non-mapped reads were extracted and BLASTX was used to search the reads against non-redundant protein database to identify possible contamination for future improvement of library preparation.</p><p>DESeq2 (Love et al., 2014) was used to conduct differential expression analysis between treated groups (H2O2 treated group and yeast fed group) and control groups within both sexes, respectively. Pre-filtering was conducted to keep genes that have at least ten total reads across all the samples for better software performance. Transcripts with adjusted P-value < 0.05 and absolute values of fold change greater than 2 were determined as differentially expressed genes (DEGs) for subsequent analysis. Magnitudes of fold change of all DEGs under both treatments were used for comparisons between sexes by Mann-Whitney U test in R. Blast2GO v4.0.7 was used with default settings for the over-representation analysis of gene ontology (GO) terms among the DEGs using the whole genome as the reference.</p><!><p>We formulated a model derived from previously published studies (Giger et al., 2006; Storey et al., 2007; Hughes et al., 2015) to estimate the proportion of gene expression variation explained by factors of interest. Briefly, a hierarchical analysis of variance (ANOVA) was performed on the normalized gene expression levels for each gene separately. The nested ANOVA model for each gene was assumed as:
Yijk=μ+Ti+Sj(i)+Rk(ij)+eijk
where Yijk is the normalized expression level for the given gene in the kth replicate (R) which is nested within jth sex (S), which is nested within ith treatment (T). μ represents the mean expression level for the given gene. The total variance in gene expression for the given gene was calculated as σtotal2=σT2+σS2+σR2, and the proportion of variance explained by each variable was measured as the ratio of the variance due to the corresponding variable to the total variance: σ2/σtotal2.</p><p>After normalization by DESeq2, genes with half of the individuals having less than 1 read per million mapped reads were filtered out. A total of 9,036 genes were kept and their normalized expression values were transformed by log2 for further partition of variance by R package variancePartition version 1.12.1 (Hoffman and Schadt, 2016). A weighted average proportion of variance due to each effect over all genes was calculated.</p><!><p>RNA-seq was conducted for 24 individuals (4 replicates * (2 treatments + 1 control) * 2 sexes), yielding a total of ~354 million raw read pairs with the length of 150 bp. Low-quality bases (base quality < 20), short reads (trimmed read length < 36 bp) and ambiguous nucleotides (marked with N) were removed. After trimming, 4.6–18.5 million read pairs per individual were retained (Table S2). All the clean reads were mapped to annotated genes in the modified reference genome (See detailed descriptions of modification in Materials and Methods 2.3), and the alignment rate ranges from 19% to 54% (Table S2). Multiple contaminations in the sequencing reads were identified from bacteria, fungi, algae and amplified genomic DNA from other species studied in the lab (Table S3). The read sequences were deposited at the National Center for Biotechnical Information (NCBI) Sequence Read Archive (SRA) under the accession numbers SRR8870100-SRR8870123.</p><!><p>Principal component analysis (PCA) was used to visualize relationships among the 24 samples in two dimensions (Figure 1A). The two sexes were separated by the first principal component (PC1, ~45%), except for one yeast fed female and one H2O2 treated male, while the addition of the second principal component (PC2, ~14%) was able to successfully separate both sexes under each condition. ANOVA revealed that, on average, 62.3%, 26.85% and 10.85% of the gene expression variation was explained by variation between sexes, among individuals and among treatments, respectively, consistent with the pattern observed in the PCA. Hierarchical clustering with average expression values (Figure 1B) also shows that sex differences are greater than treatment differences, as the three female treatments form a single cluster. In addition, for both sexes the H2O2 treatment clusters with the control to the exclusion of the yeast fed treatment.</p><p>The two sexes were compared in the control and both treatment groups. In the control group, 775 genes were found to be sex-biased (> 2-fold expression difference), accounting for 4.95% of all genes (15,646 genes; Figure 2A, Table S4). The number of sex-biased genes was 2,706 (17.30%) in the H2O2 treated group (Figure 2B, Table S5) and 416 (2.66%) in the yeast fed group (Figure 2C, Table S6). Of these sex-biased genes in the treatment groups, only 475 (H2O2) and 122 (yeast diet) were also sex-biased in the control group, while the two treatment groups shared 279 genes with the same sex-bias direction. A total of 110 genes (0.7%) were found to be sex-biased across all three groups, which indicated that 169 genes became sex-biased because of the sex-specific transcriptomic responses to oxidative stress. Additionally, as Figure 2 shows, greater than 89% of the sex-biased genes are female-biased in each of the three groups.</p><p>To understand sex-specific and stress-specific responses to oxidative stress, pairwise comparisons were conducted between treatment and control groups within each sex to identify differentially expressed genes (DEGs) (Table S7–S10). Both sex differences and stress differences in gene expression were observed based on the number of DEGs. As shown in Figure 3, males differentially expressed more genes than females under both treatments (1,346 versus 44 in the H2O2 treated group, and 577 versus 24 in the yeast fed group). Within each sex, the H2O2 treated group differentially expressed more genes than the yeast fed group (44 versus 24 in females, and 1,346 versus 577 in males). In terms of the proportion of shared DEGs between contrasts, the two treatment groups shared a greater proportion of DEGs within each sex than those shared by both sexes within each treatment (Figure 4A). Taken together, sex is the major factor affecting the overall expression differences across all the groups. The magnitudes of fold change of all DEGs were collected for comparisons between males and females. Female groups have absolute means of 15.9 and 9.6, while male groups display absolute means of 6.3 and 3.6 under the H2O2 treatment (p=1.64e-14) and the yeast fed treatment (p=7.09e-04), respectively (Figure 4B). Consequently, in response to oxidative stress, males tend to employ more DEGs while females tend to differentially express fewer genes but with a greater magnitude of fold change.</p><!><p>Over-represented GO terms were determined for DEGs in each group. No enriched GO terms were detected in females, possibly because there were so few DEGs (44 under the H2O2 treatment and 24 under the yeast fed treatment). In the H2O2 treated male group, protein binding (GO: 0005515) was down-regulated, and mRNA splicing (GO: 0000398) was found to be up-regulated (Table 1). In the yeast fed male group, a series of protein synthesis, ATP synthesis and ATP generation-related GO terms were down-regulated (Table 1). No enriched GO terms were identified between sexes.</p><!><p>Hydrogen peroxide is commonly used to induce oxidative stress, and many oxidative stress responding genes have been identified (Finkel and Holbrook, 2000; Scandalios, 2005; Seo et al., 2006; Rhee et al., 2011). In the H2O2 treated male group, differential expression was found for a number of genes known to respond to hydrogen peroxide including 2 thioredoxin domain containing proteins, 2 heat shock proteins, 1 oxidoreductase, 2 mitogen activated protein kinases, 1 superoxide dismutase, 1 glutathione S-transferase and 1 thioredoxin (Table S7). In contrast, DEGs in the H2O2 treated female group did not include known candidates (Table S8). Similarly, 1 heat shock protein, 2 oxidoreductases and 2 glutathione S-transferases were found to be differentially expressed in the yeast fed male group (Table S9), while no known candidates were found in the yeast fed female group (Table S10). Those DEGs shared by both sexes in each group have no obvious direct role in responding to oxidative stress (Table S11 and S12).</p><p>DEGs shared by all the contrasts indicate a central role in responding to oxidative stress and may be used as biomarkers or indicators for oxidative stress-related studies. Only one gene (TCALIF_08006 on chromosome 5) was found to be differentially expressed in all the contrasts and this gene was up-regulated under both treatments within each sex (Figure S1). Since this gene has no annotation information and the only hit in the NCBI database is an uncharacterized protein from the copepod Eurytemora affinis, de novo assembled transcriptomes from another two T. californicus populations (Kelly et al., 2017) were searched to validate the accuracy and completeness of this gene. The protein sequence of this gene is 100% matched to the transcripts found in those two transcriptomes, indicating the reliability and validity of our finding.</p><!><p>We used single individual RNA-seq to determine transcriptome-wide gene expression patterns responding to oxidative stress in both sexes of the copepod T. californicus. As the first study in this species to assess gene expression in individuals, our results provide new information on among-individual variance in stress response, which is the raw material for natural selection. Comparisons between sexes demonstrated sex-specific transcriptomic responses to oxidative stresses, where males and females expressed different genes as well as different patterns of fold change. These results will contribute to a better understanding of sex-specific response mechanisms to oxidative stress, particularly in the absence of sex chromosomes where interpretation is not complicated by dosage compensation and shuffling of sexually antagonistic genes on or off of sex chromosomes.</p><p>Our sequencing data analysis only yielded alignment rates of 19% - 54%, which are lower than those of published data using pooled samples. Since RNA extraction of a single individual copepod only yields ~10 ng total RNA, the LM-Seq protocol utilizes 20 cycles of PCR amplification to generate sufficient product for sequencing (Hou et al., 2015). It is noted that PCR does not amplify all sequences equally well (Parekh et al., 2016), which may add PCR noise or enlarge the proportion of contamination and hence affect the accuracy of RNA-seq analysis. In fact, contaminations from multiple sources were detected in the reads, including bacteria, fungi, algae from the culturing and amplified genomic DNA from a different lab project. Therefore, even greater attention should be paid to cleaning and disinfecting before sampling as well as during library preparation, and fewer PCR cycles should be tested to assess the effect of PCR amplification on mapping efficiencies or further data analyses.</p><p>Although contamination reduced actual sequencing depth of each sample, our power to accurately detect DEGs was strong, in part due to the use of four replicates of each group. Increased numbers of biological replicates have been proven to improve power to estimate differential expression significantly regardless of sequencing depth, especially for analysis of highly and moderately expressed genes (Liu et al., 2013; Williams et al., 2014). Also, according to the recommended sequencing depth for human RNA-seq studies (ENCODE Guidelines and Best Practices for RNA-seq, Revised December 2016), thirty million aligned reads/pairs are required for adequate quantification and evaluation of transcripts in long RNA-seq library. A rough calculation relative to human genome and transcriptome sizes suggested that two million aligned reads/pairs are needed for T. californicus RNA-seq analysis, indicating that our sequencing depth (> 3 million average mapped reads per sample) is sufficient. Last but not least, we excluded genes that have less than or equal to 30 total counts across all the samples and 13,228 genes (84.5% of total genes) remained, indicating that our sequencing depth covers the majority of the coding part of the genome.</p><p>Variation among stress treatments accounted for 10.85% of the gene expression variance in our study. Transcriptomic response to the two oxidative stressors, hydrogen peroxide and yeast diet, was fairly distinct, with 19.3% overlap in DEGs for females and 9.5% overlap in males. This modest overlap implies substantial divergence in the physiological pathways responding to prooxidant conditions (H2O2) vs. reduced antioxidant conditions (yeast diet).</p><p>Variation among individuals accounted for 26.85% of the gene expression variance in our study. This is greater than the variance found among stress treatments, underscoring the information that can be missed through the common reliance on biological averaging in pooled samples. Because natural selection works on differences among individuals, understanding among-individual variation is critical to predicting how populations adapt to oxidative stress. A total of 1,215 genes (13.4% of 9,036 genes analyzed here) displayed at least 90% of gene expression variation explained by among-individual variation and GO terms "protein phosphorylation" (GO: 0006468, FDR = 0.045) and "protein kinase activity" (GO: 0004672, FDR = 0.012) were found to be enriched. Protein phosphorylation catalyzed by kinases is a ubiquitous control mechanism by altering protein structure and function to regulate gene transcription and translation, membrane transport and multiple metabolic pathways (Johnson and Barford, 1993). Since it is a flexible type of modification on proteins to respond to environmental conditions and external signals (Johnson and Barford, 1993), this result may also indicate how among-individual variation for stress response was maintained at the transcription level.</p><p>The factor with the greatest effect was sex, accounting for 62.3% of the variance in gene expression. Within treatment groups, at least 89% of the sex-biased genes were female-biased. In contrast, comparisons between control and stress treatment groups showed at least 96% of the DEGs occurring in males. This result is similar to previous work on the same species which showed an overall female-bias in gene expression within treatments, while differential expression under thermal stress, salinity stress and copper stress was greater in males (Foley, 2017). This greater transcriptomic response to stress in males may be due to their lower stress tolerance, which has been reported for a range of stressors including temperature, salinity and marine pollutants (Willett, 2010; Kelly et al., 2012; Foley et al., 2019). Because we only measured gene expression at one time point, it could be that female transcriptomic response to oxidative stress was initially higher, and returned to control levels more rapidly over the 35-day experiment, showing stronger transcriptomic resilience (e.g. Franssen et al., 2011; Brennan et al., 2015; Seneca and Palumbi, 2015). Interestingly, we found that females differentially expressed fewer genes but with a greater magnitude of fold change, indicating a more targeted response to oxidative stress in females.</p><p>Across the three treatments, between 5% and 17% of genes were found to show sex-bias. Sex-biased gene expression has been studied most extensively in taxa with heteromorphic sex chromosomes, where it can occur in up to 90% of the transcriptome (Ingleby et al., 2015). Because sex chromosomes are often predicted to be enriched for sexually antagonistic genes (e.g. Rice, 1984; Ellegren and Parsch, 2007), sex-biased expression might be expected to be reduced in species with non-chromosomal sex determination. Although studies of such taxa are too scarce to test this hypothesis, studies to date in species lacking heteromorphic sex chromosomes have found moderate frequencies of sex-biased gene expression (41% in two bivalve species, Ghiselli et al., 2018; 17–42% in three Daphnia species, Molinier et al., 2018). To our knowledge, our study is the first to report sex-specific gene expression in a species with polygenic sex determination, and the frequency of sex-biased genes was found to be modest. Importantly, our estimate of sex-bias could be inflated since we used whole animals, rather than specific tissues, and males and females can be expected to differ in the proportion of tissue types in addition to the expression of genes within tissue types.</p><p>Functional analysis in this study was largely focused on males, due to the small number of DEGs found in females. In yeast fed males, protein synthesis, ATP synthesis and ATP generation processes were found to be down-regulated. As a major site of energy source, mitochondria use oxygen and sugars to generate ATP along with endogenous ROS during electron transfer through the respiratory chain (Derick et al., 2001; Li et al., 2013). For these yeast fed males with compromised antioxidant defense it may be advantageous to decrease ROS production by slowing down the process of oxidative phosphorylation. Additionally, the excessive ROS is demonstrated to be able to damage ribosomal RNA inside the ribosome, and subsequently inhibit the process of translation (Willi et al., 2018). This was consistent with our findings in the yeast fed male group, especially downregulation of the GO term "structural constituent of ribosome".</p><p>In the H2O2 treated male group, although several known oxidative stress-responding genes were found to be differentially expressed, no directly-related GO terms were identified. One candidate gene was identified to be up-regulated across all four groups but with different expression levels between sexes. According to our sequence analysis, this candidate gene may be copepod-specific. Upon validation with more experiments, this gene may be used as a biomarker for sex-specific oxidative stress detection and observation in copepods.</p><!><p>Individual T. californicus were used to test sex-specific transcriptional responses to oxidative stress induced by pro-oxidant conditions (H2O2 exposure) and reduced antioxidant conditions (yeast diet). Our study is the first to apply single individual RNA-seq to copepods. Expression variation among individuals was found to be even greater than that among stress treatments, illustrating the downside of traditional pooling strategies. This study is also the first one to report sex-specific gene expression in this species, and sex differences accounted for the majority of expression variance found in our study. While gene expression was female-biased within treatments, males differentially expressed a substantially greater number of genes in response to oxidative stress. Female response to oxidative stress was found to be more targeted, differentially expressing a smaller number of genes but with a larger magnitude of fold change. Together, our results underscore the importance of incorporating individual and sex-specific responses into experimental designs aiming at understanding transcriptomic stress responses.</p>