Although the existence of a cellular heat shock response is nearly universal, its relationship to organismal thermal tolerance is not completely understood. Many of the genes involved are known to be regulated by the highly conserved heat shock transcription factor-1 (HSF-1), yet the regulatory network is not fully characterized. Here, we investigated the role of HSF-1 in gene expression following thermal stress using knockdown of HSF-1 by RNA interference in the intertidal copepod Tigriopus californicus. We observed some evidence for decreased transcription of heat shock protein genes following knockdown, supporting the widely acknowledged role of HSF-1 in the heat shock response. However, the majority of differentially expressed genes between the control and HSF-1 knockdown groups were upregulated, suggesting that HSF-1 normally functions to repress their expression. Differential expression observed in genes related to chitin and cuticle formation lends support to previous findings that these processes are highly regulated following heat stress. We performed a genome scan and identified a set of 396 genes associated with canonical heat shock elements. RNA-seq data did not find those genes to be more highly represented in our HSF-1 knockdown treatment, indicating that requirements for binding and interaction of HSF-1 with a given gene are not simply predicted by the presence of HSF-1 binding sites. Further study of the pathways implicated by these results and future comparisons among populations of T. californicus may help us understand the role and importance of HSF-1 in the heat shock response and, more broadly, in organismal thermal tolerance.

As the climate warms and organisms are pushed to their thermal limits, understanding the mechanistic basis underlying thermal tolerance is of increasing interest. Thermal stress elicits a near-universal cellular ‘heat shock response’, which has been characterized in a diversity of organisms over the past 40 years, ranging from slime molds (Loomis and Wheeler, 1980) to Drosophila (Ashburner and Bonner, 1979; Hoffmann et al., 2003) to humans (Kregel, 2002). Many of the genes involved are conserved across taxa (Lindquist, 1986; Feder and Hofmann, 1999), including several families of heat shock proteins (HSPs) involved in processing cellular proteins that have lost functional conformation owing to elevated temperature; HSPs prevent aggregation of denatured proteins, re-fold proteins in correct conformations, or target misfolded proteins for degradation (Morimoto, 1998; Feder and Hofmann, 1999; Sørensen et al., 2003). However, despite generally contributing to thermal tolerance, there is some variation in the expression of the heat shock response in different systems (Morimoto, 1998; Feder and Hofmann, 1999). For example, some systems (e.g. corals and marine snails) show a ‘frontloading’ of heat shock proteins, whereby thermally tolerant organisms demonstrate constitutively elevated expression of HSPs and other genes involved in the heat shock response (Barshis et al., 2013; Gleason and Burton, 2015). In contrast, many other organisms tend to show a more pronounced upregulation of these genes in the face of stress (e.g. mussels in Lockwood et al., 2010; Tigriopuscalifornicus copepods in Schoville et al., 2012; Drosophila in Sørensen et al., 2005).

Differentially expressed genes involved in the heat shock response are often regulated by a highly conserved heat shock transcription factor-1 (HSF-1) (Wu, 1995; Morimoto, 1998; Westerheide et al., 2012). In addition to its role in thermal response, HSF-1 has been shown to be important in a number of other processes across taxa, including development (Jedlicka et al., 1997; Xiao et al., 1999; Anckar and Sistonen, 2007; Åkerfelt et al., 2010), longevity (Steinkraus et al., 2008; Åkerfelt et al., 2010; Seo et al., 2013) and a range of stress responses (Sorger, 1991; Jedlicka et al., 1997; Zhong et al., 1998; Walker et al., 2003; Åkerfelt et al., 2010; Morton and Lamitina, 2013; Furuhashi and Sakamoto, 2015; Brunquell et al., 2016; Labbadia et al., 2017; Takii et al., 2017; Mazin et al., 2018). Following heat shock, HSF-1 monomers trimerize and accumulate in the nucleus (Westwoodt and Wu, 1993; Morimoto, 1998), where they bind to highly conserved heat shock elements (HSEs), which are short inverted repeat sequences in the promoter region of genes. HSF-1 binds the canonical HSE sequence motif; substitutions within HSEs can decrease expression of the associated gene (Sorger, 1991; Fernandes et al., 1994; Åkerfelt et al., 2010; Tangwancharoen et al., 2018).

An important aspect of cellular response to thermal stress is the variety of genes involved. Though HSPs are clearly important, there can be anywhere from hundreds to thousands of genes with differential expression following a given thermal stress, and many of them are not HSPs (e.g. Brunquell et al., 2016; Lockwood et al., 2010; Schoville et al., 2012; Stillman and Tagmount, 2009). For example, several studies have found cuticle protein genes or collagen-associated genes to be involved in the heat shock response; Brunquell et al. (2016) found cuticle genes to be the top category of genes upregulated by HSF-1 independently of heat stress, while Schoville et al. (2012) found that cuticle genes were upregulated or downregulated by populations of T. californicus that differed in their heat tolerance.

Though the mechanism of HSF-1 action seems straightforward, previous research has found that having a ‘perfect’, canonical HSE does not guarantee binding by HSF (Guertin and Lis, 2010), so further investigation into HSF binding and gene activation in vivo is necessary to understand the role of this transcription factor in the heat shock response. In the present study, we used T. californicus, an intertidal copepod with a well-characterized heat shock response (e.g. Pereira et al., 2017; Tangwancharoen et al., 2018). Previous studies have shown that T. californicus upregulates HSPs in response to stress, though the degree to which upregulation occurs depends on the thermal tolerance of the population being studied. Although there are several HSF family members in many organisms, only HSF-1, the most conserved and active HSF in the heat shock response (Sarge et al., 1993; Pirkkala et al., 2001; Brunquell et al., 2016), has been found in the T. californicus draft genome (Barreto et al., 2018). Thus, understanding the heat shock response in different populations necessitates first understanding its regulation by HSF-1. Gene knockdown via RNA interference (RNAi) has previously been successfully employed to knock down expression of a specific small HSP in T. californicus (Barreto et al., 2015); here, we used this approach for knockdown of HSF. Following validated knockdown and subsequent heat stress, we used RNA sequencing (RNA-seq) to identify genes that are differentially expressed following heat stress with and without knockdown of HSF.

Copepod collection and maintenance

Tigriopus californicus (Baker 1912) copepods were collected from tidepools in San Diego, CA, USA (32°45′N, 117°15′W), and maintained in 400-ml beakers containing 200 ml of 0.2 µm-filtered seawater. Beakers were kept at 20°C with a 12 h:12 h light:dark cycle. Copepods were fed powdered spirulina ad libitum and maintained in laboratory conditions for at least 1 month prior to use in this study.

dsRNA synthesis

Because T. californicus HSF-1 lacks introns (Tangwancharoen et al., 2018), genomic DNA was used as template for synthesis of dsRNA. Genomic DNA was extracted by incubating individual copepods in 15 µl lysis buffer (Willett and Burton, 2001) at 55°C for 90 min followed by 95°C for 15 min. Following lysis, PCR was used to amplify 300–600 bp regions of the HSF gene using primers with T7 promoter sequences (TAATACGACTCACTATAGGG) added on the 5′ end (see Table 1 for primer sequences). PCR reactions were performed in 25-µl reactions and contained 1X GoTaq Green Master Mix (Promega), 0.8 µmol l−1 primers and 1 µl DNA extract. Thermal cycling parameters were 95°C for 3 min, 34 cycles of 95°C for 20 s, 58°C for 1 min and 72°C for 1 min, and finally 72°C for 5 min. The molecular weight of the PCR product was checked on a 1.5% agarose gel.

Table 1.

Primer sequences for dsRNA generation (excluding T7 tails) and qPCR

Primer sequences for dsRNA generation (excluding T7 tails) and qPCR
Primer sequences for dsRNA generation (excluding T7 tails) and qPCR

Following DNA amplification, PCR product was purified using Sephadex G-50 columns (GE Healthcare) and 1 µg product was used as a template for transcription in the HiScribe T7 High Yield RNA Synthesis Kit (New England Bioscience) following the manufacturer's protocol, which included overnight incubation for dsRNA <300 bp and 2 h incubation for dsRNA >300 bp. Transcription product was treated with Turbo DNase (37°C, 1 h; Thermo Fisher Scientific) to remove DNA template and purified with an RNeasy column (Qiagen). Following the method of Barreto et al. (2015), we also produced a control dsRNA by amplifying a 392-bp region of the pCR4-TOPO plasmid (Invitrogen).

Electroporation

For dsRNA delivery, groups of 100 adult copepods were placed into 0.4-cm gap electroporation cuvettes (Bio-Rad) containing 140 µl artificial seawater and 60 µl target dsRNA (50 µg), control dsRNA or DEPC-H2O. Barreto et al. (2015) did not see significant differences between control dsRNA and DEPC-H2O electroporation. Here, we used both methods and also did not see differences in control gene expression. Copepods were electroporated in a Gene Pulser Xcell (Bio-Rad) with square-wave pulses (1–100 V and 10–50 V) and subsequently transferred back to Petri dishes in a 20°C incubator and held for 2 days before use in any further experiments.

Heat shock

Copepods from either control dsRNA or target dsRNA groups that survived electroporation were counted and split into two groups, with half serving as control animals (non-heat stressed) and half placed into 15 ml conical tubes containing 4 ml seawater and heat stressed by submerging the tubes in a 35°C water bath for 1 h. After this treatment, they were transferred back to 20°C incubators. Subsets of each group were either (1) left at 20°C for survivorship assessment 3 days later (10 animals) or (2) allowed to recover for 1 h before homogenization for western blot (five animals) or RNA extraction (remaining animals, at least 10).

RNA extraction and gene expression validation by qPCR

RNA was extracted from groups of at least 10 animals (depending on the number remaining after electroporation) using TRI Reagent (Sigma-Aldrich) in accordance with the manufacturer's instructions. cDNA was synthesized using the High Capacity RNA-to-cDNA kit (Applied Biosystems) and diluted to 3 ng µl−1. Quantitative PCR was performed in 15 µl reactions containing 9 ng template, 1X iTaq Universal SYBR Green Supermix (Bio-Rad) and 0.35 µmol l−1 of each primer (Table 1). Reactions were run on a Stratagene MX3000P (Agilent) with 95°C denaturation for 2 min, 40 cycles of 95°C for 10 s and 59°C for 20 s. Gene expression was assessed using the 2−ΔΔCt method (Schmittgen and Livak, 2008) and normalized with the GAPDH gene.

Western blot validation

Groups of five copepods were homogenized in a mixture of 20 µl 10% sodium dodecyl sulfate and 20 µl Laemmli sample buffer using polypropylene pestles (Thomas Scientific) in 1.5 ml microcentrifuge tubes. The homogenate was heated at 95°C for 5 min and then frozen at −20°C prior to use. Samples were run on a 7.5% SDS-PAGE gel, probed with an anti-HSF polyclonal antibody. For antibody production, T. californicus HSF-1 was amplified with an N-terminal 6xHis-tag, inserted into a pProEx Htb expression vector (Invitrogen) and transformed into Escherichia coli BL21 (DE3) pLysE cells (see Tangwancharoen et al., 2018). Expression of recombinant HSF-1 was induced using addition of isopropyl-β-D-1-thiogalactopyranoside (IPTG) and purified using His60 Ni Superflow Resin columns (Clontech). Purified HSF-1 was sent to GenScript (Piscataway, NJ, USA) for antibody production in rabbit. Purified antibody was prepared by GenScript. PAGE gels were electrophoretically transferred to nitrocellulose membranes and probed with the primary (HSF-1) antibody; blots were then probed with a HRP-conjugated goat anti-rabbit secondary antibody and incubated with SuperSignal West Pico Chemiluminescent Substrate (Thermo Fisher Scientific). Western blots were visualized on a ChemiDoc (Bio-Rad) and analyzed using Image Lab 4.0.1 (Bio-Rad). Analysis was performed by measuring the relative intensity of all bands on each blot as a percentage, with the darkest band being set to 100%. Relative intensity of the replicates used in this study is therefore reported as the average of the intensities of the three control (control RNAi) and three treatment (HSF-1 RNAi) bands on each of three trial gels, as well as the relative intensities of the bands corresponding to each of the three treatment samples chosen for sequencing (Table 2). See Fig. S1 for antibody validation gel images.

Table 2.

Relative band intensity measured from a western blot of three trials comparing control RNAi and HSF-1 RNAi samples (n=3 for each trial)

Relative band intensity measured from a western blot of three trials comparing control RNAi and HSF-1 RNAi samples (n=3 for each trial)
Relative band intensity measured from a western blot of three trials comparing control RNAi and HSF-1 RNAi samples (n=3 for each trial)

RNA sequencing

After confirming a potentially successful knockdown of HSF using qPCR results and western blot visualization, RNA samples were obtained from four groups: control with heat stress, control with no heat stress, knockdown with heat stress and knockdown with no heat stress (three biological replicates each). RNA remaining after use in cDNA synthesis/qPCR reactions (described above) was quantified and sent to Novogene (Davis, CA, USA) for sequencing on a 150 bp Illumina HiSeq paired-end platform. Novogene performed quality control on the RNA samples (Agilent 2100 analysis) and prepared the RNA-seq libraries using New England Biosciences Next® Ultra™ RNA Library Prep Kit.

Differential expression analysis

Sequences were trimmed for quality using CLC Genomics Workbench version 11.0.1 (Qiagen Bioinformatics). Trimming involved removal of low-quality sequences (limit=0.01), removal of ambiguous nucleotides (maximal allowed=2), removal of Novogene adapter sequences, and finally removal of sequences less than 50 nucleotides in length. Reads were then mapped to the T. californicus transcriptome (version 4; https://i5k.nal.usda.gov/Tigriopus_californicus) using CLC Genomics Workbench (minimum length fraction of read overlap=0.5, minimum sequence similarity=0.8). Uniquely mapped reads were used in the subsequent analysis. Differential gene expression analyses were performed using DESeq2 (Love et al., 2014). Read counts were normalized and filtered by DESeq2 using the average expression of each gene to establish a filtering threshold (Love et al., 2014). The false discovery rate (FDR) for differential expression was set to Padj<0.05 and ±LFC (logarithmic fold change) of 1.

Principal component analysis

The principal component analysis (PCA) was plotted using the plotPCA function of the DESeq2 package using the top 1000 most differentially expressed genes.

Heat map

Heat maps were made in R (https://cran.r-project.org) using hierarchical clustering on dissimilarities of all genes showing significant differential expression in either the control versus knockdown comparison or the heat stress versus non-heat stress comparison.

Heat shock element genome scan

In order to find genes with upstream HSEs, we scanned the T. californicus genome (Barreto et al., 2018) using the ChIPpeakAnno package in R (Zhu et al., 2010; Zhu, 2013). Parameters were set to identify genes with canonical HSEs (‘NGAANNTTCNNGAAN’) within 2000 bp upstream of their start site. Blast2Go was used to identify gene functional categories of annotated genes.

Survivorship following electroporation and heat shock

Survivorship of electroporated samples chosen for sequencing ranged from 44 to 86 animals out of 100. When maintained at 20°C, there was no significant difference in survivorship (2 days post-electroporation) between copepods electroporated with target dsRNA versus control/diH2O (Fisher's exact test, n=300 for each control and treatment, P=0.2671; Fig. 1). In contrast, 2 days post-electroporation, survivorship following a 1-h heat shock at 35°C significantly differed between groups electroporated with target dsRNA and control (Fisher's exact test, n=50 for treatment and 53 for control, P=0.0152; Fig. 1).

Fig. 1.

Percent survivorship following electroporation and heat stress (1 h at 35°C 2 days post-electroporation) with control versus HSF dsRNA. A two-tailed Fisher's exact test found no difference in survivorship based on electroporation treatment (n=300 animals in each control and treatment, P=0.3163), while survivorship significantly differed between control and HSF-1 RNAi-treated samples following heat stress (n=30 for treatment and 33 for control, P=2.2×10−16). Survivorship following electroporation is the percentage of animals surviving at 2 days post-electroporation (in the absence of thermal stress), while survivorship following heat stress is the percentage surviving 3 days post-stress. Red diamonds show mean values, and black shapes show replicates 1, 2 and 3 for control RNAi samples (C1, C2, C3) and HSF-1 RNAi samples (T1, T2, T3).

Fig. 1.

Percent survivorship following electroporation and heat stress (1 h at 35°C 2 days post-electroporation) with control versus HSF dsRNA. A two-tailed Fisher's exact test found no difference in survivorship based on electroporation treatment (n=300 animals in each control and treatment, P=0.3163), while survivorship significantly differed between control and HSF-1 RNAi-treated samples following heat stress (n=30 for treatment and 33 for control, P=2.2×10−16). Survivorship following electroporation is the percentage of animals surviving at 2 days post-electroporation (in the absence of thermal stress), while survivorship following heat stress is the percentage surviving 3 days post-stress. Red diamonds show mean values, and black shapes show replicates 1, 2 and 3 for control RNAi samples (C1, C2, C3) and HSF-1 RNAi samples (T1, T2, T3).

HSF and hspb1 expression for knockdown validation

Based on qPCR results, RNAi resulted in only a modest decrease in HSF transcript abundance. Fold change of HSF transcripts in comparison to controls was −1.2-, −1.7- and −1.8-fold lower for the three replicate knockdown samples selected for RNA sequencing (corresponding to samples T1, T2 and T3, respectively). In contrast, hspb1 showed a significant response to HSF knockdown. Log fold change for hspb1 (a presumed HSF target; Tangwancharoen et al., 2018) after heat stress was found to be significantly lower in the treated groups than in the control (t=3.261, P=0.03211; Fig. 2). Notably, the sample with the lowest degree of HSF knockdown (T1) also had the highest upregulation of hspb1, validating the connection between the knockdown degree and hspb1 expression.

Fig. 2.

Log fold change of hspb1 in control and knockdown groups following heat stress (1 h at 35°C) and recovery (1 h at 20°C). There was a significantly lower upregulation of hspb1 in knockdown groups (n=3, P=0.03211). Red diamonds show mean values, and black shapes show replicates 1, 2 and 3 for control RNAi samples (C1, C2, C3) and HSF-1 RNAi samples (T1, T2, T3).

Fig. 2.

Log fold change of hspb1 in control and knockdown groups following heat stress (1 h at 35°C) and recovery (1 h at 20°C). There was a significantly lower upregulation of hspb1 in knockdown groups (n=3, P=0.03211). Red diamonds show mean values, and black shapes show replicates 1, 2 and 3 for control RNAi samples (C1, C2, C3) and HSF-1 RNAi samples (T1, T2, T3).

Protein expression

Protein expression was used in conjunction with qPCR to validate knockdown and select samples to be sent for RNA sequencing. Western blots of the control and HSF-1 RNAi samples show a decrease in HSF protein levels for the samples that were chosen for sequencing in all three cases (Table 2, Fig. S2).

Heat shock element genome scan

Our scan identified 396 unique genes containing at least one canonical HSE within 2000 bp upstream of their start sites. Of these genes, 18 had two or three perfect HSEs. These included genes encoding: SACM1L, a phosphatidylinositide phosphatase; ALDH16A1, an aldehyde dehydrogenase; DNAJA1, a heat shock protein; eIF3a, a component of the eukaryotic translation initiation factor (eIF-3) complex that helps initiate cell proliferation; GLRX3, an iron-sulfur cluster assembly factor; Hr4, a nuclear hormone receptor; INTS10, a component of a complex involved in snRNA transcription; NDUFA9, a subunit of Complex I in the electron transport chain; NKAIN, a sodium ion transport regulator; osa, which can be involved in chromatin regulation and development; sip-1 or hsp beta-1, a small heat shock protein previously shown to be involved thermal tolerance in T. californicus (Barreto et al., 2015; Tangwancharoen et al., 2018); and ttv, a protein involved in development and biosynthesis. For a full list of genes with HSEs, see Table S2.

Illumina sequencing and differential expression analysis

Sequencing resulted in ∼30–50 million reads per sample (Table S1) and an average of 80.5% of the reads mapped to the T. californicus transcriptome. For analysis of differential expression, read counts were normalized in DESeq2 (median of ratios method). Notably, we expect the HSF-1 RNAi knockdown treatment to directly impact the abundance of HSF transcripts; knockdown resulted in a 9% reduction in HSF-1 transcripts in the non-heat-stressed group (5221±309 control and 4740±88 knockdown, n=3 per group) and 7% reduction in the heat-stressed group (5965±489 control and 5491±352 knockdown). Although the data are consistent with the expectation of reduced transcripts for the RNAi target, the magnitude of the knockdown was not significant via a two-way ANOVA, most likely because of the small sample size (knockdown versus control P=0.1996; heat stress versus no heat stress P=0.0598; interaction P=0.9925).

A PCA of the top 1000 differentially expressed genes (Fig. 3) shows general separation along PC1 between the six control and six knockdown samples with the exception of one control (C3) and one knockdown sample (T2). Similarly, there is clear separation between heat stress and control samples along PC2. Owing to the small sample size, we pooled samples across treatments and performed two separate differential expression analyses in order to increase power. First, we analyzed differential expression based on RNAi treatment alone (control versus knockdown). In this comparison, 32 genes were identified as differentially expressed between groups (heat map of differentially expressed genes, see Fig. 4). Of these, only one had a validated canonical HSE within 2000 bp and that gene (TCALIF_02782) was un-annotated. The list of 32 differentially expressed genes in response to knockdown included five probable transcription factors, all of which were downregulated in response to knockdown (TBX20, HLHmbeta, lin-11, fd96Ca and sox14) and one unknown protein, which was upregulated and identified as another probable transcription factor. Additionally, differentially expressed genes included four genes involved in chitin and cuticle formation, including CHIA and three unidentified genes with characteristics of cuticle-associated proteins. All cuticle-associated genes were upregulated in the knockdown. For the full list of genes, see Table S3.

Fig. 3.

Principal component analysis (PCA) of thermal stress and knockdown versus control treatment of the top 1000 differentially expressed genes from DESeq2 analysis. Labels correspond to the three control RNAi (C1, C2, C3) and three HSF-1 RNAi (T1, T2, T3) replicates.

Fig. 3.

Principal component analysis (PCA) of thermal stress and knockdown versus control treatment of the top 1000 differentially expressed genes from DESeq2 analysis. Labels correspond to the three control RNAi (C1, C2, C3) and three HSF-1 RNAi (T1, T2, T3) replicates.

Fig. 4.

Differential expression of control versus knockdown RNA-seq samples. A total of 32 genes were identified as differentially expressed between control and knockdown. C samples, control RNAi treatment; T samples, HSF-1 RNAi treated; HS, heat stress, noHS, no heat stress.

Fig. 4.

Differential expression of control versus knockdown RNA-seq samples. A total of 32 genes were identified as differentially expressed between control and knockdown. C samples, control RNAi treatment; T samples, HSF-1 RNAi treated; HS, heat stress, noHS, no heat stress.

Next, we compared differential expression between heat-stressed and non-heat-stressed groups, with the PCA (Fig. 3) showing separation of thermal stress samples along a horizontal axis. In response to heat stress, 335 genes were differentially expressed (shown in a heat map, see Fig. 5). Of these, 19 had canonical HSEs, and these included heat shock proteins, ion channel and transport proteins, and a transcription factor. See Table S4 for the full list of genes with differential expression in response to heat stress. No genes were differentially expressed in both analyses (i.e. using RNAi treatment as the main factor or using heat as the main factor). However, of the 19 genes with HSEs, eight showed decreased expression (ranging from 2.7 to 51.8% reduction averaged over the three replicates) upon comparison of normalized read counts between heat-stressed control and heat-stressed knockdown groups (see Table S5 for list of differentially expressed genes with HSEs). Though not of sufficient magnitude to be significant in the transcriptome-wide differential expression analysis with DESeq2, these changes in gene expression provide some evidence for an effect that may be attributable to the HSF-1 knockdown. The eight genes include: genes encoding four heat shock proteins, hspb1 and DNAJA1 (an Hsp40), both of which have two heat shock elements, as well as hsp70 and hsp83; STIP1, which coordinates HSP70 function; Slc34a2, a sodium-dependent phosphate transfer protein; lin-11, a transcription factor; and Chrna3, an acetylcholine receptor subunit. Note that in preliminary analyses, we tested for interaction between the two variables (heat stress and HSF-1 RNAi) but found only one significant gene (perhaps owing to the relatively small sample size), which was un-annotated.

Fig. 5.

Differential expression between heat stress and non-heat stress samples. A total of 334 genes showed differential expression. C samples, control RNAi treatment; T samples, RNAi treated; HS, heat stress; noHS, no heat stress.

Fig. 5.

Differential expression between heat stress and non-heat stress samples. A total of 334 genes showed differential expression. C samples, control RNAi treatment; T samples, RNAi treated; HS, heat stress; noHS, no heat stress.

In this study, we used RNAi in T. californicus to knockdown expression of HSF-1, a key transcription factor in the cellular response to thermal stress. HSF knockdown was verified at the transcript level by qPCR and at the protein level by western blot. We then used RNAseq to investigate the role of HSF-1 in the regulation of the entire transcriptome. Though the degree of knockdown was limited, we saw a corresponding decrease in survivorship to heat stress in copepods electroporated with HSF-1 dsRNA, consistent with HSF-1's putative role in heat shock response. Additionally, we saw evidence for differential expression between knockdown and control groups, both independent of heat stress and in a dampened upregulation of genes responding to heat shock. These changes in gene expression provide possible candidates for genes involved in HSF-1 regulation under both stressful and non-stressful conditions.

Interestingly, the majority of genes showing differential expression between knockdown and control were upregulated in the knockdown, indicating that HSF may normally act to repress these genes under both stress and non-stress conditions. Multiple genes involved in chitin and cuticle metabolism showed changes in expression, including a chitinase gene (CHIA), a chitin binding protein gene, a flexible cuticle precursor protein gene and N-acylglucosamine 2-epimerase. The upregulation of a chitinase in response to HSF knockdown indicates that under normal conditions, HSF acts to downregulate this gene, possibly preventing it from breaking down the chitin that makes up the cuticle. During heat stress, cuticle genes have been found to be upregulated in a heat-tolerant population of T. californicus (Schoville et al., 2012); similarly, in Caenorhabditis elegans, cuticle structure genes were identified as the top functional category under HSF-1 control (Brunquell et al., 2016). Owing to their upregulation during the heat shock response and the repression of chitinase under normal conditions, cuticle gene activity is evidently important in Tigriopus, though the function of these changes in gene expression remains to be elucidated. Genes associated with aging have been linked to cuticle structure genes in C. elegans (Brunquell et al., 2016), suggesting that processes related to longevity may be involved in the regulation of cuticle genes observed in T. californicus.

Of the 32 genes that were differentially expressed between the knockdown and control groups, eight were downregulated in the knockdown. These downregulated genes include several probable transcription factors: TBX20, HLHmbeta, lin-11, fd96Ca and sox14. Their downregulation in the knockdown group suggests that they are normally upregulated by HSF-1. HSF-1 has been shown to interact with other transcription factors in a range of organisms (e.g. Brunquell et al., 2016; Stephanou et al., 1999), including other forms of HSF in some organisms (Takii et al., 2017). These types of interactions with other transcription factors complicate the role of HSF-1 and may explain why we do not see as clear a relationship between knockdown and heat shock response as we might have expected.

In addition to RNAseq analysis, we also performed a genome scan to identify genes with potential HSF-1 binding sites (HSEs) within 2000 bp of their start sites. While this list of genes provides a starting point for future analysis, we did not find genes with HSEs to be significantly over-represented in our differential expression analysis comparing knockdown and control samples. This appears to support the idea, as other studies have found, that HSF-1 does not necessarily need canonical HSEs for binding and may even bind more successfully to degenerate HSEs (Guertin et al., 2012). Li et al. (2016) found evidence that HSF-1 is more likely to bind to canonical HSEs following heat shock, while preferring degenerate HSEs during development. This distinction could perhaps explain why many genes that were differentially expressed in our knockdown did not have canonical HSEs, as the genes found in this analysis were separate from the genes that were differentially expressed between heat-stressed and non-heat-stressed groups. Similarly, canonical HSEs are not always sufficient for binding HSF-1 even after heat shock; using chromatin immunoprecipitation followed by sequencing (ChIP-seq) analysis, Guertin and Lis (2010) found that many HSE motifs were not occupied by HSF. Furthermore, even when HSF is strongly bound to an HSE, it does not lead to gene activation in some cases, and may require additional downstream promoter features for induction (Guertin and Lis, 2010; Vihervaara et al., 2017).

Upon comparison of heat shock and non-heat shock, regardless of knockdown or control treatments, 334 genes were differentially expressed, which is similar to the DEG count of 356 observed in a previous heat shock study in T. californicus (Schoville et al., 2012). Of these, 107 were upregulated following heat stress. Additionally, 19 genes of the 334 had associated canonical HSEs within 2 kb of their transcription start sites, with nine showing downregulation in response to heat stress and 10 showing upregulation. Though none of the genes showing differential expression following heat shock (comparing all heat shock samples with all non-heat shock samples) were genes that also showed differential expression following knockdown (comparing all knockdown samples with all control samples), we did see evidence for an effect of HSF-1 knockdown on gene expression when we compared the level of upregulation between knockdown and control groups in genes with HSEs: of the 19 genes with canonical HSEs in our analysis of differential expression between heat-stressed and non-heat-stressed groups, five showed decreased upregulation with knockdown treatment. All five are involved in the heat shock response: STIP1, which acts as a co-chaperone for HSP90; hsp70 and hsp83; DNAJA1, an HSP40; and sip-1, or hspb1. While the decreased expression of these genes was not strong enough to be detected by the transcriptome-wide DESeq2 analysis of knockdown versus control groups (likely because of their high upregulation following heat stress regardless of knockdown), their decreased read counts in the knockdown samples suggests that HSF-1 plays a role in their upregulation following heat shock and that the presence of canonical HSEs may be important in driving their upregulation.

Despite the compelling examples (mentioned above) that exceptions to the rule of HSF-1 interacting with canonical HSEs exist, previous studies in T. californicus have supported our result that canonical HSEs may be important in mounting an effective heat stress response. Tangwancharoen et al. (2018) showed that in F1 hybrids between a thermally tolerant (San Diego) and thermally sensitive (Santa Cruz) population, expression of the San Diego HSPB1 allele was 2- to 3-fold higher than expression of the Santa Cruz allele in the same animals. Because the same trans elements (HSF-1) exist in each case, the difference in expression is likely due to the substitutions in the HSEs preceding hspb1 between the two populations (i.e. cis regulation). Because HSPB1 has been shown to be important in heat stress survivorship (Barreto et al., 2015; Tangwancharoen et al., 2018), this allele-specific expression is a good indication that canonical HSEs are important in the heat shock response of T. californicus. Our findings may also support this conclusion for additional heat shock protein genes. Furthermore, both DNAJA1 and hspb1 have two canonical HSEs in their promoters, which may increase their level of control from HSF-1. In support of this result, Li et al. (2016) similarly found that genes involved in heat shock were more likely to have tandem canonical HSEs that bind HSF-1. Taken together, our findings that many differentially expressed genes following knockdown (that may not be involved in heat stress) lack a canonical HSE and that there is decreased upregulation after heat shock following knockdown are both supported by previous studies.

In sum, despite the limited magnitude of HSF-1 knockdown achieved in these experiments, we identified some genes with differential expression following a decrease of HSF-1. More complete knockdown of HSF would likely expand this list. Brunquell et al. (2016) used RNAi feeding with larval C. elegans, which appears to be an effective method of RNAi delivery in that system; although electroporation of T. californicus adults has been optimized and validated for hspb1 (Barreto et al., 2015), this gene shows very strong upregulation (∼100-fold) following thermal stress (Schoville et al., 2012). In contrast, upregulation of HSF-1 following heat treatment is far less pronounced (∼1.3-fold), and knockdown by RNAi can be difficult to validate as continued low-level expression may obscure physiological consequences of the knockdown.

One clear takeaway from the results of our RNA-seq is the complexity of HSF-1 interactions. Despite generating a list of genes with perfect HSEs, we did not see evidence for differential expression based on HSF-1 knockdown alone in any of these genes. Although we did see decreased upregulation following heat shock in five HSE-adjacent genes, the reduction was relatively small. Using ChIP-seq to identify HSF binding locations in Drosophila, Guertin and Lis (2010) found that only a small proportion of potential HSE motifs were occupied by HSF, and that even binding of HSF to HSEs did not always lead to gene activation. From our results, it seems likely that HSF may also bind to non-canonical HSEs and activate adjacent genes. It is also possible that some of the other genes it acts on, such as the transcription factors we identified, have their own control cascades, further complicating the interactions. In the future, a ChIP-seq approach in T. californicus could help clarify these interactions.

There are a variety of factors that may have impacted the low degree of knockdown and subsequent effects on gene expression we observed. First, our electroporation method utilizes adult copepods, and therefore likely does not deliver dsRNA to all cells within the organism. Therefore, we are only seeing effects of knockdown in the cells in which it has entered. Second, our timing of sample collection may have resulted in a dampening of response. While we found that heat stress sooner than 2 days following electroporation was too stressful and resulted in high animal mortality, even at low temperatures, this may indicate that the effects of knockdown are not as strong at the time point. Similarly, we euthanized animals for RNA extraction after a 1-h heat stress and 1-h recovery. Though this time point may have been a good balance between HSF response at the mRNA level (which has been suggested to happen rapidly in T. californicus in recent studies in the laboratory) and gene response, recent evidence from the laboratory suggests that a longer heat stress and recovery time may be needed in order to observe higher upregulation of genes in response to heat stress (e.g. Harada and Burton, 2019). Thus, we may not be seeing the maximal effect of knockdown on these genes.

Studying the complex interactions between HSF and the genes it controls can bring us closer to understanding the role of gene regulation in thermal tolerance, in addition to other processes related to stress, development and longevity. From an evolutionary perspective, it is interesting to note that there is substantial inter-population variation in the HSF protein sequence: the HSF-1 allele from the northern Santa Cruz T. californicus population has 10 amino acid substitutions compared with the southern San Diego population (Tangwancharoen et al., 2018). It is possible that these substitutions change the efficacy of HSF's regulation of transcription during heat stress or in non-stress scenarios in the northern populations. Further investigation into these pathways may shed light on the functional consequences of these substitutions in the future.

The authors thank the members of the Burton Lab for their help and support, especially Dr Sumaetee Tangwancharoen, Dr Timothy Healy and Gary Moy for assistance and advice.

Author contributions

Conceptualization: A.E.H., R.S.B.; Methodology: A.E.H., R.S.B.; Software: A.E.H., R.S.B.; Validation: A.E.H., R.S.B.; Formal analysis: A.E.H., R.S.B.; Investigation: A.E.H., R.S.B.; Resources: A.E.H., R.S.B.; Data curation: A.E.H., R.S.B.; Writing - original draft: A.E.H., R.S.B.; Writing - review & editing: A.E.H., R.S.B.; Visualization: A.E.H., R.S.B.; Supervision: R.S.B.; Project administration: R.S.B.; Funding acquisition: R.S.B.

Funding

Funding for this study was provided by a National Science Foundation grant to R.S.B. (DEB1551466).

Data availability

Data from this study are available in BioProject under accession number PRJNA603100.

Åkerfelt
,
M.
,
Morimoto
,
R. I.
and
Sistonen
,
L.
(
2010
).
Heat shock factors: integrators of cell stress, development and lifespan
.
Nat. Rev. Mol. Cell Biol.
11
,
545
-
555
.
Anckar
,
J.
and
Sistonen
,
L.
(
2007
).
Heat shock factor 1 as a coordinator of stress and developmental pathways
. In
Molecular Aspects of the Stress Response: Chaperones, Membranes, and Networks
(ed.
P.
Csermely
and
L.
Vígh
), pp.
96
-
106
.
New York
:
Springer
.
Ashburner
,
M.
and
Bonner
,
J. J.
(
1979
).
The induction of gene activity in Drosophila by heat shock
.
Cell
17
,
241
-
254
.
Barreto
,
F. S.
,
Schoville
,
S. D.
and
Burton
,
R. S.
(
2015
).
Reverse genetics in the tide pool: knock-down of target gene expression via RNA interference in the copepod Tigriopus californicus
.
Mol. Ecol. Resources
15
,
868
-
879
.
Barreto
,
F. S.
,
Watson
,
E. T.
,
Lima
,
T. G.
,
Willett
,
C. S.
,
Edmands
,
S.
,
Li
,
W.
and
Burton
,
R. S.
(
2018
).
Genomic signatures of mitonuclear coevolution across populations of Tigriopus californicus
.
Nat. Ecol. Evol.
2
,
1250
-
1257
.
Barshis
,
D. J.
,
Ladner
,
J. T.
,
Oliver
,
T. A.
,
Seneca
,
F. O.
,
Traylor-Knowles
,
N.
and
Palumbi
,
S. R.
(
2013
).
Genomic basis for coral resilience to climate change
.
Proc. Natl Acad. Sci. USA
110
,
1387
-
1392
.
Brunquell
,
J.
,
Morris
,
S.
,
Lu
,
Y.
,
Cheng
,
F.
and
Westerheide
,
S. D.
(
2016
).
The genome-wide role of HSF-1 in the regulation of gene expression in Caenorhabditis elegans
.
BMC Genomics
17
,
1
-
18
.
Feder
,
M. E.
and
Hofmann
,
G. E.
(
1999
).
Heat-shock proteins, molecular chaperones, and the stress response: evolutionary and ecological physiology
.
Annu. Rev. Physiol.
61
,
243
-
282
.
Fernandes
,
M.
,
Xiao
,
H.
and
Lis
,
J. T.
(
1994
).
Fine structure analyses of the Drosophila and Saccharomyces heat shock factor - heat shock element interactions
.
Nucleic Acids Res.
22
,
167
-
173
.
Furuhashi
,
T.
and
Sakamoto
,
K.
(
2015
).
Heat shock factor 1 prevents the reduction in thrashing due to heat shock in Caenorhabditis elegans
.
Biochem. Biophys. Res. Commun.
462
,
190
-
194
.
Gleason
,
L. U.
and
Burton
,
R. S.
(
2015
).
RNA-seq reveals regional differences in transcriptome response to heat stress in the marine snail Chlorostoma funebralis
.
Mol. Ecol.
24
,
610
-
627
.
Guertin
,
M. J.
and
Lis
,
J. T.
(
2010
).
Chromatin landscape dictates HSF binding to target DNA elements
.
PLoS Genet.
6
,
e1001114
.
Guertin
,
M. J.
,
Martins
,
A. L.
,
Siepel
,
A.
and
Lis
,
J. T.
(
2012
).
Accurate prediction of inducible transcription factor binding intensities in vivo
.
PLoS Genet.
8
,
e1002610
.
Harada
,
A. E.
and
Burton
,
R. S.
(
2019
).
Ecologically relevant temperature ramping rates enhance the protective heat shock response in an intertidal ectotherm
.
Physiol. Biochem. Zool.
92
,
152
-
162
.
Hoffmann
,
A. A.
,
Sørensen
,
J. G.
and
Loeschcke
,
V.
(
2003
).
Adaptation of Drosophila to temperature extremes: Bringing together quantitative and molecular approaches
.
J. Therm. Biol.
28
,
175
-
216
.
Jedlicka
,
P.
,
Mortin
,
M. A.
and
Wu
,
C.
(
1997
).
Multiple functions of Drosophila heat shock transcription factor in vivo
.
EMBO J.
16
,
2452
-
2462
.
Kregel
,
K. C.
(
2002
).
Heat shock proteins: modifying factors in physiological stress responses and acquired thermotolerance
.
J. Appl. Physiol.
92
,
2177
-
2186
.
Labbadia
,
J.
,
Brielmann
,
R. M.
,
Neto
,
M. F.
,
Lin
,
Y.-F.
,
Haynes
,
C. M.
and
Morimoto
,
R. I.
(
2017
).
Mitochondrial stress restores the heat shock response and prevents proteostasis collapse during aging
.
Cell Rep.
21
,
1481
-
1494
.
Li
,
J.
,
Chauve
,
L.
,
Phelps
,
G.
,
Brielmann
,
R. M.
and
Morimoto
,
R. I.
(
2016
).
E2F coregulates an essential HSF developmental program that is distinct from the heat-shock response
.
Genes Dev.
30
,
2062
-
2075
.
Lindquist
,
S.
(
1986
).
The heat-shock response
.
Annu. Rev. Biochem.
55
,
1151
-
1191
.
Lockwood
,
B. L.
,
Sanders
,
J. G.
and
Somero
,
G. N.
(
2010
).
Transcriptomic responses to heat stress in invasive and native blue mussels (genus Mytilus): molecular correlates of invasive success
.
J. Exp. Biol.
213
,
3548
-
3558
.
Loomis
,
W. F.
and
Wheeler
,
S.
(
1980
).
Heat shock response of Dictyostelium
.
Dev. Biol.
79
,
399
-
408
.
Love
,
M. I.
,
Huber
,
W.
and
Anders
,
S.
(
2014
).
Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2
.
Genome Biol.
15
,
550
.
Mazin
,
P. V.
,
Shagimardanova
,
E.
,
Kozlova
,
O.
,
Cherkasov
,
A.
,
Sutormin
,
R.
,
Stepanova
,
V. V.
,
Stupnikov
,
A.
,
Logacheva
,
M.
,
Penin
,
A.
,
Sogame
,
Y.
, et al.
(
2018
).
Cooption of heat shock regulatory system for anhydrobiosis in the sleeping chironomid Polypedilum vanderplanki
.
Proc. Natl Acad. Sci. USA
115
,
E2477
-
E2486
.
Morimoto
,
R. I.
(
1998
).
Regulation of the heat shock transcriptional response: cross talk between a family of heat shock factors, molecular chaperones, and negative regulators
.
Genes Dev.
12
,
3788
-
3796
.
Morton
,
E. A.
and
Lamitina
,
T.
(
2013
).
Caenorhabditis elegans HSF-1 is an essential nuclear protein that forms stress granule-like structures following heat shock
.
Aging Cell
12
,
112
-
120
.
Pereira
,
R. J.
,
Sasaki
,
M. C.
and
Burton
,
R. S.
(
2017
).
Adaptation to a latitudinal thermal gradient within a widespread copepod species: the contributions of genetic divergence and phenotypic plasticity
.
Proc. R. Soc. B
284
,
20170236
.
Pirkkala
,
L.
,
Nykänen
,
P.
and
Sistonen
,
L.
(
2001
).
Roles of the heat shock transcription factors in regulation of the heat shock response and beyond
.
FASEB J.
15
,
1118
-
1131
.
Sarge
,
K. D.
,
Murphy
,
S. P.
and
Morimoto
,
R. I.
(
1993
).
Activation of heat shock gene transcription by heat shock factor 1 involves oligomerization, acquisition of DNA-binding activity, and nuclear localization and can occur in the absence of stress
.
Mol. Cell. Biol.
13
,
1392
-
1407
.
Schmittgen
,
T. D.
and
Livak
,
K. J.
(
2008
).
Analyzing real-time PCR data by the comparative CT method
.
Nat. Protoc.
3
,
1101
-
1108
.
Schoville
,
S. D.
,
Barreto
,
F. S.
,
Moy
,
G. W.
,
Wolff
,
A.
and
Burton
,
R. S.
(
2012
).
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.
12
,
170
.
Seo
,
K.
,
Choi
,
E.
,
Lee
,
D.
,
Jeong
,
D.-E.
,
Jang
,
S. K.
and
Lee
,
S.-J.
(
2013
).
Heat shock factor 1 mediates the longevity conferred by inhibition of TOR and insulin/IGF-1 signaling pathways in C. elegans
.
Aging Cell
12
,
1073
-
1081
.
Sørensen
,
J. G.
,
Kristensen
,
T. N.
and
Loeschcke
,
V.
(
2003
).
The evolutionary and ecological role of heat shock proteins
.
Ecol. Lett.
6
,
1025
-
1037
.
Sørensen
,
J. G.
,
Nielsen
,
M. M.
,
Kruhøffer
,
M.
,
Justesen
,
J.
and
Loeschcke
,
V.
(
2005
).
Full genome gene expression analysis of the heat stress response in Drosophila melanogaster
.
Cell Stress Chaperones
10
,
312
-
328
.
Sorger
,
P. K.
(
1991
).
Heat shock factor and the heat shock response
.
Cell
65
,
363
-
366
.
Steinkraus
,
K. A.
,
Smith
,
E. D.
,
Davis
,
C.
,
Carr
,
D.
,
Pendergrass
,
W. R.
,
Sutphin
,
L.
,
Kennedy
,
B. K.
and
Kaeberlein
,
M.
(
2008
).
Dietary restriction suppresses proteotoxicity and enhances longevity by an hsf-1-dependent mechanism in Caenorhabditis elegans
.
Aging Cell
7
,
394
-
404
.
Stephanou
,
A.
,
Isenberg
,
D. A.
,
Nakajima
,
K.
and
Latchman
,
D. S.
(
1999
).
Signal transducer and activator of transcription-1 and heat shock factor-1 interact and activate the transcription of the Hsp-70 and Hsp-90β gene promoters
.
J. Biol. Chem.
274
,
1723
-
1728
.
Stillman
,
J. H.
and
Tagmount
,
A.
(
2009
).
Seasonal and latitudinal acclimatization of cardiac transcriptome responses to thermal stress in porcelain crabs, Petrolisthes cinctipes
.
Mol. Ecol.
18
,
4206
-
4226
.
Takii
,
R.
,
Fujimoto
,
M.
,
Matsuura
,
Y.
,
Wu
,
F.
,
Oshibe
,
N.
,
Takaki
,
E.
,
Katiyar
,
A.
,
Akashi
,
H.
,
Makino
,
T.
,
Kawata
,
M.
, et al.
(
2017
).
HSF1 and HSF3 cooperatively regulate the heat shock response in lizards
.
PLoS ONE
12
,
e0180776
.
Tangwancharoen
,
S.
,
Moy
,
G. W.
and
Burton
,
R. S.
(
2018
).
Multiple modes of adaptation: regulatory and structural evolution in a small heat shock protein gene
.
Mol. Biol. Evol.
35
,
2110
-
2119
.
Vihervaara
,
A.
,
Mahat
,
D. B.
,
Guertin
,
M. J.
,
Chu
,
T.
,
Danko
,
C. G.
,
Lis
,
J. T.
and
Sistonen
,
L.
(
2017
).
Transcriptional response to stress is pre-wired by promoter and enhancer architecture
.
Nat. Commun.
8
,
255
.
Walker
,
G. A.
,
Thompson
,
F. J.
,
Brawley
,
A.
and
Scanlon
,
T.
(
2003
).
Heat shock factor functions at the convergence of the stress response and developmental pathways in Caenorhabditis elegans
.
FASEB J.
17
,
1960
-
1962
.
Westerheide
,
S. D.
,
Raynes
,
R.
,
Powell
,
C.
,
Xue
,
B.
and
Uversky
,
V. N.
(
2012
).
HSF transcription factor family, heat shock response, and protein intrinsic disorder
.
Curr. Protein Pept. Sci.
13
,
86
-
103
.
Westwoodt
,
J. T.
and
Wu
,
C.
(
1993
).
Activation of Drosophila heat shock factor: conformational change associated with a monomer-to-trimer transition
.
Mol. Cell. Biol.
13
,
3481
-
3486
.
Willett
,
C. S.
and
Burton
,
R. S.
(
2001
).
Viability of cytochrome c genotypes depends on cytoplasmic backgrounds in Tigriopus californicus
.
Evolution
55
,
1592
-
1599
.
Wu
,
C.
(
1995
).
Heat shock transcription factors: structure and regulation
.
Annu. Rev. Cell Dev. Biol.
11
,
441
-
469
.
Xiao
,
X.
,
Zuo
,
X.
,
Davis
,
A. A.
,
McMillan
,
D. R.
,
Curry
,
B. B.
,
Richardson
,
J. A.
and
Benjamin
,
I. J.
(
1999
).
HSF1 is required for extra-embryonic development, postnatal growth and protection during inflammatory responses in mice
.
EMBO J.
18
,
5943
-
5952
.
Zhong
,
M.
,
Orosz
,
A.
and
Wu
,
C.
(
1998
).
Direct sensing of heat and oxidation by Drosophila heat shock transcription factor
.
Mol. Cell
2
,
101
-
108
.
Zhu
,
L.
(
2013
).
Integrative analysis of ChIP-chip and ChIP-seq dataset
. In
Methods in Molecular Biology (Methods and Protocols): Tiling Arrays
, Vol.
1067
(ed.
T.
Lee
and
A.
Shui Luk
), pp.
105
-
124
.
Totowa, NJ
:
Humana Press
.
Zhu
,
L. J.
,
Gazin
,
C.
,
Lawson
,
N. D.
,
Pagès
,
H.
,
Lin
,
S. M.
,
Lapointe
,
D. S.
and
Green
,
M. R.
(
2010
).
ChIPpeakAnno: a Bioconductor package to annotate ChIP-seq and ChIP-chip data
.
BMC Bioinformatics
11
,
237
.

Competing interests

The authors declare no competing or financial interests.

Supplementary information