Blubber transcriptome responses to repeated ACTH administration in a marine mammal

Chronic physiological stress impacts animal fitness by catabolizing metabolic stores and suppressing reproduction. This can be especially deleterious for capital breeding carnivores such as marine mammals, with potential for ecosystem-wide effects. However, the impacts and indicators of chronic stress in animals are currently poorly understood. To identify downstream mediators of repeated stress responses in marine mammals, we administered adrenocorticotropic hormone (ACTH) once daily for four days to free-ranging juvenile northern elephant seals (Mirounga angustirostris) to stimulate endogenous corticosteroid release, and compared blubber tissue transcriptome responses to the first and fourth ACTH administrations. Gene expression profiles were distinct between blubber responses to single and repeated ACTH administration, despite similarities in circulating cortisol profiles. We identified 61 and 12 genes that were differentially expressed (DEGs) in response to the first ACTH and fourth administrations, respectively, 24 DEGs between the first and fourth pre-ACTH samples, and 12 DEGs between ACTH response samples from the first and fourth days. Annotated DEGs were associated with functions in redox and lipid homeostasis, suggesting potential negative impacts of repeated stress on capital breeding, diving mammals. DEGs identified in this study are potential markers of repeated stress in marine mammals, which may not be detectable by endocrine profiles alone.

suppression of energetically expensive processes, such as reproduction, growth, and immune responses 9 . While short-term (acute) stress responses are adaptive [14][15][16] , repeated (chronic) stress exposure may deplete energy stores and impair immunity and reproduction, reducing fitness. This could be detrimental to threatened animal populations that already experience nutritional stress and low fecundity, such as some marine mammals 17 . However, little is known about the indicators and downstream effects of repeated stress in wildlife, hindering the ability of conservation biologists to identify stressed animals and predict the effects of repeated stress on population stability 18,19 .
Marine mammals rely on endogenous energy reserves to support life history stages characterized by high energy expenditure and nutrient limitation. For example, semiaquatic marine mammals (e.g. pinnipeds -seals and sea lions) reproduce on land using a capital breeding strategy, in which animals rely solely on energy stores to fuel metabolic demands of reproduction 20,21 . Many fully aquatic marine mammals (e.g. cetaceans) are also capital breeders that migrate long distances between foraging and breeding grounds 22 . Both examples necessitate energy accumulation and storage during feeding and energy mobilization during fasting, a process that occurs for several weeks to several months as part of marine mammal life history. The main energy depot in cetaceans and pinnipeds is blubber, a specialized type of subcutaneous adipose tissue that plays significant roles in metabolism (similarly to white adipose tissue, WAT, in other mammals) and thermoregulation. Blubber is vertically stratified into at least two layers: an outer layer used for thermoregulation, and an inner, metabolically active layer 23 . WAT is a major target of stress hormones in mammals 24 and studies in marine mammals have shown that acute HPA axis activation mobilizes stored lipids 10,25 . It is plausible that chronic or repeated stress would deregulate fasting metabolism and deplete lipid stores, among other impacts, limiting the ability of capital breeders to meet the energetic demands of reproduction and leading to declines in fitness.
In this study, we evaluated the effects of repeated HPA axis activation on blubber gene expression in a fasting-adapted marine mammal, the northern elephant seal (Mirounga angustirostris). Elephant seals are an established marine mammal study system in which variation in baseline hormone levels 26 , responses to acute ACTH administration 14,25,27,28 , and metabolic adjustments during fasting 29,30 have been described. In addition, anesthesia procedures have been shown to minimize handling stress in this species, enabling measurements of true baseline hormone and gene expression levels 25 . We simulated repeated physiological stress exposure by administering synthetic ACTH to juvenile elephant seals once daily for 4 days, inducing significant elevations in corticosteroid levels that we described in a previous publication 12 . To evaluate downstream mechanisms by which repeated HPA axis activation may impact metabolic homeostasis, we compared transcriptome responses of elephant seal blubber tissue to the first and last ACTH administration by RNA sequencing (RNAseq). We identified 27 genes that were differentially expressed in blubber during the response to the fourth versus the first ACTH administration, many of which were associated with lipid storage and mobilization. These genes are potential biomarkers for discriminating between acute and repeated ACTH responses in marine mammals, which otherwise cannot be distinguished by cortisol profiles alone. All procedures involving animals were conducted in accordance with the relevant guidelines and regulations of IACUC and BUMED protocols and the NOAA permit. Repeated stress was simulated by administering ACTH to each study animal once every 24 hours for 4 consecutive days. Blood and blubber samples were collected before ACTH administration ("pre-ACTH") and 4 hours after ACTH administration ("ACTH response") on the first and fourth days (Fig. 1). Details of the experimental manipulation and hormone response have been reported previously 12 . In summary, ACTH administration induced significant elevation in circulating cortisol and aldosterone levels and decline in total triiodothyronine (T3) levels. Cortisol levels recovered to baseline within 24 hours of each ACTH administration and cortisol responses to ACTH did not vary in magnitude between the first and fourth administrations. However, aldosterone responses to ACTH were higher on day 4 than on day 1 and did not recover to baseline between administrations (for details, see McCormley et al. 12 ). Four of the 7 animals that participated in the repeated ACTH administration experiment described in McCormley et al. 12 (seals 2, 4, 6, and 7) were selected for the current study based on amounts of sequencing quality RNA that could be obtained from blubber samples.

Materials and
ACtH administration and blubber sampling. Study animals were chemically immobilized and blood samples were collected as previously described 12 . Blubber samples ("day 1 pre-ACTH") were collected from the posterior flank region of each animal using a 6.0 mm diameter biopsy punch (Miltex, USA). The inner half (closest to musculature) of each biopsy sample was isolated, minced, and placed in cryovials containing RNAlater ™ Stabilization Solution (1.5 mL per ~300 mg tissue; Invitrogen, USA). Samples were incubated for 24 hours at 4 °C, after which RNAlater ™ solution was removed and samples were stored at −80 °C until RNA isolation. After baseline sampling, animals received 20 U of corticotropin LA gel (Wedgewood Pharmacy, USA) via intramuscular injection (mass specific dose: 0.17 ± 0.02 U/kg). A paired set of blubber and blood samples was collected 4 hours after ACTH administration ("day 1 ACTH response"); biopsies were collected from the contralateral side of the animal to the ACTH injection. On the second and third days of the experiment, ACTH was administered approximately 24 hours after the injection from the previous day, but no post-ACTH sampling was conducted. On the www.nature.com/scientificreports www.nature.com/scientificreports/ fourth day, pre-ACTH sampling ("day 4 pre-ACTH"), ACTH administration, and post-ACTH sampling ("day 4 ACTH response") were conducted as described for the first day.
RNA isolation. Blubber samples were homogenized by bead beating in Qiazol (Qiagen, USA; ~100 mg tissue per 0.5 mL) using Bullet Blender Storm 24 (Next Advance, USA; Speed 12, two 2-minute cycles). After bead beating, another 0.5 mL of Qiazol was added to each sample and incubated for 5 minutes at room temperature with occasional vortexing. Lysates were further disrupted using QiaShredder tubes (Qiagen, USA) to shear genomic DNA and centrifuged for 15 minutes to separate lipids and cellular components. RNA was isolated from homogenates using the Lipid RNeasy Tissue Kit (Qiagen, USA) after chloroform (VWR Life Sciences, USA) phase extraction according to manufacturer's protocol. Genomic DNA was removed using a 20-minute on-column DNase I digest (Qiagen, USA). RNA concentration was quantified using Qubit 3.0 Fluorometer Broad Range RNA Assay (Life Technologies, USA). RNA integrity was evaluated using the Total RNA 6000 Pico kit on the 2100 Bioanalyzer (Agilent Technologies, USA). Mean (±SD) RNA integrity number (RIN) and 28S/18S rRNA abundance ratios were 6.81 ± 0.32 and 1.03 ± 0.27, respectively. We used Qiagen RNAlater ™ for preservation of blubber samples and isolation of RNA with RIN ≥ 7 in a previous study 14 . However, Ambion RNAlater ™ reagent used in this study did not adequately prevent RNA degradation. Due to the logistical challenges of the experiment, blubber sampling could not be repeated and we proceeded with RNA sequencing, using ribosomal RNA depletion instead of polyA mRNA enrichment during library preparation to reduce 3′ degradation bias 31 . Adequate amounts sequencing-quality RNA were obtained from 14 samples: the entire sample set (day 1 pre-ACTH, day 1 ACTH response, day 4 pre-ACTH, day 4 ACTH response) from seals 4, 6 and 7, and two samples (day 1 ACTH response, day 4 ACTH response) from seal 2.
Library preparation and RNA sequencing. cDNA library preparation and Illumina HiSeq 4000 sequencing were conducted at the University of California, Berkeley QB3 Vincent J. Coates Genomics Sequencing Laboratory. Ribosomal RNA (rRNA) was depleted using Ribo-Zero rRNA Removal Kit (Human/Mouse/Rat; Illumina, Inc., USA). Yield and quality of mRNA were assessed using fluorometric methods and an Agilent 2100 Bioanalyzer, respectively (Agilent Technologies, USA). PrepX directional RNAseq library kits (WaferGen Bio-systems Inc., USA) were used for library preparation and libraries were quantified with KAPA Library Quantification kits for Illumina platforms (Kapa Biosystems, USA) using a Bio-Rad CFX Connect (Bio-Rad Laboratories, USA). Thirteen cycles of indexing PCR using KAPA High Fidelity Hotstart Amplification Kits (Kapa Biosystems, USA) were used before library quantification and validation. Individual libraries were barcoded, indexed as TruSeq single index per library, pooled in equimolar concentrations (400 ng), and sequenced as 100 base-pair paired-end reads on the Illumina HiSeq 4000. This generated a total of 150,055 sequenced megabases and a mean (±SD) of 26.9 ± 4.2 million paired-end 100-bp reads per sample. Sequencing data were demultiplexed and converted from bclfile to fastq file with Illumina's bcl2Fastq software v2.18 (Illumina, USA). Raw data were uploaded to NCBI Sequence Read Archive (SRA accession: SRP157071).
De novo transcriptome assembly. All computational analyses were conducted using the Extreme Science and Engineering Discovery Environment (XSEDE) 32 Bridges Large High Performance Computing Cluster at the Pittsburg Supercomputing Center through allocation TG-IBN150010. Transcriptome assembly was conducted using Trinity v2.4.0 33 . Sequencing adapters and poor quality bases were trimmed from reads using Trimmomatic  12 ) and blubber tissue sampling design. ACTH was administered to juvenile elephant seals once every 24 hours for 4 days. Blubber biopsies were collected immediately before ACTH administration ("pre-ACTH") and 4 hours after ACTH administration ("ACTH response") on the first and fourth days. RNA isolated from blubber tissue was used for transcriptome sequencing. Seal (color modified, from https://pixabay.com), syringe (from https://pixabay. com), adipose (enhanced contrast, from https://openclipart.org), and microcentrifuge tube (from http://www. clker.com) images were obtained under CC0 1.0 Universal: CC0 1.0 Public Domain Dedication. License: https:// creativecommons.org/publicdomain/zero/1.0.
www.nature.com/scientificreports www.nature.com/scientificreports/ run in Trinity with default settings. Quality of raw and trimmed reads was evaluated using FastQC v0.11.3 34 . De novo assembly was conducted after in silico read normalization (50X coverage, kmer size 25 bp) using Trinity with all default settings (paired-end, strand-specific FR mode). Read mapping metrics were obtained using Bowtie2 v2.2.7 35 and transcriptome completeness was analyzed using BUSCO v1.22 36 (metazoan BUSCO dataset downloaded on 8/24/2017). transcriptome annotation. The assembly was annotated by aligning translated DNA query sequences against the UniProtKB/SwissProt protein sequence database (downloaded on 9/29/2017; "BLASTX") using DIAMOND v.0.8.31 37 with e-value threshold for significant matches of 1e-3 and the "more-sensitive" option. Putative protein sequences encoded in the transcriptome were predicted using TransDecoder v3.0.1 38 and annotated against the UniProtKB/SwissProt protein sequence database ("BLASTP") using DIAMOND with e-value cutoff of 1e-3. Differentially expressed genes that did not have hits to the SwissProt database were manually annotated against the NCBI RefSeq database using NCBI BLASTX v2.7.1 with an e-value cutoff of 1e-5. Functional annotation of KEGG categories and GO enrichment was performed using DAVID Bioinformatics Resources v6.8 39 with the human genome as background. Categories were considered significantly enriched at p < 0.05 (adjusted for multiple comparisons using Benjamini-Hochberg correction 40 ).

Gene expression analyses.
All gene expression analyses were conducted using the Trinity pipeline.
Transcript abundance was estimated using Kallisto v0.43.0 41 . Differential expression analysis was performed at the gene level using the DESeq2 package in Bioconductor v3.5 42 in R v3.4.1 with false discovery rate cutoff of 0.05 and log2 fold-change cutoff of 1.0. Pairwise comparisons between sampling conditions are shown in Table 3. Libraries from seal 2 were included only in the ACTH response comparison (day 4 ACTH response vs day 1 ACTH response). All other comparisons were conducted using libraries from seals 4, 6, and 7 only. Protein-protein interaction network prediction was conducted using STRING v10.5 43 . Network data were imported into Cytoscape v3.5.0 44 and filtered by experimental interactions using the edge weighted spring embedded layout. Unconnected nodes were removed from the analysis. Cytoskape network statistics were calculated using the network analyzer tool with an undirected analysis.

Results
Blubber transcriptome assembly. Blubber biopsies were collected from 4 juvenile northern elephant seals immediately prior to ACTH administration ("pre-ACTH") and 4 hours following ACTH administration ("ACTH response") on the first and fourth days of the experiment (Fig. 1). Morphometric data, ACTH doses, and corticosteroid concentrations measured at the four sampling points are presented in Table 1 (from McCormley et al. 12 ). Fourteen cDNA libraries from day 1 pre-ACTH (seals 4, 6, and 7), day 1 ACTH response (seals 2, 4, 6, and 7), day 4 pre-ACTH (seals 4, 6, and 7), and day 4 ACTH response (seals 2, 4, 6, and 7) were sequenced using two Illumina HiSeq 4000 lanes. The libraries from seal 2 were used for assembly but were included only in the ACTH response comparison to increase power for detection of differentially expressed genes (see Methods). Raw sequencing reads were uploaded to NCBI Short Read Archive (BioProject ID: PRJNA485363, SRA accession: SRP157071). A single reference transcriptome was assembled de novo using Trinity software. We assembled 1.36 billion bases into 2,031,456 contigs (or "transcripts") in 1,216,779 gene clusters ( Table 2). The high number of transcripts is common for de novo assemblies 45 and is likely a result of sequence polymorphism, alternative splicing, and variability in individuals' responses to ACTH. The transcriptome assembly is available at figshare: (https://figshare.com/s/37cba07b2e9877f25b50).
To evaluate read representation in the reference transcriptome, we mapped sequenced reads back to the assembly. Ninety-two percent of reads mapped back as proper pairs to the assembly, suggesting that it accurately represents sequenced reads (Table 2). To evaluate assembly completeness, we performed a search for metazoan Benchmarking Universal Single-Copy Orthologs (BUSCOs). Out of 978 metazoan BUSCOs in the OrthoDB catalog, 85.8% were present as complete orthologs in the seal assembly, 51.6% had matches to a single seal transcript, 34.1% were duplicated, 13.8% were fragmented, and only 0.4% were missing from the assembly ( Table 2). This suggests that the blubber assembly is fairly complete as it contains the majority of orthologs predicted to be expressed in all metazoans. Due to the large number of short and potentially redundant transcripts in the www.nature.com/scientificreports www.nature.com/scientificreports/ assembly, we filtered the transcriptome using Transdecoder to retain only transcripts with putative protein-coding regions. The in silico-translated proteome contained 266,916 predicted proteins, some of which were predicted from multiple open reading frames of the same transcripts.
To identify seal homologs of vertebrate genes with known functions, the raw transcriptome assembly was annotated by translated nucleotide search (BLASTX) and the Transdecoder-predicted peptides were annotated by protein search (BLASTP) against the UniProt SwissProt protein database (e-value threshold 1e-3). BLASTX homologs were detected for 372,783 elephant seal transcripts and BLASTP homologs were detected for 164,557 Transdecoder-predicted elephant seal proteins. Transcriptome annotation data are available at figshare (https:// figshare.com/s/37cba07b2e9877f25b50). To identify gene functional categories (KEGG pathways) overrepresented in the blubber transcriptome, we used the DAVID functional annotation tool. Fifty-two KEGG pathways were enriched in the assembly relative to the human genome background (adjusted p < 0.05). The largest category was metabolic pathways, which contained 882 elephant seal homologs. Other pathways of interest included the child categories purine metabolism (135 genes), insulin signaling pathway (112 genes), Wnt signaling pathway (102 genes), AMPK signaling pathway (101 genes), insulin resistance (96 genes), thyroid hormone signaling pathway (94 genes), carbon metabolism (91 genes), adipocytokine signaling pathway (58 genes), fatty acid metabolism (44 genes), and fatty acid degradation (39 genes), among others. These KEGG enrichment data were similar to those from other elephant seal blubber transcriptomes 14,46 . Differential gene expression. Differences in transcript abundance between sampling conditions were identified using DESeq2. The numbers of differentially expressed (DEG) and annotated genes (adjusted p-value < 0.05 and FDR < 0.05) identified in each pairwise comparison between conditions are shown in Table 3. The magnitude and significance of gene expression differences in each comparison are presented in Figs 2 and 3. We first examined changes in gene expression in response to the first ACTH administration, fourth ACTH administration, and over the course of the entire experiment (day 1 pre-ACTH vs. day 4 ACTH response; "overall ACTH response;" Fig. 2). We then compared expression between day 1 and day 4 pre-ACTH conditions ("pre-ACTH comparison") and between day 1 and day 4 ACTH responses ("ACTH response comparison, " Fig. 3).
We identified 61 annotated, unique (i.e. all transcript isoforms collapsed into a single gene) DEGs in response to the first ACTH administration (Supplementary File 1)  www.nature.com/scientificreports www.nature.com/scientificreports/ lists identified 2 Gene Ontology (GO) biological process categories enriched in first ACTH response (positive regulation of fat cell differentiation, negative regulation of transcription from RNA polymerase II promoter), 18 GO categories enriched in overall ACTH response (including fatty acid transport, long-chain fatty acid metabolic process, long-chain fatty acid-CoA ligase activity, positive regulation of fat cell differentiation, response to hypoxia, transcription factor binding) and 1 GO category enriched in the ACTH response comparison (zinc ion binding). KEGG categories were significantly enriched only in the overall ACTH response dataset (ECM-receptor interaction, focal adhesion, pathways in cancer, and small cell lung cancer).
We next manually categorized unique annotated DEGs based on known function using literature search to understand their potential roles in repeated stress responses and metabolic homeostasis. Genes that were upregulated (n = 58) during the first ACTH response included many known GC targets and genes of interest involved in oxidative stress defenses (GPX3), fatty acid oxidation (ACADM), immune function (TLR4), adipogenesis (DKK1), lipid and monocarboxylate transport (SLC16A9), ketogenesis (HMGCS2), insulin resistance and obesity (ZBTB16), NADP biosynthesis and saturated fatty acid metabolism (NADK2), circadian rhythm (PER1), and inactive steroid receptor complexes (FKBP5). Other differentially expressed genes were associated with sphingolipid metabolism, DNA replication and repair, cytoskeleton and ECM remodeling, protein synthesis and posttranslational modification, protein transport, and embryonic development and cell differentiation. Genes downregulated (n = 3) during the first ACTH response were associated with ECM remodeling and insulin resistance (ADAMTS9), cell growth and migration (ABI2), and negative regulation of tyrosine kinase activity (PTPN12) (Supplementary File 1).
Genes upregulated (n = 10) during the fourth ACTH response included those associated with protein degradation (CACYBP), inhibition of apoptosis (DYRK3), and transcriptional regulation (NCOR1), as well as those involved in cell adhesion and migration, sphingolipid biosynthesis, and nucleotide metabolism. Genes downregulated (n = 2) during the fourth ACTH response included an activin receptor associated with adipogenesis Figure 2. Volcano plots showing log2 fold change (x-axis) and significance (−log10 * adjusted p-value; y-axis) of genes differentially expressed during the overall ACTH response (a), first ACTH response (b), and fourth ACTH response (c). Significantly upregulated and downregulated genes are shown in orange and purple, respectively. Log2 fold-change and significance cutoffs for differential expression were |1| and 0.05, respectively. www.nature.com/scientificreports www.nature.com/scientificreports/ (ACVR2A) and a protein glycosylation enzyme (Supplementary File 2). Only two DEGs (ABI2 and SPTSSB) were common between day 1 and day 4 ACTH responses.
Lastly, we examined the potential for interaction between proteins encoded by annotated DEGs by conducting protein-protein interaction (PPI) network analysis using the overall ACTH response dataset. The predicted PPI network contained 20 nodes, 25 edges, network density of 0.132, network heterogeneity of 0.687, network centralization of 0.263, average local clustering coefficient of 0.135, and PPI network enrichment p-value of 0.00386 (Fig. 4), suggesting that DEGs are likely functioning as a coordinated network to regulate cellular and metabolic responses to repeated ACTH administration.

Discussion
Overall, this study shows that gene expression differences between single and repeated ACTH responses in marine mammal blubber are detectable despite similarities in cortisol responses, providing more sensitive indication of stress response states than endocrine profiles alone. We administered ACTH once daily for four days to juvenile northern elephant seals and previously reported significant elevation in corticosteroid levels in response to the first and fourth administrations 12 . However, while the aldosterone response was facilitated by repeated ACTH administration, the cortisol response was not, highlighting the unreliability of glucocorticoid measurements alone for discriminating between acute and repeated HPA axis activation. Therefore, we used RNAseq to examine downstream changes in blubber gene expression in response to repeated ACTH administration. Such changes in blubber gene activity reflect activation of genomic and cell-surface hormone receptors, including glucocorticoid receptor (GR), mineralocorticoid receptor (MR), and the adrenocorticotropic hormone receptor (MC2R) in multiple cell types, capturing the complexity of target tissue responses [47][48][49] . We identified DEGs during the responses to both the first and fourth ACTH administration, although the number of DEGs identified in response to the latter was much smaller than the former. We also identified transcriptional differences between pre-ACTH and post-ACTH response states from the first and fourth days. Annotated DEGs were associated with functions such as lipid and carbohydrate metabolism, adipogenesis, adipokine regulation of energy balance, insulin resistance, oxidative stress, transcriptional regulation, and ECM remodeling in other mammals and provide insights into mechanisms by which repeated HPA axis activation may affect marine mammal physiology.
Genes upregulated in response to the first ACTH challenge were associated with protection from oxidative damage, fatty acid oxidation, immune response, adipogenesis, nutrient transport, sphingolipid and ketone metabolism, cell differentiation, cytoskeletal and ECM remodeling, and insulin resistance and obesity in other animals, while downregulated genes were associated with ECM remodeling and cell migration. Several identified DEGs are primarily known for their roles in other biological processes (e.g. embryonic development) and have only recently been implicated in lipid metabolism and adipogenesis. These include GDF3 and DKK1, which promote adipogenesis and metabolic homeostasis 50,51 , and the transcriptional regulator ZBTB16, which plays a role in obesity and brown adipose energetics 52 . This is consistent with pro-adipogenic effects of GCs in other mammals 24 . Other genes upregulated in response to the first ACTH adminstration included those associated with immune response, which is consistent with studies in other mammals in which acute stress enhanced immune capacity,  www.nature.com/scientificreports www.nature.com/scientificreports/ while chronic stress had immunosuppressive effects 53 . DKK1 and similar functional classes of genes associated with adipogenesis, immune response, and extracellular matrix function were also upregulated in blubber 2 hours after a single ACTH administration 14 . This suggests that transcriptional changes measured 2 and 4 hours after ACTH administration are capturing the same biological response in blubber. Upregulation of genes involved in lipid metabolism is consistent with known lipolytic effects of corticosteroids and the increase in circulating free fatty acid levels that has been measured in response to ACTH in this species 25,28,54 . Interestingly, one upregulated gene, FKBP5, is a GR antagonist 55 , suggesting a potential negative feedback mechanism to limit acute responses to ACTH. Genes downregulated in response to the first ACTH administration were assocated with ECM remodeling and cell migration, suggesting suppression of some energy-intensive functions during the ACTH response. However, many upregulated genes were also involved in ECM remodeling, highlighting the complexity of this specialized connective tissue.
Genes upregulated in response to the fourth ACTH administration included those associated with regulation of transcription, protein degradation, inhibition of apoptosis, and cell migration, while downregulated genes included those associated with adipogenesis and post-translational protein modification. Significantly fewer DEGs were identified in response to the fourth ACTH administration compared to the first, potentially because many cellular responses were already engaged by the previous ACTH responses. Upregulated genes included two that may protect blubber tissue from deleterious effects of repeated corticosteroid exposure: DYRK3, a kinase that phosphorylates and activates the sirtuin SIRT1 and promotes cell survival during stress 56 and NCOR1, a nuclear corepressor that mediates ligand-induced downregulation of the GR gene 57 . DYRK3 was also upregulated in response to a single ACTH administration in a previous study 14 . Surprisingly, the type II BMP receptor ACVR2A, which promotes adipogenesis by GDF6 signaling 58 , was downregulated during the fourth ACTH response, whereas pro-adipogenic factors were upregulated in response to a single ACTH administration in this and a previous study 14 . These results suggest that the transcriptional response to a fourth sequential ACTH administration may serve to promote cell survival and limit adipogenesis and GC sensitivity. However, it is unclear whether this potential negative feedback is sufficient to prevent deleterious consequences of repeated HPA axis activation, as there were dramatic gene expression differences between the baseline state on day 1 and the ACTH response state on day 4.
The largest number of DEGs was identified in the overall ACTH response comparison, as this comparison captured all gene expression differences between pre-ACTH samples on day 1 and post-ACTH samples on day 4. Upregulated genes of interest suggest that repeated ACTH administration impacted lipid and ketone metabolism, adipogenesis, redox homeostasis, and thyroid hormone signaling. Upregulation of genes encoding leptin and the ketogenesis enzyme HMGCS2 suggests a potential increase in energy expenditure and ketone production to support metabolic demands of repeated HPA axis activation. However, leptin is also involved in modulating immune responses 59 , and its role in elephant seal metabolism is currently unclear 60 . Genes associated with lipolysis and fatty acid oxidation were upregulated concomitant with those involved in fatty acid and triglyceride synthesis (e.g. ACSM1). ACSM1 was also upregulated 24 hours after a single ACTH administration in a previous study 14 . Glucorticoid-associated lipogenesis has been described in obese humans 24,61 , and may be mediated by MR 62 . Upregulation of GPD1, which promotes triglyceride synthesis in response to excess glucose, and AGT, a component of the renin-angiotensin (RAS) system associated with obesity, suggest that repeated ACTH administration may potentially increase lipogenesis in a naturally obese and insulin-resistant mammal 30 . Downregulation of the iodothyronine deiodinase DIO2, which converts thyroxine to bioactive triiodothyronine (T3), is consistent www.nature.com/scientificreports www.nature.com/scientificreports/ with suppression of circulating T3 levels after repeated ACTH administration that was shown in the same animals 12 . Therefore, these data suggest that fasting-adapted mammals may simultaneously metabolize and maintain lipid stores during repeated HPA axis perturbation. However, suppression of T3 and increased lipogenesis may decrease energy expenditure and impair the ability of fasting-adapted marine mammals to participate in energetically demanding life history stages.
Our previous studies have shown that while cortisol levels returned to baseline within 24 hours of each ACTH administration 12 , changes in gene expression in response to ACTH could be detected even when cortisol was no longer elevated 14 . To identify whether transcriptional changes alone can be used to detect repeated ACTH exposure, we compared gene expression profiles in blubber samples collected immediately prior to the first and fourth ACTH administrations. Genes upregulated in pre-ACTH samples from day 4 relative to day 1 were associated with regulation of glycolysis, ketogenesis, fatty acid biosynthesis, adipokine signaling, oxidative stress, and adipogenesis. Downregulated genes were associated with genetic and epigenetic regulation of transcription, mitochondrial dynamics, and immune cell differentiation. These results suggest that recovery from repeated ACTH administration may involve metabolic and cellular adjustments such as increased glycolysis and fatty acid oxidation, synthesis of fatty acids and ketones, antioxidant responses, and increased adipogenesis. Upregulation of antioxidant enzymes GPX3 and MGST1 may serve to protect seal tissues from oxidative damage associated with exposure to GCs 63,64 . Two genes upregulated in day 4 pre-ACTH samples (24 hours after the third ACTH administration) -ACSM1 and CDO1 -were also upregulated 24 hours after a single ACTH administration in a previous study 14 . However, genes associated with ketogenesis and antioxidant defenses were not altered in response to a single ACTH challenge and may therefore serve as markers of recently experienced repeated HPA axis activation.
To identify transcriptional differences between blubber responses to single and repeated ACTH administrations, we compared the ACTH response samples from days 1 and 4. We identified nine DEGs that were upregulated and three that were downregulated in day 4 ACTH response samples compared to those from day 1, and suggest them as potential markers of repeated HPA axis stimulation in marine mammal blubber. Upregulated genes were associated with lipid droplet formation, energy homeostasis, and lipolysis, while downregulated genes were associated with adipogenesis and lipid uptake and oxidation. Three upregulated genes encoded the lipid droplet-associated proteins CIDEA and the perilipins PLIN1 and PLIN4, which are involved in regulation of lipolysis and lipid droplet expansion 65 . While GCs have been shown to upregulate genes encoding some lipid droplet proteins 66,67 , they mainly affect perilipins by inhibitory phosphorylation that facilitates lipolysis 68 . Transcriptional upregulation of perilipins in response to repeated ACTH administration may function to restrict lipolysis in an attempt to conserve lipid stores in a fasting-adapted mammal. Other genes upregulated during the fourth ACTH response compared to the first included the glycoprotein AZGP1 and glycine receptor GLRA2, which promote lipolysis 69 , ADAMTS16, which is associated with regulation of blood pressure 70 , PIAS4, which inhibits AMPK and SIRT1 (and thus derepresses lipid and protein synthesis 71,72 ), and CNTNAP3, a cell adhesion protein that is a marker of brown adipocytes 73 . Downregulated genes included THBS1, which is associated with adipocyte proliferation and fatty acid uptake 74 , C1QTNF3, an adipokine that increases hepatic lipid oxidation and decreases lipid synthesis, gluconeogenesis and inflammation 75 , and TGFBI, a TGF beta family member associated with type 2 diabetes and adipogenesis 76 . While interpretation of DEG roles in ACTH responses is complicated by their pleiotropic (or poorly described) functions and the heterogeneity of cell types in blubber tissue, these data suggest that repeated ACTH administration primarily increases lipid catabolism and dysregulates lipid storage. This is consistent with significant declines in circulating triglyceride levels detected after 4 ACTH administrations 77 . Therefore, repeated HPA axis activation may impact the ability of marine mammals to maintain prolonged fasts and high levels of energy expenditure required to participate in key life history stages. For example, depletion of energy reserves by repeated stress would impact reproductive expenditure (e.g. milk production, male-male competition), decreasing fitness 78 .
To our knowledge, this study is the first to examine changes in gene expression in wild marine mammals in response to repeated ACTH administration. Despite the limitations of small sample size, high degree of individual variability in ACTH responses, and small number of ACTH administrations, we identified dozens of DEGs that were differentially expressed in response to repeated ACTH administration. Most significally, we identified gene expression differences between single and repeated ACTH responses in blubber despite similarities in cortisol responses, providing more sensitive indication of response states than endocrine profiles alone. Further work is necessary to correlate transcriptional changes with metabolic and redox profiles and validate the ability of these biomarkers to discriminate baseline stress states in marine mammals. Nevertheless, these findings provide some of the first insights into the physiological impacts of repeated HPA axis stimulation in free-ranging marine mammals, and suggest potential molecular markers of repeated stress exposure that can complement other approaches used by conservation biologists to evaluate the effects of anthropogenic activity on threatened wildlife populations.