Skip to main content

Sex differences in early transcriptomic responses to oxidative stress in the copepod Tigriopus californicus



Patterns of gene expression can be dramatically different between males and females of the same species, in part due to genes on sex chromosomes. Here we test for sex differences in early transcriptomic response to oxidative stress in a species which lacks heteromorphic sex chromosomes, the copepod Tigriopus californicus.


Male and female individuals were separately exposed to control conditions and pro-oxidant conditions (hydrogen peroxide and paraquat) for periods of 3 hours and 6 hours. Variance partitioning showed the greatest expression variance among individuals, highlighting the important information that can be obscured by the common practice of pooling individuals. Gene expression variance between sexes was greater than that among treatments, showing the profound effect of sex even when males and females share the same genome. Males exhibited a larger response to both pro-oxidants, differentially expressing more than four times as many genes, including up-regulation of more antioxidant genes, heat shock proteins and protease genes. While females differentially expressed fewer genes, the magnitudes of fold change were generally greater, indicating a more targeted response. Although females shared a smaller fraction of differentially expressed genes between stressors and time points, expression patterns of antioxidant and protease genes were more similar between stressors and more GO terms were shared between time points.


Early transcriptomic responses to the pro-oxidants H2O2 and paraquat in copepods revealed substantial variation among individuals and between sexes. The finding of such profound sex differences in oxidative stress response, even in the absence of sex chromosomes, highlights the importance of studying both sexes and the potential for developing sex-specific strategies to promote optimal health and aging in humans.


Oxygen, as a necessity to aerobic organisms, acts as a terminal site for accepting electrons in respiration, and at the same time, works as a generator of dangerous side products – reactive oxygen species or ROS [1]. Aerobic organisms are frequently exposed to various endogenous and environmental conditions which give rise to the accumulation of ROS. Aquatic organisms are faced with a particularly high number of ROS sources, including variation in temperature, salinity and oxygen level, as well as other factors such as heavy metal ions, pesticides, herbicides and oil-related pollutants [2]. Accumulated ROS may damage proteins, lipids, DNA, RNA and other cellular components, thereby resulting in severe physiological dysfunction [3, 4]. Consequently, organisms have developed a series of ROS scavenging systems, including antioxidant enzymes such as superoxide dismutase and catalase and non-enzymatic antioxidants such as glutathione and ascorbic acid [1, 2].

Usually ROS levels are maintained at a dynamic and steady state, with the continuous generation of ROS being balanced by continuous elimination [2]. This steady state concentration is favored since some kinds of ROS (e.g. H2O2) have been found to play important roles in regulation of cellular processes as signaling molecules and also protect against infection [5,6,7]. However, when ROS levels are enhanced, a disturbance of redox status, termed oxidative stress, will be triggered. The antioxidant system will then be provoked to scavenge the excessive ROS to prevent oxidative stress-induced cellular damage, and eventually the ROS level will decrease to a steady-state level, during which the expression of hundreds or thousands of genes may be affected [7]. Oxidative stress not only induces ROS detoxification-related gene expression, such as genes encoding antioxidant enzymes, but also influences the expression of genes involved in signal transduction, protein or lipid metabolism and transcriptional or translational regulation (e.g. [8,9,10,11]).

The effects of oxidative stress have been proven to differ substantially between sexes. Under exposure to ethanol, Drosophila melanogaster females exhibited higher resistance to oxidative stress as evident from lower mortality and ROS levels as well as higher activities of antioxidant enzymes [12]. Similarly, studies in rats have found higher protection from antioxidant enzymes and lower oxidative damage in females [13,14,15]. Furthermore, as reviewed in multiple studies, the sex-specific effects of oxidative stress might play a critical role in sex differences in development of various diseases and lifespan [16, 17]. Sex-biased responses to oxidative stress may be particularly relevant to sex-specific aging, due to the relationship between oxidative stress and aging and the correlation between oxidative stress tolerance and long lifespan [12, 18,19,20,21]. Thus, measuring transcriptional responses to oxidative stress in both sexes will not only help elucidate sex-specific oxidative stress response systems, but also shed light on the molecular mechanisms of sex differences in longevity.

This study is focused on the intertidal copepod Tigriopus californicus, an ideal lab model due to its short generation time, easy culturing and extensive genetic and genomic resources [22]. The species is also an alternative model for sex-specific studies due to the absence of sex chromosomes [22,23,24,25,26], where sex differences will not be complicated by dosage compensation or sexual conflict [27,28,29,30,31]. Instead, sex determination is polygenic with a small environmental component [23,24,25,26]. Despite the absence of sex chromosomes in this species, substantial sex differences have been documented, with males showing lower tolerance to a wide range of stressors [32,33,34]. Our previous study successfully applied single-individual RNA-seq to this species for the first time, in an analysis of sex-specific transcriptomic responses to long-term oxidative stress induced by 35 days of exposure to hydrogen peroxide (H2O2) and a yeast diet [35]. Our results showed that less tolerant male copepods differentially expressed more genes in response to the long-term oxidative stress treatments while females differentially expressed fewer genes but with larger magnitudes of fold change. Here we test for sex differences in response to short-term oxidative stress. While this study did not directly measure oxidative stress (such as by ROS levels), it used two chemicals, hydrogen peroxide and paraquat that are well-established as pro-oxidants. Because transcriptomic response to oxidative stress in other species has shown dramatic temporal changes [11, 36], sex-specific responses in T. californicus in the first few hours of exposure may be very different from those after a month of exposure.

Sex-specific response may also differ substantially between exposures to different pro-oxidants. In Drosophila melanogaster, Pomatto et al. [37] found that females but not males acclimated to H2O2 stress, whereas males but not females acclimated to paraquat (N,N′-dimethyl-4,4′-bipyridinium dichloride) stress. H2O2 is a commonly used stressor for oxidative stress exposure since it is produced by the mitochondrial electron-transport chain and directly increases the ROS level [1, 38]. Paraquat is an herbicide, and its toxicity comes from the generation of superoxide anion which can form more ROS such as hydrogen peroxide and hydroxyl radical [39, 40]. Given previous evidence for contrasting sex-specific responses to H2O2 and paraquat, we assessed sex differences in transcriptomic responses to these same chemicals. While most transcriptomic studies in T. californicus have pooled hundreds of individuals to increase the total amount of RNA [41,42,43,44,45,46], this can greatly increase pooling bias and false discovery rate [47,48,49]. Here we instead used single-individual RNA-seq to assess individual variation in stress response, which was shown to be substantial in our previous work [35]. To assess temporal dynamics of an initial stress response, gene expression was characterized at two different time points (3 h and 6 h).


Single-individual RNA-seq

To assess sex-specific transcriptomic responses to oxidative stress, we performed single-individual RNA-seq on both sexes after exposure to two stressors (H2O2 and paraquat) for 3 h and 6 h. Four biological replicates were used for all treatments at each time point within each sex. A total of 48 individual samples were sequenced on two lanes, yielding ~ 917.7 million 150 bp read pairs. After trimming adapters, low quality bases (base quality < 20) and short reads (trimmed length < 36 bp), ~ 824.8 million read pairs were retained ranging from 9.9 to 27.4 for each individual (Table S1). Clean reads were mapped to annotated genes in the reference genome, resulting in 39.5% mapping efficiency on average (Table S1). The read sequences were deposited at the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) under the accession numbers SRR10444974 – SRR10445021.

Alignment of the non-mapped reads revealed contamination from bacteria, fungi, algae and amplified genomic DNA from other species studied in the lab. This contamination is likely aggravated by the library preparation protocol, which requires twenty cycles of PCR. Despite the contamination, the quality of the mapped sequence data is high. The average number of mapped read pairs per sample is just under 7 million (Table S1), which is higher than our previous transcriptome study and well above the sequencing depth recommended for human RNA-seq studies, after adjusting for genome size (ENCODE Guidelines and Best Practices for RNA-Seq: Revised December 2016). Our power to identify DEGs is further supported by the use of four replicates per treatment. Together, the sequencing depth and replication in this study are more than sufficient for analyses of the moderately or highly expressed protein coding genes [50, 51]. In addition, we filtered out genes with less than or equal to 50 total counts across all the individuals and 13,739 genes (87.8% of total genes) remained, indicating the relatively satisfying coverage of the genome coding regions.

Principal component analysis (PCA) was conducted to show relationships among the 48 samples and identify possible outliers. One female and one male under the control condition at 6 h, and one male under the paraquat condition at 3 h were removed as outliers from further analyses. The resultant PCA (Fig. 1a) showed that 18.0% of the variance among samples is attributed to PC1, and 10.0% to PC2. Males and females were separated by both PC1 and PC2 (meaning that a diagonal line can be drawn to separate the two sexes), regardless of treatments and time points, except for one female under the H2O2 treatment at 3 h. No obvious effect of treatment or time point was observed in the PCA plot. Sample clustering analysis (Fig. 1b) also showed that sex differences were greater than treatment differences with the two sexes forming two separate clusters. Additionally, the two treatment groups were clustered together within each sex, while separate from control groups (Fig. 1b). ANOVA revealed that, on average, at 3 h 37.0, 54.0 and 9.0% of the gene expression variation, and at 6 h 51.1, 42.9 and 6.0% of the gene expression variation were explained by variations between sexes, among individuals and among treatments, respectively.

Fig. 1
figure 1

Expression patterns across sexes, treatments and time points. a Principal component analysis (PCA) of the top 500 genes in terms of variance across individuals. Information on sex and treatment is indicated in the plot. Samples with a black border were collected at the 3 h time point, and samples without a black border were collected at the 6 h time point. The proportion of variance explained by the first two principal components is indicated beside the two axes. b Heatmap of average transformed expression values across all groups. The two time points were combined for each group to show the effects of sex and treatment because they did not show overall apparent differences in gene expression. Expression values were transformed by the vst function in DESeq2. Hierarchical clustering was conducted using one minus Pearson correlation in Morpheus. For each gene, high relative expression is indicated by warmer colors (red), and low relative expression is indicated by colder colors (blue)

Differential expression analysis between sexes within treatments

The two sexes were compared under each treatment at each time point to assess sex-bias in gene expression. Of the 14,230 genes analyzed here, 8.5–24.3% of the genes showed sex-biased expression within treatments (Fig. 2). Within each treatment, 78.7–88.6% of the genes with sex-biased expression are female-biased (For the same gene, females have higher expression levels than males; Fig. 2). Under oxidative stress, less than 37.0% of the sex-biased genes were shared by both treatments (445 at 3 h, 421 at 6 h), indicating that sex-biased gene expression is stressor-specific.

Fig. 2
figure 2

Volcano plots showing the sex-biased gene expression within each group. The number and percentage (relative to all the genes analyzed in this study) of sex-biased genes are indicated in the titles of each plot. Red dots represent highly expressed genes in females (female-biased gene expression) and blue dots represent highly expressed genes in males (male-biased gene expression). The number of both female and male biased genes is indicated with the same color in the plots

Differential expression analysis between treatments within sexes

The analyses focused on comparisons to the control at the same time point to avoid complications due to batch effects of sequencing. To understand sex differences in gene expression responding to short-term oxidative stress, comparisons between treated groups and control groups were performed within each sex at each time point. Sex differences were observed with males differentially expressing more genes (between 81.0 and 91.0% of the total DEGs) than females under each stress treatment at each time point (Fig. 3). Analysis of the percentage of shared DEGs between contrasts at each time point (Fig. 4) showed that the two treatment groups shared a greater proportion of DEGs within sex than that shared by two sexes within treatment. Both findings indicated that the sex difference is a major source of variation in gene expression, consistent with the results from the PCA and the sample clustering analysis. At each time point, the fold changes of genes differentially expressed under either H2O2 or paraquat within each sex were collected for comparisons. At the 3 h time point, the females had average absolute fold change values of 16.9 and 10.8 while the males had means of 14.5 and 11.8 under H2O2 (Mann-Whitney-Wilcoxon test between sexes, p = 4.99e – 05) and paraquat treatments (p = 0.21), respectively. Also, at the 6 h time point, the females displayed mean values of 17.1 and 18.2 while the males averaged 13.5 and 13.2 under H2O2 (p = 2.2e – 16) and paraquat treatments (p = 2.2e – 16), respectively. Thus, under short-term oxidative stress, females differentially expressed fewer genes in all four comparisons but with greater magnitudes of fold change in three out of four comparisons.

Fig. 3
figure 3

Volcano plots showing differentially expressed genes (DEGs) within each sex under each treatment at each time point. a DEGs under H2O2 treatment at 3 h. b DEGs under paraquat treatment at 3 h. c DEGs under H2O2 treatment at 6 h. d DEGs under paraquat treatment at 6 h. The number and percentage (relative to all the genes analyzed in this study) of DEGs are indicated in the plots. Red dots represent DEGs in females and blue dots represent DEGs in males

Fig. 4
figure 4

Venn diagrams showing the number of differentially expressed genes shared between contrasts at 3 h (a), at 6 h (b) and between the two time points (c)

Further, more DEGs were identified at 6 h than those found at 3 h within each sex under each treatment (Fig. 3). In terms of percentages of up-regulation and down-regulation, individuals at 6 h up-regulated more genes than those at 3 h under the same treatment within the same sex (Table 1). Within time points (Fig. 4a and b), females shared a smaller proportion of DEGs between stress treatments than males did, indicating a more stressor-specific response. Between time points (Fig. 4c), females also shared a smaller proportion of DEGs than males did, indicating a more temporally dynamic response.

Table 1 The number and percentage of differentially expressed genes within each comparison

GO enrichment analysis

Up-regulated and down-regulated genes within each contrast were used for GO enrichment analyses separately (Tables S2, S3, S4, S5). The GO enrichment results were further investigated to show enriched GO terms shared between groups (Table S6). Both sexes up-regulated peptide catabolic process (GO:0043171) and peptidase activity (GO:0004175) at 3 h. At 6 h, peptidase activity (GO:0070011) was still up-regulated in males but down-regulated in females. Hydrolase activity (GO:0016787) was found to be up-regulated in females at 3 h and in males at 6 h. Also, both sexes up-regulated synaptic signaling (GO:0098916, GO:0099536, GO:0099537) only at 6 h. In females, pyruvate and lactate transport related GO terms (GO:0006848, GO:1901475, GO:0015727, GO:0035873, GO:0015129) were up-regulated, while renal function (GO:0072127, GO:0072128, GO:0072129, GO:0072130), gonadotropin secretion (GO:0032274) and luteinizing hormone secretion (GO:0032275) were down-regulated at both time points. However, no GO terms were shared by both time points in males, while at the later time point, more oxidative stress responding GO terms, such as autophagy (GO:0010508), protein repair (GO:0030091) and catalytic activity (GO:0003824) were up-regulated.

Identification of oxidative stress responding genes

Oxidative stress has been well studied, and a wide range of genes have been identified to respond to or scavenge ROS. Based on the available gene annotation information for T. californicus [22], we searched our DEG lists for antioxidant genes identified in previous work [52,53,54]. We searched in particular for heat shock proteins, which were known to respond to oxidative stress in Tigriopus (e.g. [55]) and for protease genes, since T. californicus has been shown to exhibit sex differences in proteolytic activity following stress exposure [34]. Table 2 listed 12 gene categories found to be differentially regulated. Notably, catalase, glutathione reductase and thioredoxin reductase were not found to be differentially expressed in any group.

Table 2 The number of differentially expressed oxidative stress responding genes

For the antioxidants and heat shock proteins (Table 2), the vast majority of DEGs were up-regulated. More up-regulated genes were found in male groups than female groups and more in 6 h groups than 3 h groups. No gene up-regulation was found in female groups at 3 h. For protease genes (Tables 2 and S7), at the 3 h time point no proteases were differentially expressed in females, while six and three genes (two genes were shared) were down-regulated in males under H2O2 and paraquat, respectively. However, at the later time point, females up-regulated the same two genes to respond to both stressors, while 11 (2 down-regulated and 9 up-regulated) and 16 (3 down-regulated and 13 up-regulated) genes were differentially expressed in males under H2O2 and paraquat respectively with 9 genes shared between stressors.


Our previous work [35] used single-individual RNA-seq in the copepod T. californicus to assess sex differences in transcriptomic responses after long-term exposure (35 days) to oxidative stressors. The current study took a similar approach to assess sex differences after short-term exposure (3 h and 6 h) to oxidative stressors. Despite the absence of sex chromosomes in this species, both studies found substantial sex differences in gene expression. Both studies showed that sex-biased expression was female-biased within each group, but the proportion of genes with sex biased gene expression (8.5–24.3%) is much lower than the 34.7–42.7% found in another copepod species Lepeophtheirus salmonis that possesses sex chromosomes [56]. Since extractions used whole animals instead of specific tissues, sex-biased expression within treatments may be partially driven by differences in the tissue types (such as gonads) and proportions in the two sexes. For both studies, pairwise comparisons between stress treated groups and control groups showed that males differentially expressed more genes than females. These differences may be associated with the strong stress tolerance in females and the weak tolerance in males, reflected by their responses to a wide range of stressors including temperature, salinity and marine pollutants [32,33,34]. In addition, both studies found that while females have fewer DEGs in response to oxidative stress, the average magnitudes of fold change were generally greater, suggesting a more targeted response in females. These results show that despite the dynamic temporal nature of stress-induced gene expression, consistent sex differences can be found across widely divergent exposure times.

Our finding of a greater number of stress DEGs in the less tolerant type (male copepods) is similar to those found in a number of other systems including corals [57], snails [58] and abalones [59]. Importantly we cannot say whether the lower number of stress DEGs in females, by itself, indicates that females have a more efficient response, or simply experienced a lower level of stress. Our finding of a more targeted response (fewer DEGs but with greater fold change) in the more tolerant type (female copepods) is particularly similar to results of Chen et al. [59], where the more heat-tolerant abalone line differentially expressed fewer genes yet showed greater magnitudes of change for certain key genes including those in the heat shock protein family. Notably, these other studies compared more and less tolerant populations, while our study is the first, to our knowledge, to show such patterns in comparisons between sexes.

Two time points (3 h and 6 h) were used in this study to capture early responses in transcriptomes. The results demonstrated temporal differences with more DEGs and more up-regulation at 6 h than at 3 h. Interestingly, sex differences were also observed in comparisons between the two time points, with females sharing a smaller fraction of DEGs between time points, thereby inducing different gene sets to cope with oxidative stress over time. Moreover, a search of well-known oxidative stress responding genes in each contrast DEG list showed that females at 3 h did not up-regulate any of these known genes, and both males and females up-regulated more of these known genes at 6 h than at 3 h. Results suggested that glutathione S-transferases (GSTs) may be particularly good marker genes for oxidative stress responses since they were highly expressed in both sexes under both stress treatments at 6 h. Notably, GSTs have been identified as important biomarkers for oxidative stress response to H2O2, heavy metals and high temperature in studies of the congeners T. japonicus and T. kingsejongensis [60,61,62], with stress-induced changes in GST levels paralleling changes in ROS levels [60].

While our previous study [35] used one pro-oxidant treatment and one reduced antioxidant treatment, the present study used two pro-oxidants: H2O2 and paraquat. These chemicals were chosen based on a Drosophila study [37], in which females but not males acclimated to H2O2, while males but not females acclimated to paraquat. Our study instead found that in response to both stressors, males differentially expressed more than four times as many genes as females. Further, the two stressors were found to induce distinct responses within each sex, with only 18.5 and 41.7% overlap in DEGs for females and 27.2 and 49.0% overlap in DEGs for males at 3 h and 6 h respectively. In terms of shared genes between sexes for each stressor, 7.2 and 19.5% overlap were observed at 3 h and 6 h, respectively, suggesting divergence in responding mechanisms induced by H2O2 versus paraquat.

Sex difference was found to be a major factor affecting gene expression, accounting for 37.0 and 51.1% of the total expression variance at 3 h and 6 h, respectively. Variation among individuals was also substantial, accounting for 54.0 and 42.9% of the total expression variance at the two time points. And finally, differences among treatments were relatively small, making up only 9.0 and 6.0% of the total variance at the two time points. Expression variance among individuals in this study was notably higher than the 27.0% found in our previous study of long-term exposure to oxidative stress [35]. Since no individuals died during the short-term exposure, we hypothesize that the previous long-term exposure eliminated susceptible individuals and thereby reduced individual variation. In both studies, expression variance among individuals was greater than that among treatments highlighting the information that is obscured through the common practice of combining individuals into pooled samples. Since both studies used outbred populations, the observed expression variance may be due to genetic variation within populations. Such variation is important to document, as it is the fuel for natural selection to promote adaptive change in oxidative stress tolerance. Individual variation in gene expression could also be caused by variation in mating and/or reproductive status, by exposure time differences or by individual variation in handling stress, including the removal of eggs sacs from females. While sources of individual variation cannot be partitioned in the current study, the magnitude of individual variance found argues for future work explicitly designed to partition sources of individual variation.

GO enrichment analyses were conducted for DEGs in all the contrasts. Female responses to oxidative stress at 3 h demonstrated up-regulation of hydrolytic activity and peptidolytic activity dealing with damaged cellular components, but these activities seemed to cease at 6 h. Additionally, pyruvate and lactate transport were active across the two time points in females, protecting cells from oxidative stress by activating cellular defense and antioxidant actions [63, 64]. Some effects on neural system were also observed at 6 h in both sexes, suggested by the up-regulation of synaptic signaling. In contrast to females, males continually up-regulated peptidase activity through both time points, and at the later time point activated hydrolase activity, protein repair and autophagy to fix or clean up dysfunctional components, indicating the less efficient responses to oxidative stress in males.


In summary, early transcriptomic responses to the pro-oxidants H2O2 and paraquat were profoundly different between males and females. The gene expression variance between sexes was even greater than that among treatments for both time points. In response to oxidative stress, the less tolerant male copepods differentially expressed many more genes, including more up-regulation of antioxidants and heat shock proteins at both time points, and more up-regulation of protease genes at the later time point. While females had far fewer DEGs, their average magnitudes of fold change were greater. Between stress treatments, females shared a smaller proportion of DEGs but more similar expression patterns of antioxidant and protease genes. Between time points, females shared more GO terms but a smaller fraction of DEGs. These overall patterns were consistent with a larger early transcriptomic response in the less tolerant male copepods, with sex differences in the time course of the response and the level of specificity between stressors.


Culture maintenance and generation of age-matched copepods

Copepods were collected from supralittoral splash pools in San Diego, CA (32°74′60″N, 117°25′51″W) during summer 2018. They were confirmed to be T. californicus based on habitat (no other copepod species inhabit this extreme environment [65]), morphology and the following transcriptome sequencing data, and were maintained in a 20 °C incubator prior to the start of the experiment. Sea water for culturing was collected from Catalina Island and salinity ranged from 32 to 34 parts per thousand (ppt). Generation of age-matched animals was facilitated by the well-known reproductive biology of this species [66], with year-round reproduction and no evidence of seasonality. Male T. californicus clasp immature females (copepodid stages II-V) prior to the female’s terminal molt, and the female is inseminated and released upon maturation (stage VI) [66]. Females only mate once during their lifetime and produce successive fertilized egg sacs from stored sperm [66]. To generate cohorts of copepods with known age, 90 gravid and healthy females with orange egg sacs were isolated and held at 20 °C overnight to allow eggs to hatch. Females with orange eggs were targeted because the orange coloration is indicative of late-stage development prior to hatching. Upon hatching, larvae were incubated at 20 °C on a 12 h:12 h light:dark cycle and grown into adults after six naupliar stages and five copepodid stages. They were allowed to develop for 4 weeks or until gravid females were observed and then male and female adults (distinguishable by their antennule structures) were separated into two dishes. Notably, a subset of animals is expected to have mated prior to separation of two sexes. All adult animals were transferred to a new dish weekly to remove newly hatched larvae and thereby prevent mixing generations. Throughout, animals were fed weekly using a mixture of powdered Spirulina (Nutrex Hawaii, USA) and ground Tetramin flakes (Tetra, Germany) at a concentration of 0.1 g of each food per 1 L filtered seawater (37 μm).

Oxidative stress treatment and sample collection

Three media were prepared for oxidative stress treatments: (1) 1.2 mM H2O2 treatment: 12.3 μL 30.0% H2O2 solution (EMD Millipore, Germany) was added to 100 mL filtered seawater; (2) 80 ppm (parts per million, 1 ppm = 0.0001%) paraquat treatment: 0.1 g paraquat (Sigma-Aldrich, USA) was added to 1 mL deionized water to make 100,000 ppm solution and then 80 μL was transferred to 100 mL filtered seawater; and (3) control condition: filtered seawater only. The concentrations of H2O2 and paraquat were determined by averaging the mean 24 h LC50 values (the lethal concentration required to kill 50.0% of the population) for both males and females from the San Diego population (unpublished data). Treatments used eight-week old adults from the same dish, with 30 individuals of each sex used for each of the two stress treatments, and 15 individuals of each sex used for the control group. Higher sample sizes were used in the stress treatments to compensate for potential mortality. No food was provided during the acute stress treatment, and cultures were kept in a 20 °C incubator. Previous time-series exposure experiments to H2O2 in T. japonicus showed that expression of the oxidative stress responding biomarker glutathione S-transferase Sigma-class (GST-S) was down-regulated within the first hour and then increased until 6 h or later post exposure [60,61,62]. To capture the early transcriptomic response in T. californicus, two time points (3 h and 6 h) after the start of the experiment were chosen. Four individuals of each sex were sampled from each group. Each individual was dried on a piece of filter paper for one to two seconds and then transferred into a 2 mL nuclease-free tubes with 30–50 1.1 mm diameter zirconia/silica beads (BioSpec Products, USA). The egg sacs of gravid females were removed on the filter paper using dissecting needles and excluded from further sampling. Samples were flash-frozen in liquid nitrogen and stored at − 80 °C until RNA extraction. No individuals died during the exposure experiments.

Transcriptome sequencing

RNA-seq was performed on individual animals to avoid pooling bias (e.g. [47]) and to assess expression variance among individuals. RNA extraction was started with the addition of 300 μL TRIzol Reagent (Ambion, USA) and homogenization by TissueLyser (Qiagen, USA). The Direct-zol RNA MicroPrep Plus kit (Zymo Research, USA) was utilized for further extraction steps including in-column DNase I treatment, and finally 15 μL nuclease-free water was added to elute RNA from the column. RNA was stored at − 80 °C until preparation of RNA sequencing libraries. Forty-eight libraries were constructed according to the LM-Seq protocol [67]. For details of single-individual RNA-seq library preparation, please refer to our previous study [35]. One modification was made in the protocol to use 6-digit index primers (Illumina RNA PCR Index Primers RPI1-RPI24) to replace the 10-digit index primers due to the small number of libraries in each sequencing lane and the common settings for demultiplexing libraries. As in our previous work [35], 20 cycles of PCR were used for library preparation, since preliminary assays with fewer cycles did not generate sufficient products for construction of sequencing libraries. Additional dual size-selection was conducted by AMPure XP beads (Beckman Coulter, USA). Final libraries were quantified using a Qubit dsDNA HS Assay Kit (Invitrogen, USA) and quality-assessed by Bioanalyzer 2100 system (Agilent, USA). Twenty-four libraries from each time point were pooled together based on their molar concentration for sequencing in one lane. All 48 libraries were sequenced by 150 bp paired-end reads in two lanes of the Illumina HiSeq 4000 platform at Fulgent Genetics (CA, USA).

Differential expression analysis

Adapter removal, quality trimming and length trimming were performed using Trimmomatic v0.38 with default parameters [68]. After evaluation by FastQC v0.11.8, trimmed reads were mapped to the San Diego population reference genome with the accession number of GCA_007210705.1 [22] by HISAT2 v2.1.0 [69], and StringTie v1.3.4d [70] was used to estimate the count value for each annotated gene. PCA was conducted to identify a possible outlier if its PC eigenvalue was more than three standard deviations away from the mean value for at least one of the first five PCs. DESeq2 [71] was utilized for differential expression analysis between sexes as well as between oxidative stress treated groups and control groups at each time point within each sex. P-values were adjusted by the Benjamini and Hochberg (BH) procedure [72] to control the false discovery rate. Genes with adjusted P-value < 0.05 and |fold change| > 2 were considered as differentially expressed genes (DEGs). Gene ontology (GO) enrichment analyses of differentially up-regulated and down-regulated genes was conducted by topGO [73] in R.

Partitioning of gene expression variation

Genes with half of the individuals having less than 1 read per million mapped reads were filtered out, and then normalized expression values for the remaining genes were log2 transformed. A total of 6511 and 6098 genes were kept for the 3 h and 6 h time points, respectively. Variance in gene expression was partitioned by the R package variancePartition version 1.12.1 [74] to estimate the variance explained by effects including sex, treatment and individual replicate at each time point. Detailed methods were described in our previous study [35]. Briefly, the nested analysis of variance (ANOVA) model for each gene was assumed as:

$$ {\mathrm{Y}}_{\mathrm{i}\mathrm{j}\mathrm{k}}=\upmu +{\mathrm{T}}_{\mathrm{i}}+{\mathrm{S}}_{\mathrm{j}\left(\mathrm{i}\right)}+{\mathrm{R}}_{\mathrm{k}\left(\mathrm{i}\mathrm{j}\right)}+{\mathrm{e}}_{\mathrm{i}\mathrm{j}\mathrm{k}} $$

where μ represents the mean expression level for the given gene, and Yijk is the normalized expression level in the kth replicate (R) which is nested within jth sex (S), which is nested within ith treatment (T). 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 due to all three variables: σ2 / (σT2 + σS2 + σR2). An average percentage of variance due to each effect was calculated for each time point.

Availability of data and materials

The sequences generated during this study were deposited at the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) under the accession numbers SRR10444974 – SRR10445021. The reference genome used in our analyses is available with the GenBank assembly accession of GCA_007210705.1, and the genome assembly version is Tcal_SD_v2.1.



Reactive Oxygen Species


Glutathione S-Transferase


Ligation Mediated RNA sequencing


Principal Component Analysis


Differentially Expressed Gene


Gene Ontology


  1. Halliwell B, Gutteridge JM. Free radicals in biology and medicine. USA: Oxford University Press; 2015.

    Book  Google Scholar 

  2. Lushchak VI. Environmentally induced oxidative stress in aquatic animals. Aquat Toxicol. 2011;101(1):13–30.

    Article  CAS  Google Scholar 

  3. Lesser MP. Oxidative stress in marine environments: biochemistry and physiological ecology. Annu Rev Physiol. 2006;68:253–78.

    Article  CAS  Google Scholar 

  4. Kong Q. Lin C-lG. Oxidative damage to RNA: mechanisms, consequences, and diseases. Cell Mol Life Sci. 2010;67(11):1817–29.

    Article  CAS  Google Scholar 

  5. Giorgio M, Trinei M, Migliaccio E, Pelicci PG. Hydrogen peroxide: a metabolic by-product or a common mediator of ageing signals? Nat Rev Mol Cell Bio. 2007;8(9):722.

    Article  CAS  Google Scholar 

  6. Bigarella CL, Liang R, Ghaffari S. Stem cells and the impact of ROS signaling. Development. 2014;141(22):4206–18.

    Article  CAS  Google Scholar 

  7. Paiva CN, Bozza MT. Are reactive oxygen species always detrimental to pathogens? Antioxid Redox Signal. 2014;20(6):1000–37.

    Article  CAS  Google Scholar 

  8. Morel Y, Barouki R. Repression of gene expression by oxidative stress. Biochem J. 1999;342(3):481–96.

    Article  CAS  Google Scholar 

  9. Mostertz J, Scharf C, Hecker M, Homuth G. Transcriptome and proteome analysis of Bacillus subtilis gene expression in response to superoxide and peroxide stress. Microbiology. 2004;150(2):497–512.

    Article  CAS  Google Scholar 

  10. Vanderauwera S, Zimmermann P, Rombauts S, Vandenabeele S, Langebartels C, Gruissem W, Inzé D, Van Breusegem F. Genome-wide analysis of hydrogen peroxide-regulated gene expression in Arabidopsis reveals a high light-induced transcriptional cluster involved in anthocyanin biosynthesis. Plant Physiol. 2005;139(2):806–21.

    Article  CAS  Google Scholar 

  11. Vandenbroucke K, Robbens S, Vandepoele K, Inzé D, Van de Peer Y, Van Breusegem F. Hydrogen peroxide–induced gene expression across kingdoms: a comparative analysis. Mol Biol Evol. 2008;25(3):507–16.

    Article  CAS  Google Scholar 

  12. Niveditha S, Deepashree S, Ramesh S, Shivanandappa T. Sex differences in oxidative stress resistance in relation to longevity in Drosophila melanogaster. J Comp Physiol B. 2017;187(7):899–909.

    Article  CAS  Google Scholar 

  13. Borrás C, Sastre J, García-Sala D, Lloret A, Pallardó FV, Viña J. Mitochondria from females exhibit higher antioxidant gene expression and lower oxidative damage than males. Free Radic Biol Med. 2003;34(5):546–52.

    Article  CAS  Google Scholar 

  14. Kayali R, Çakatay U, Tekeli F. Male rats exhibit higher oxidative protein damage than females of the same chronological age. Mech Ageing Dev. 2007;128(5–6):365–9.

    Article  CAS  Google Scholar 

  15. Gómez-Pérez Y, Gianotti M, Lladó I, Proenza AM. Sex-dependent effects of high-fat-diet feeding on rat pancreas oxidative stress. Pancreas. 2011;40(5):682–8.

    Article  CAS  Google Scholar 

  16. Kander MC, Cui Y, Liu Z. Gender difference in oxidative stress: a new look at the mechanisms for cardiovascular diseases. J Cell Mol Med. 2017;21(5):1024–32.

    Article  Google Scholar 

  17. Tower J. Sex-specific gene expression and life span regulation. Trends Endocrinol Metab. 2017;28(10):735–47.

    Article  CAS  Google Scholar 

  18. Finkel T, Holbrook NJ. Oxidants, oxidative stress and the biology of ageing. Nature. 2000;408(6809):239–47.

    Article  CAS  Google Scholar 

  19. Dröge W. Oxidative stress and aging. In: Roach RC, Wagner PD, Hackett PH, editors. Hypoxia: Through the Lifecycle. Springer Science & Business Media; 2003. p. 191–200.

  20. Vermeulen C, Van De Zande L, Bijlsma R. Resistance to oxidative stress induced by paraquat correlates well with both decreased and increased lifespan in Drosophila melanogaster. Biogerontology. 2005;6(6):387–95.

    Article  CAS  Google Scholar 

  21. Neki N. Oxidative stress and aging. Bangladesh J Med Sci. 2015;14(3):221–7.

    Article  Google Scholar 

  22. Barreto FS, Watson ET, Lima TG, Willett CS, Edmands S, Li W, Burton RS. Genomic signatures of mitonuclear coevolution across populations of Tigriopus californicus. Nat Ecol Evol. 2018;2(8):1250–7.

    Article  Google Scholar 

  23. Ar-rushdi AH. The cytology of achiasmatic meiosis in the female Tigriopus (Copepoda). Chromosoma. 1962;13(5):526–39.

    Article  Google Scholar 

  24. Voordouw MJ, Anholt BR. Heritability of sex tendency in a harpacticoid copepod, Tigriopus californicus. Evolution. 2002;56(9):1754–63.

    Article  Google Scholar 

  25. Alexander H, Richardson J, Anholt B. Multigenerational response to artificial selection for biased clutch sex ratios in Tigriopus californicus populations. J Evol Biol. 2014;27(9):1921–9.

    Article  CAS  Google Scholar 

  26. Alexander H, Richardson J, Edmands S, Anholt B. Sex without sex chromosomes: genetic architecture of multiple loci independently segregating to determine sex ratios in the copepod Tigriopus californicus. J Evol Biol. 2015;28(12):2196–207.

    Article  CAS  Google Scholar 

  27. Isensee J, Witt H, Pregla R, Hetzer R, Regitz-Zagrosek V, Noppinger PR. Sexually dimorphic gene expression in the heart of mice and men. J Mol Med. 2008;86(1):61–74.

    Article  Google Scholar 

  28. Spigler RB, Lewers KS, Ashman TL. Genetic architecture of sexual dimorphism in a subdioecious plant with a proto-sex chromosome. Evolution. 2011;65(4):1114–26.

    Article  Google Scholar 

  29. Disteche CM. Dosage compensation of the sex chromosomes. Annu Rev Genet. 2012;46:537–60.

    Article  CAS  Google Scholar 

  30. Hill GE. Sex linkage of nuclear-encoded mitochondrial genes. Heredity (Edinb). 2013;112:469–70.

    Article  CAS  Google Scholar 

  31. Dean R, Zimmer F, Mank JE. The potential role of sexual conflict and sexual selection in shaping the genomic distribution of Mito-nuclear genes. Genome Biol Evol. 2014;6(5):1096–104.

    Article  CAS  Google Scholar 

  32. Willett CS. Potential fitness trade-offs for thermal tolerance in the intertidal copepod Tigriopus californicus. Evolution. 2010;64(9):2521–34.

    Article  Google Scholar 

  33. Kelly MW, Sanford E, Grosberg RK. Limited potential for adaptation to climate change in a broadly distributed marine crustacean. Proc R Soc B. 2012;279(1727):349–56.

    Article  Google Scholar 

  34. Foley HB, Sun PY, Ramirez R, So BK, Venkataraman YR, Nixon EN, Davies KJA, Edmands S. Sex-specific stress tolerance, proteolysis, and lifespan in the invertebrate Tigriopus californicus. Exp Gerontol. 2019;119:146–56.

    Article  CAS  Google Scholar 

  35. Li N, Arief N, Edmands S. Effects of oxidative stress on sex-specific gene expression in the copepod Tigriopus californicus revealed by single individual RNA-seq. Comp Biochem Physiol D. 2019;31:100608.

  36. Lucau-Danila A, Lelandais G, Kozovska Z, Tanty V, Delaveau T, Devaux F, Jacq C. Early expression of yeast genes affected by chemical stress. Mol Cell Biol. 2005;25(5):1860–8.

    Article  CAS  Google Scholar 

  37. Pomatto LC, Carney C, Shen B, Wong S, Halaszynski K, Salomon MP, Davies KJ, Tower J. The mitochondrial Lon protease is required for age-specific and sex-specific adaptation to oxidative stress. Curr Biol. 2017;27(1):1–15.

    Article  CAS  Google Scholar 

  38. Abele-Oeschger D, Sartoris FJ, Pörtner HO. Hydrogen peroxide causes a decrease in aerobic metabolic rate and in intracellular pH in the shrimp Crangon crangon. Comp Biochem Physiol C. 1997;117(2):123–9.

    Google Scholar 

  39. Suntres ZE. Role of antioxidants in paraquat toxicity. Toxicology. 2002;180(1):65–77.

    Article  CAS  Google Scholar 

  40. Pomeroy A. Biochemical mechanisms of paraquat toxicity. London: Academic; 2012.

  41. Graham AM, Barreto FS. Loss of the HIF pathway in a widely distributed intertidal crustacean, the copepod Tigriopus californicus. Proc Natl Acad Sci U S A. 2019;116(26):12913–8.

    Article  CAS  Google Scholar 

  42. DeBiasse MB, Kawji Y, Kelly MW. Phenotypic and transcriptomic responses to salinity stress across genetically and geographically divergent Tigriopus californicus populations. Mol Ecol. 2018;27(7):1621–32.

    Article  CAS  Google Scholar 

  43. Lima TG, Willett CS. Locally adapted populations of a copepod can evolve different gene expression patterns under the same environmental pressures. Ecol Evol. 2017;7(12):4312–25.

    Article  Google Scholar 

  44. Pereira RJ, Barreto FS, Pierce NT, Carneiro M, Burton RS. Transcriptome-wide patterns of divergence during allopatric evolution. Mol Ecol. 2016;25(7):1478–93.

    Article  CAS  Google Scholar 

  45. Schoville SD, Barreto FS, Moy GW, Wolff A, Burton RS. Investigating the molecular basis of local adaptation to thermal stress: population differences in gene expression across the transcriptome of the copepod Tigriopus californicus. BMC Evol Biol. 2012;12(1):170.

    Article  CAS  Google Scholar 

  46. Barreto FS, Moy GW, Burton RS. Interpopulation patterns of divergence and selection across the transcriptome of the copepod Tigriopus californicus. Mol Ecol. 2011;20(3):560–72.

    Article  Google Scholar 

  47. Rajkumar AP, Qvist P, Lazarus R, Lescai F, Ju J, Nyegaard M, Mors O, Børglum AD, Li Q, Christensen JH. Experimental validation of methods for differential gene expression analysis and sample pooling in RNA-seq. BMC Genomics. 2015;16(1):548.

    Article  CAS  Google Scholar 

  48. Mary-Huard T, Daudin JJ, Baccini M, Biggeri A, Bar-Hen A. Biases induced by pooling samples in microarray experiments. Bioinformatics. 2007;23(13):i313–8.

    Article  CAS  Google Scholar 

  49. Kendziorski C, Irizarry R, Chen KS, Haag J, Gould M. On the utility of pooling biological samples in microarray experiments. Proc Natl Acad Sci U S A. 2005;102(12):4252–7.

    Article  CAS  Google Scholar 

  50. Liu Y, Zhou J, White KP. RNA-seq differential expression studies: more sequence or more replication? Bioinformatics. 2013;30(3):301–4.

    Article  CAS  Google Scholar 

  51. Williams AG, Thomas S, Wyman SK, Holloway AK. RNA-seq data: challenges in and recommendations for experimental design and analysis. Curr Protoc Hum Genet. 2014;83(1):11–3.

    Google Scholar 

  52. Scandalios J. Oxidative stress: molecular perception and transduction of signals triggering antioxidant gene defenses. Braz J Med Biol Res. 2005;38(7):995–1014.

    Article  CAS  Google Scholar 

  53. Birben E, Sahiner UM, Sackesen C, Erzurum S, Kalayci O. Oxidative stress and antioxidant defense. World Allergy Organ J. 2012;5(1):9.

    Article  CAS  Google Scholar 

  54. Kurutas EB. The importance of antioxidants which play the role in cellular response against oxidative/nitrosative stress: current state. Nutr J. 2015;15(1):71.

    Article  CAS  Google Scholar 

  55. Kim BM, Rhee JS, Jeong CB, Seo JS, Park GS, Lee YM, Lee JS. Heavy metals induce oxidative stress and trigger oxidative stress-mediated heat shock protein (hsp) modulation in the intertidal copepod Tigriopus japonicus. Comp Biochem Physiol C. 2014;166:65–74.

    CAS  Google Scholar 

  56. Poley JD, Sutherland BJ, Jones SR, Koop BF, Fast MD. Sex-biased gene expression and sequence conservation in Atlantic and Pacific salmon lice (Lepeophtheirus salmonis). BMC Genomics. 2016;17(1):483.

    Article  CAS  Google Scholar 

  57. Barshis DJ, Ladner JT, Oliver TA, Seneca FO, Traylor-Knowles N, Palumbi SR. Genomic basis for coral resilience to climate change. Proc Natl Acad Sci U S A. 2013;110(4):1387–92.

    Article  CAS  Google Scholar 

  58. Gleason LU, Burton RS. RNA-seq reveals regional differences in transcriptome response to heat stress in the marine snail Chlorostoma funebralis. Mol Ecol. 2015;24(3):610–27.

    Article  CAS  Google Scholar 

  59. Chen N, Huang Z, Lu C, Shen Y, Luo X, Ke C, You W. Different Transcriptomic responses to thermal stress in heat-tolerant and heat-sensitive Pacific abalones indicated by cardiac performance. Front Physiol. 2019;9:1895.

    Article  Google Scholar 

  60. Han J, Jeong CB, Byeon E, Lee JS. Effects of temperature changes on the generation of reactive oxygen species and the expression and activity of glutathione-S transferases in two congeneric copepods Tigriopus japonicus and Tigriopus kingsejongensis. Fisheries Sci. 2018;84(5):815–23.

    Article  CAS  Google Scholar 

  61. Lee KW, Raisuddin S, Rhee JS, Hwang DS, Yu IT, Lee YM, Park HG, Lee JS. Expression of glutathione S-transferase (GST) genes in the marine copepod Tigriopus japonicus exposed to trace metals. Aquat Toxicol. 2008;89(3):158–66.

    Article  CAS  Google Scholar 

  62. Lee YM, Lee KW, Park H, Park HG, Raisuddin S, Ahn IY, Lee JS. Sequence, biochemical characteristics and expression of a novel sigma-class of glutathione S-transferase from the intertidal copepod, Tigriopus japonicus with a possible role in antioxidant defense. Chemosphere. 2007;69(6):893–902.

    Article  CAS  Google Scholar 

  63. Wang X, Perez E, Liu R, Yan LJ, Mallet RT, Yang SH. Pyruvate protects mitochondria from oxidative stress in human neuroblastoma SK-N-SH cells. Brain Res. 2007;1132:1–9.

    Article  CAS  Google Scholar 

  64. Tauffenberger A, Fiumelli H, Almustafa S, Magistretti PJ. Lactate and pyruvate promote oxidative stress resistance through hormetic ROS signaling. Cell Death Dis. 2019;10(9):1–16.

    Article  CAS  Google Scholar 

  65. Dethier MN. Tidepools as refuges: predation and the limits of the harpacticoid copepod Tigriopus californicus (baker). J Exp Mar Biol Ecol. 1980;42(2):99–111.

    Article  Google Scholar 

  66. Burton RS. Mating system of the intertidal copepod Tigriopus californicus. Mar Biol. 1985;86(3):247–52.

    Article  Google Scholar 

  67. Hou Z, Jiang P, Swanson SA, Elwell AL, Nguyen BKS, Bolin JM, Stewart R, Thomson JA. A cost-effective RNA sequencing protocol for large-scale gene expression studies. Sci Rep. 2015;5:9570.

    Article  CAS  Google Scholar 

  68. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30(15):2114–20.

  69. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12(4):357–60.

    Article  CAS  Google Scholar 

  70. Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015;33(3):290–5.

    Article  CAS  Google Scholar 

  71. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.

    Article  CAS  Google Scholar 

  72. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Series B Stat Methodol. 1995;57(1):289–300.

    Google Scholar 

  73. Alexa A, Rahnenfuhrer J. topGO: enrichment analysis for gene ontology. R Package Version. 2010;2:2010.

    Google Scholar 

  74. Hoffman GE, Schadt EE. variancePartition: interpreting drivers of variation in complex gene expression studies. BMC Bioinformatics. 2016;17(1):483.

    Article  Google Scholar 

Download references


Computation for the work described in this paper was supported by the University of Southern California’s Center for High-Performance Computing (


This work was supported by the National Institute on Aging of the U.S. National Institutes of Health [grant R21AG055873 to SE] and the U.S. National Science Foundation [grant DEB-1656048 to SE]. EJH was supported by a University Research Associates Program Fellowship from the University of Southern California.

Author information

Authors and Affiliations



NL and SE conceived and designed experiments. BAF, MP and EJH performed experiments. NL analyzed data and prepared the draft manuscript. BAF and SE revised the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Ning Li.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

All authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Li, N., Flanagan, B.A., Partridge, M. et al. Sex differences in early transcriptomic responses to oxidative stress in the copepod Tigriopus californicus. BMC Genomics 21, 759 (2020).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI:


  • RNA-seq
  • Hydrogen peroxide
  • Paraquat
  • Tolerance
  • Antioxidant