A blubber gene expression index for evaluating stress in marine mammals

Blood hormone levels are often used to measure stress in animals, but blood cannot be sampled from many large, wild marine mammals. Gene expression in blubber, which can be sampled more easily, is correlated with blood cortisol and thyroid hormones, providing an alternative way to measure stress in marine mammals.


Introduction
Many marine mammal populations are in decline, which may make them especially vulnerable to anthropogenic disturbance (Davidson et al., 2012). Human activity in the world's oceans continues to increase and is associated with multiple stressors (e.g. sound, chemical and physical pollution, fishing, climate change), which may have complex interactive effects on marine wildlife (Halpern et al., 2015, Maxwell et al., 2013, Simmonds, 2018. The mammalian physiological response to stress is mediated by the sympathetic nervous system via catecholamines, and the hypothalamic-pituitary-adrenal (HPA) axis via glucocorticoids (GCs). While adaptive in the short term, repeated or chronic stress exposure has been associated with adverse effects on reproduction, immune function, metabolism and oxidative load, among other functions (Chatelain et al., 2020, Costantini et al., 2011, Romero et al., 2016, Sapolsky et al., 2000. In marine mammals, disturbance may also cause alterations in diving and foraging behavior, potentially influencing foraging effort and success (Erbe et al., 2018, Nabi et al., 2018. Disturbance mitigation efforts require an understanding of the transfer functions that link disturbance to physiology and the ability to identify affected animals. Therefore, recent research priorities have focused on developing methods for assessing the physiological impacts of stressors on marine mammals (National Academies of Sciences, 2017).
Evaluating the incidence and impacts of stress in wildlife is notoriously challenging (Dickens et al., 2013. Commonly used markers, such as circulating GCs, are extremely variable and do not fully capture the complexity of the physiological response (MacDougall-Shackleton et al., 2019. For example, marine mammal responses to various stressors can be associated with significant elevation of the mineralocorticoid aldosterone and decrease in thyroid hormone levels , DeRango et al., 2019a, DeRango et al., 2019b, McCormley et al., 2018, St Aubin et al., 1992. Furthermore, circulating hormone measurements are not practical for many marine mammal species because they are difficult to access for blood sampling (e.g. large size, aquatic lifestyle) and, like other wildlife, are subject to the artefacts of capture stress. To address these challenges, methods for remote biopsy sampling (e.g. via a dart attached to a pole or deployed by crossbow) of marine mammals and hormone detection in a variety of other matrices, such as blubber, respiratory vapor and feces, have been actively developed (Boggs et al., 2017, Burgess et al., 2017, Burgess et al., 2018, Hunt et al., 2019, Noren et al., 2012. However, hormone measurements often lack the sensitivity to identify recent exposure to a stressor and do not provide additional information on downstream consequences. Corticosteroid and thyroid hormones exert their physiological effects by binding to receptors in target tissues and altering the activity of cellular proteins and genes (Brent, 2012, Sapolsky et al., 2000. Glucocorticoid (GR), mineralocorticoid (MR) and thyroid (TR) hormone receptors regulate expression of hundreds of genes involved in metabolism, inflammation and endocrine signaling in mammalian adipose, a tissue that is abundant and accessible in many marine mammals (Lee et al., 2014, Obregon, 2014. Changes in expression of stress-responsive genes in blubber, the specialized subcutaneous adipose tissue in marine mammals, have the potential to serve as a highly sensitive indicator of recent (within 24 h; Khudyakov et al., 2017) exposure to a stressor. As marine mammal blubber is vertically stratified by fatty acid composition, cellular morphology, rate of lipid accumulation and mobilization and gene expression (Ball et al., 2017, Guerrero et al., 2017, Struntz et al., 2004, the impacts of stress hormones are likely to be more pronounced in the inner, more metabolically active layer of blubber. We previously identified a cluster of genes that were differentially expressed in inner blubber in response to acute and repeated HPA axis stimulation in a model marine mammal, the northern elephant seal (Mirounga angustirostris; Deyarmin et al., 2019, McCormley et al., 2018, Khudyakov et al., 2017. These included upregulated genes encoding metabolic and antioxidant enzymes (long-chainfatty-acid CoA ligase 1, ACSL1; hydroxymethylglutaryl-CoA synthase 2, HMGCS2; cysteine dioxygenase 1, CDO1; E3-SUMO-protein ligase PIAS4, PIAS4; glutathione peroxidase 3, GPX3; microsomal glutathione S-transferase 1, MGST1), adipokines (adiponectin, ADIPOQ; leptin, LEP), proadipogenic factors (dickkopf-related protein 1, DKK1; zincalpha-2-glycoprotein 1, AZGP1), lipid droplet proteins (perilipin 1, PLIN1; cell death activator CIDE-A, CIDEA) and other metabolism-regulating proteins (glycine receptor subunit alpha-2, GLRA2), and a gene associated with extracellular matrix remodeling (transforming growth factor beta-induced, TGFBI), which was downregulated (Deyarmin et al., 2019). Gene expression changes were correlated with other markers, including increase in cortisol, aldosterone and reverse triiodothyronine (rT3) levels and decrease in total triiodothyronine (tT3) levels (McCormley et al., 2018).
Many of the hormones that are responsive to stressors serve important functions in the maintenance of metabolism and show natural variations associated with daily and seasonal activities in animals. Genetic markers that co-vary with these hormones under normal (e.g. predictable and/or cyclical) activity may also be useful indicators of the response to a stressor. In this study, we evaluated the potential of these gene markers to discriminate hormonal conditions in marine mammals without the addition of exogenous stressors (i.e. natural variation within the population). We measured their expression in inner blubber along with circulating hormone levels in a cross-section of juvenile northern elephant seals of varying physiological states and identified a gene expression index comprised of 10 genes that varied with cortisol and thyroid hormone levels. We propose that baseline inner blubber expression levels of these genes may serve as sensitive and informative indicators of potential stress exposure in marine  mammals and that the approach described here may be used to identify markers in other marine mammal tissues, such as outer blubber and skin.

Study subjects
All animal handling procedures were conducted under National Marine Fisheries Service permit 19108 and were approved by Sonoma State University and University of the Pacific Institutional Animal Care and Use Committees and the Department of the Navy's Bureau of Medicine and Surgery. A total of 30 juvenile (0.8-1.8 year old) northern elephant seals of varying apparent body condition were sampled for the study (17 females, 12 male, one unknown). One male subject was excluded from the analyses due to an insufficient amount of sample for gene expression analyses.

Sample collection
Baseline sample collection was conducted as previously described (Deyarmin et al., 2019, McCormley et al., 2018 during September to November of 2017 at Año Nuevo State Reserve, San Mateo County, CA, USA. This period is generally a short haul out and is not associated with extensive fasting, molting or reproduction. Seals were chemically immobilized using an intramuscular injection of ∼1 mg/kg tiletaminezolazepam HCl (Telazol, Fort Dodge Animal Health, USA). Sedation was maintained with intravenous doses of ketamine (0.25-1 mg/kg) (Fort Dodge Animal Health, IA, USA). Blood samples were collected within 15 ± 5 min of initial immobilization from the extradural vein using an 18 G, 3.25inch spinal needle into vacutainer tubes (BD Life Sciences, USA) and chilled until further processing. Serum or plasma were isolated by centrifugation for 15 min at 3000 g and stored at −80 • C. Blubber samples were collected within approximately 30 min of initial immobilization from the posterior flank of the animal using a 6.0-mm diameter biopsy punch (Miltex, USA). Two blubber biopsies were collected from each animal, one of which was immediately frozen in liquid nitrogen while the other was minced and incubated in RNAlater TM Stabilization Solution (∼300 mg tissue per 1.5 mL; Invitrogen, USA) with tube rocking for 24 h at 4 • C, which has been shown to decrease RNA degradation in adipose tissue samples from other mammals (Hemmrich et al., 2010). Blubber samples were stored at −80 • C until further processing (RNAlater TM was removed from samples prior to freezing). After collection of morphometric measurements of standard length (SL) and axillary girth (AG), seals received rear flipper tags (Dalton, Oxon, UK) and were allowed to recover from anesthesia and resume normal activity.

RNA isolation
The inner half of each blubber sample (closest to muscle) was used for RNA extraction as this tissue layer is considered more metabolically active than outer blubber in marine mammals (Ball et al., 2017, Struntz et al., 2004. In addition, inner blubber was used for elephant seal transcriptomes that identified the markers used in this study (Deyarmin et al., 2019, Khudyakov et al., 2017. Blubber tissue (∼100 mg) was minced with a sterile scalpel on dry ice and homogenized in 1 mL of Qiazol (Qiagen, USA) by bead beating in a Bullet Blender Storm 24 instrument (Next Advance, USA) for two 2min cycles. Tissue homogenates were further disrupted using QIAshredder tubes (Qiagen, USA). RNA was purified using the RNeasy Lipid Mini Kit (Qiagen, USA) following the manufacturer's protocol with a 15-min on-column DNase I digest. RNA quantity and integrity were determined using the Broad Range RNA Assay on the Qubit 3.0 Fluorometer (Life Technologies, USA) and the RNA Pico 6000 Assay on the 2100 Bioanalyzer (Agilent Technologies, USA), respectively. RNA samples had integrity values (RINs) between 7.0 and 8.6. No differences in RINs were observed between samples that were flash-frozen in liquid nitrogen and those that were preserved in RNAlater TM solution. Therefore, tissues preserved using either method were used interchangeably for RNA isolation.

RT-qPCR
Complementary DNA (cDNA) was synthesized from 500 ng of total RNA by reverse transcription (RT) using SuperScript IV VILO Master Mix with ezDNase (Thermo Fisher, USA). cDNA samples were diluted 1:10 and 2 μL were used in each 20 μL real-time polymerase chain reaction (qPCR) with PowerUp SYBR Green Master Mix (Thermo Fisher, USA). qPCR was performed on a QuantStudio 5 Real-Time PCR System instrument (Thermo Fisher, USA) using the following program: 2 min at 50 • C, 2 min at 95 • C and 40 cycles of 15 s at 95 • C followed by 60 s at 60 • C. All samples were run in triplicate with all intra-assay and inter-assay CVs < 0.5%. No-template and no-RT enzyme controls were included in each run and showed no amplification.
Primers for qPCR assays were designed to target genes that were previously identified as differentially expressed in elephant seal blubber in response to repeated ACTH administration (Deyarmin et al., 2019). Candidate genes were selected based on three characteristics: (i) significant BlastX hits to a protein with known function in the UniProt SwissProt database, (ii) high transcript abundance (number of transcripts per million, TPM ≥ 20) and (iii) relatively consistent expression levels between biological replicates (mean CV for replicate TPM = 24.4%). Primers were designed to target highly conserved regions of differentially expressed genes using NCBI Primer-Blast. Primer sequences are shown in Table 2. All primers were used at 400 nM final concentration. Primer efficiencies (Table 2) were determined using standard curves of five 1:2 dilutions of pooled cDNA. Primer specificity was confirmed using melt curve analysis, gel electrophoresis and Sanger sequencing of qPCR products.
YWHAZ and NONO, which were used as reference genes in a previous blubber gene expression study (Khudyakov et al., 2017), were evaluated for stability in expression across samples using BestKeeper algorithm (Pfaffl, 2004). YWHAZ was the most stable gene across all samples (SD = 0.47, CV = 2.06%) and was thus used as the reference gene in this study. Normalized gene expression values (delta C T ) were calculated by subtracting the C T of the gene of interest from the C T of the reference gene (Huggett et al., 2005, Schmittgen et al., 2008. Normalized gene expression ratios and primer efficiencies were calculated using Microsoft Excel 2017.

Statistical analyses
All statistical analyses were conducted using R v3.6.0 in RStudio v1.1.453 (R Core Team, 2016. Association between variables was evaluated using Spearman correlation (r S ) analysis with Benjamini and Hochberg correction for multiple hypothesis testing (FDR = 0.05) using the 'psych' package (Revelle, 2019). Differences in putative stress markers between sexes were evaluated by one-way analysis of variance (ANOVA). Levene's and Shapiro-Wilk tests were used to determine whether variables met equal variances and normality assumptions, respectively. Nonparametric ANOVA (Kruskal-Wallis χ 2 test) was used for variables that did not meet ANOVA assumptions.
Body condition index (CI) was calculated as CI = (AG/SL) × 100, which has been correlated with adiposity in several species of phocid seals (Pitcher, 1986, Rea et al., 2016, Ryg et al., 1990. Large CI values indicate animals with higher AG, and thus potentially higher adiposity. Therefore, animals with higher CI are likely to be healthy, while those with lower CI may be experiencing nutritional stress (Le Boeuf and Crocker, 2005).
Principal components analysis (PCA) was conducted using the principal function in the 'psych' package (Revelle, 2019). Initial analysis was used to obtain eigenvectors for each component in the data. Three components had eigenvalues above Kaiser's criterion of 1 and in combination explained 76% of the variance in the data; these were extracted for the final analyses. The mean of the communalities for variables after extraction was 0.76, suggesting that Kaiser's criterion was appropriate for this dataset (communalities >0.7 for analyses with <30 samples; Field et al., 2012). Factor interpretability was improved using varimax rotation. Rotated factor loadings are presented in Table 3.
Rotated factor scores were used for linear regression as predictors of hormone levels. All three components were initially included in the models, after which components that did not have significant effects were removed and more parsimonious models were run. Model residuals were visually assessed for normality and homoscedasticity and variables were logtransformed if necessary. Specifically, cortisol and aldosterone were the only variables that required log transformation. All final models were assessed for presence of outliers and influential cases, which were not detected.
PCA was used to reduce dimensionality of the gene expression dataset. This approach was deemed appropriate using several diagnostic criteria (Field et al., 2012). First, examination of the correlation matrix of dC T values for 14 genes identified four (MGST1, GLRA2, PIAS4, LEP) that were removed from subsequent analyses due to having two or fewer significant (P < 0.05) correlations with other gene markers (Fig. 2). Bartlett's test of sphericity indicated that the correlations between the remaining 10 genes were suf-   Table 2: Primer sequences and amplification efficiencies (E) of qPCR assays. All sequences are written in the 5 to 3 direction. Primers were designed using gene sequences from elephant seal blubber transcriptomes (Deyarmin et al., 2019, Khudyakov et al., 2017 (Table 3).
There were no significant associations between the components and aldosterone (P = 0.25) or CI (P = 0.36). Aldosterone was not associated with any of the gene markers, while CI was negatively associated with DKK1 (r S = −0.49, P = 0.036). RC3 varied with sex and was higher in males than females (ANOVA: F 2,26 = 9.47, P = 0.0008). RC1 and RC2 were not associated with sex (ANOVA; RC1: P = 0.69, RC2: P = 0.31).

Discussion
In this study, we evaluated whether circulating levels of aldosterone, tT3 and rT3 and expression of a suite of previously identified stress-responsive genes (Deyarmin et al., 2019) may be more informative than baseline cortisol levels alone in assessing stress state in marine mammals. We measured levels of these markers in a cross-sectional sample of elephant seals of differing physiological condition evidenced by variability in body condition and hormone levels. These stress-responsive genes included several known targets of glucocorticoid and thyroid hormones and encoded adipokines, lipid metabolism enzymes and factors that regulate adipogenesis. We found that that baseline blubber expression of 10 genes-ADIPOQ, ACSL1, HMGCS2, AZGP1, PLIN1, CIDEA, GPX3, CDO1, DKK1 and TGFBI-was predictive of circulating levels of cortisol and thyroid hormones but not aldosterone or body condition. Gene expression data were collapsed into three principal components that significantly varied with circulating hormone levels. We propose that RC1 represents an index of GR (and potentially TR) target gene expression that reflects increased cortisol and decreased tT3 levels and potential consequent alterations in adipogenesis, lipid metabolism, insulin sensitivity and redox homeostasis.   index of pro-adipogenic gene expression that is associated with elevated rT3. Finally, RC3 represents an index of gene expression associated with reduced tT3 levels and consequent potential alterations in lipid homeostasis in marine mammals.
RC1, which represented expression of the antioxidant enzyme GPX3 and pro-adipogenic factor DKK1, was positively associated with cortisol and negatively associated with tT3. Cortisol induces differentiation of pro-adipocytes into mature adipocytes, stimulates lipolysis and promotes insulin resistance (Lee et al., 2014), while T3 stimulates proliferation and inhibits differentiation of pro-adipocytes and enhances both lipogenesis and lipolysis in mammals (Carmean et al., 2014). Suppression of thyroid hormone levels during stress and a negative relationship between cortisol and T3 have been reported in several terrestrial and marine mammal species (Atkinson et al., 2015, Helmreich et al., 2011. Expression of GPX3 and DKK1, which are both known GR target genes, is associated with adipogenesis in rodents and humans (Christodoulides et al., 2006, Lee et al., 2008. GPX3 is also associated with insulin sensitivity and protection from oxidative damage and inflammation in adipose tissue (Christodoulides et al., 2006). A recent study found that DKK1 was upregulated in human adipocytes expressing a mutant form of thyroid receptor (TRβ), suggesting that it may be directly inhibited by T3 signaling (Liu et al., 2015). Upregulation of RC1, in conjunction with elevated cortisol and decreased tT3, may promote maturation of adipocytes, potentially in order to sustain elevated rates of lipid metabolism during stress.
RC2, representing expression of CIDEA, CDO1 and TGFBI, was positively associated with rT3, a product of thyroxine (T4) deiodination by type III deiodinase (D3) that is thought to be biologically inactive. Reverse T3 was elevated and correlated with cortisol during responses to capture stress in belugas and ACTH administration in elephant seals (Champagne et al., 2015, Ensminger et al., 2014, St Aubin et al., 1992. Interestingly, D3 was upregulated in proliferating rat adipocytes, suggesting that rT3 production may also play a role in adipogenesis (Carmean et al., 2014). CIDEA encodes a lipid droplet protein that is associated with insulin sensitivity and adipose expansion in humans (Abreu-Vieira et al., 2015). CDO1 encodes an enzyme in the taurine biosynthesis and cysteine degradation pathways that has been implicated in adipogenesis (Deng et al., 2015). TGFBI encodes an extracellular matrix protein that was recently identified as a negative regulator of adipogenesis in human cells (Zhang et al., 2019). TGFBI was downregulated in response to repeated HPA axis stimulation in elephant seals (Deyarmin et al., 2019) and was negatively associated with rT3 in this experiment. Another gene, HMGCS2, which encodes a ketogenic enzyme that regulates fatty acid oxidation in a variety of cell types (Sikder et al., 2018), loaded onto both RC2 and RC1. It is possible that the association between RC2 gene expression and rT3 levels is indirect, being driven primarily by the coincident decrease in tT3 and increase in cortisol. However, RC2 explained only 7.8% of the variance in tT3 levels and HMGCS2 was the only gene that loaded onto both cortisol-associated RC1 and rT3associated RC2. Recent research suggests that rT3 may bind to non-genomic TR isoforms and alter cytoskeletal dynamics (Davis et al., 2008); it is possible that it may have other biological functions. An increase in RC2 expression and rT3 levels may serve to enhance adipogenesis and fat catabolism during stress.

RC3
, which encompassed expression of AZGP1, ACSL1 and PLIN1, was negatively associated with tT3. A decrease in levels of tT3, which is associated with lipogenesis (Mullur et al., 2014), coincident with increased cortisol levels suggests an enhanced capacity for lipid catabolism during stress. AZGP1 promotes lipolysis and inhibits lipogenesis in mice and humans (Bing et al., 2010). ACSL1 is an enzyme involved in fatty acid uptake, beta-oxidation and reesterification (Ellis et al., 2010). PLIN1 regulates fat droplet homeostasis and facilitates hormone-stimulated lipolysis (Kimmel et al., 2016). ADIPOQ, which loaded onto both RC1 and RC3, is associated with insulin sensitivity, adipose tissue expansion and fatty acid oxidation (Stern et al., 2016). Expression of PLIN1 and ADIPOQ is indirectly regulated by thyroid hormones (Liu et al., 2012, Liu et al., 2015. Upregulation of RC3 and decrease in tT3 may potentially contribute to increased rates of lipolysis and fatty acid oxidation during stress.
The association of RC3 with sex was surprising given that tT3 levels measured in this study did not vary by sex. However, previous studies of this age class of elephant seals found that T3 and T4 were higher in males compared to females and were associated with differences in field metabolic rate between sexes (Jelincic et al., 2017, Kelso et al., 2012. It is possible that the small sample size and variability in physiological state among individuals used in this study masked observable sex differences in hormone levels. However, RC3 values were higher in males in this study and were significantly associated with tT3, suggesting that blubber gene expression measurements, which are more sensitive than circulating hormone measurements, may reflect true underlying sex differences in thyroid hormone production and signaling in marine mammals. We found no association between gene expression and aldosterone in this study. Aldosterone, a secretagogue of ACTH, was one of the most significantly altered hormones in the repeated ACTH administration experiment (McCormley et al., 2018) and has been shown to be elevated in response to external stressors in several marine mammal species besides elephant seals , DeRango et al., 2019a, DeRango et al., 2019b. Aldosterone is primarily known for its osmoregulatory and vascular functions, but recent studies have suggested that it may have additional roles in adipogenesis and insulin sensitivity (Marzolla et al., 2012). However, the biological effects of aldosterone are mediated by MR, which also binds GCs with high affinity, suggesting that MR signaling in adipose tissue may be mediated in part by GCs (Gomez-Sanchez et al., 2014). Moreover, mean baseline aldosterone concentrations measured in this and other elephant seal studies were 200-to 300-fold lower than cortisol concentrations (McCormley et al., 2018). Therefore, any changes in blubber gene expression that may reflect aldosterone signaling are likely masked by the effects of cortisol.
Body CI was negatively associated with cortisol and rT3, which is consistent with a previous study of juvenile elephant seals that found that T4 levels significantly decreased over a fasting period during which animals lost 17% of their body fat (Kelso et al., 2012). Cortisol was also negatively associated with adiposity in lactating adult female elephant seals (Fowler et al., 2018), but not in juveniles (Kelso et al., 2012), highlighting the importance of evaluating stress state in the context of life history stage, especially in fasting-adapted marine mammals. Body condition was also negatively associated with DKK1, but not with any of the other gene markers, suggesting that this gene may be a novel molecular marker of adiposity in marine mammals. Serum levels of DKK1 were also correlated with adiposity in diabetic humans, although the relationship was positive (Ali et al., 2019). This difference may be due to a mismatch between mRNA expression in adipose and circulating protein levels, or the fact that high adiposity and insulin resistance are considered pathological in humans, but healthy in fasting-adapted mammals (Houser et al., 2013). While CI is an approximation rather than a direct measurement of adiposity, our data suggest a potentially informative relationship between this metric and stress-associated hormones and one of their target genes. More accurate determination of adiposity (e.g. using the isotopic dilution method) will be necessary to validate these associations.
In summary, we have shown that baseline expression of 10 genes, which can be measured in marine mammal blubber by RT-qPCR, was strongly associated with other markers of stress in elephant seals. These genes encode proteins that regulate adipogenesis, fatty acid oxidation, redox homeostasis and insulin sensitivity, providing information about the functional consequences of stress-induced hormone changes in fasting-adapted marine mammals. Such consequences may include alterations in lipid homeostasis that disrupt the ability of capital breeding marine mammals to complete key life history stages that require fasting. Significantly, baseline gene expression levels in inner blubber predicted circulating cortisol and thyroid hormone levels, suggesting that these genes may be sensitive indicators of stress exposure in marine mammal species that are inaccessible for blood hormone measurements. Furthermore, the approach described here can be used to identify species-specific stress markers in other matrices, such as skin and outer blubber, which are more commonly sampled by remote biopsy in many free-ranging marine mammals.