A novel RT‐qPCR health assay reveals differential expression of stress and immunoregulatory genes between the seasonal migrations of humpback whales (Megaptera novaeangliae)

Health information is essential for the conservation management of whale species. However, assessing the health of free‐ranging whales is challenging as samples are primarily limited to skin and blubber tissue. Reverse transcription quantitative polymerase chain reaction (RT‐qPCR) offers a method to measure health from blubber RNA, providing insights into energetic status, stress and immune activity. To identify changes in health, natural differences in baseline gene expression linked to an individual's sex, reproductive status and life‐history stage must first be quantified. This study aimed to establish baseline gene expression indices of health in migrating humpback whales (Megaptera novaeangliae). To do this, we developed an assay to quantify seven health‐related gene transcripts (Leptin, Leptin Receptor, Adiponectin, Aryl Hydrocarbon Receptor, Tumour Necrosis Factor‐α, Interleukin‐6, Heat Shock Protein‐70) and used Bayesian mixed effect models to assess differential baseline expression based on sex, lactation status and migration stage (northbound to and southbound from the annual breeding grounds). Results showed no significant contribution of sex to differential baseline expression. However, lactating individuals exhibited downregulated AhR and HSP‐70 compared to non‐lactating conspecifics. Additionally, southbound individuals demonstrated significantly upregulated HSP‐70 and downregulated TNF‐alpha, suggesting a relationship between these inflammation‐linked transcripts and migratory fasting. Our results suggest that baseline differences due to migratory stage and lactation status should be considered in health applications of this assay. Future monitoring efforts can use our baseline measurements to better understand how gene expression is tied to population‐level impacts, such as reduced prey availability or migratory stressors.

Stranded animals can provide health information from impacted individuals in a population; however, free-ranging controls allow researchers to parse the differences between baseline and pathological health indices.For an adequate comparison, data from healthy individuals must also account for factors such as an individual's reproductive status, sex, age and natural temporal cycles of health biomarkers, given the potential impact of these features on an individual's health outcome.
The large body size and aquatic obligation of baleen whales make it difficult to handle individuals and collect health information.In other animals, health indices are often developed from samples such as blood (including serum), muscle, urine, hair, teeth or faeces (e.g.Kersey & Dehnhard, 2014;Sanchez et al., 2018;Weltring et al., 2012;Young et al., 2004).Some of these methodologies have been applied to small marine mammals, using techniques such as caudal blood draws (e.g.Houser et al., 2021;Kunisue et al., 2011;Lehnert et al., 2017).
However, replicating these studies in large cetaceans is not currently possible due to the logistical challenges of restraining large baleen whales for sample collection procedures.Additionally, the involuntary collection of muscle or visceral tissue is not optimal for free-ranging whales as these require extremely invasive and dangerous procedures.Thus, the primary sources of in vivo samples from large whales are skin and blubber.The exterior of a whale consists of thin epidermal and dermal layers and a thick layer of subcutaneous blubber underneath.The blubber layer is comprised of a complex of connective and adipose tissue that is specialized for thermoregulation and energy storage in marine environments (Iverson, 2009).Although blubber is vascularized, it is unlikely to be a viable source of blood or serum as there are only small volumes present in the tissue.Thus, health information from blubber is likely to be primarily from the adipose tissue itself.Typically, skin and blubber are collected using remote biopsies, where a projectile cutting device is shot at the whale and used to take a small sample from the whale's flank.This method is safe for both the researcher and the animal, and provides access to a range of genetic and physiological information.
Between the two tissue types collected, blubber is particularly useful for health information due to its multifunctionality in thermoregulation, energy storage and as an endocrine organ (Ball et al., 2017;Iverson, 2009;Pujade Busqueta et al., 2020).Blubber also acts as a reservoir for lipid-soluble molecules and has been used to quantify hormones related to stress (Kellar et al., 2017;Mingramm et al., 2020), reproduction (Pallin et al., 2018;Vu et al., 2015) and to monitor pollutant consumption (Gauthier et al., 1997).Remote biopsies have also been used to some extent to determine animals' energy reserves through features such as adipose cell size (Castrillon et al., 2017).
However, measures of lipid content and adiposity have also shown variability that does not reflect morphological body condition or other health features, meaning that these measures may not provide reliable information about energy reserves without additional body condition indices (Christiansen et al., 2020;Kershaw et al., 2019).
Another source of health information is gene expression.This can be measured as the presence of tissue-specific messenger RNA (mRNA) using quantitative reverse transcription polymerase chain reaction (RT-qPCR).Quantitative PCR offers significant advantages for wildlife health studies, including high sensitivity, and the ability to measure a panel of genes simultaneously (Dale et al., 2016;Wittwer et al., 2001).Some studies have noted problems with RNA extraction from blubber tissue, which are likely due to its extremely high lipid content (even for adipose tissue) and low RNA abundance (Castrillon Posada, 2019;Van Dolah et al., 2015).However, recent advances have improved RNA extraction efficiency from whale blubber tissue so that RNA can be used to gather biological information (e.g. on individual sex and energetic status; Ball et al., 2017;Linsky et al., 2022).By expanding these methods to better understand adipose-specific cytokines (adipokines; i.e. cell signalling proteins that are translated from the measured 'blubber RNA'), cetacean health studies may be able to use RT-qPCR as a powerful diagnostic tool that is capable of identifying deviant energetic status, immune activity and stress responses (Fasshauer & Bluher, 2015).In other fasting marine mammals (Mirounga angustirostris), blubber gene expression responds to hormone administration within 2 h and returns to normal in approximately 24 h (Khudyakov et al., 2017).This means that gene expression may be more temporally precise than steroid hormone measurements, which show slow turnover in human adipose and most likely take longer in cetacean outer blubber considering the tissue's low circulation (Hughes et al., 2010).Before gene expression can be used to monitor health features, however, specific health biomarkers (i.e. mRNA transcripts) must first be identified and their baseline expression levels established.
Previous research provides a basis for the selection of specific transcripts as potential cetacean health indices.In terrestrial mammals, adipokines such as leptin, insulin and adiponectin are known to regulate appetite and energy storage (Jequier, 2002;Niswender & Schwartz, 2003).These hormones have also been found to reflect morphological body condition, with leptin levels, in particular, being correlated with the amount of adipose tissue within an individual (Buff et al., 2002;Gentry et al., 2002;Ishioka et al., 2007;Niswender & Schwartz, 2003).Among mammals, migratory baleen whales may have some of the most extreme shifts in leptin expression.For example, the absolute expression levels of blubber leptin in bowhead whales (Balaena mysticetus) at the beginning of the fasting migratory period are 10-fold the adipose leptin expression of any terrestrial mammal entering hibernation (Ball et al., 2017).However, these bowhead whale measurements used post-mortem samples from deep in the blubber layer.For health and energetic monitoring efforts, adipokines like leptin must be measurable (and biologically informative) from the outer blubber layer of biopsies collected from living whales.
Other health features, including stress and immune response, have been linked to in vivo blubber transcripts for several marine mammal species.In killer whales (Orcinus orca), interleukin 10 (IL-10), aryl hydrocarbon receptor (AhR) and several other transcripts show increased expression in response to polychlorinated biphenyl pollutants (PCB, Buckman et al., 2011).This increase in expression suggests they have the potential to be used as indices of immune response to the consumption of organic pollutants.Similar cytokines, such as interleukin-6, are also known to increase in cetaceans in response to other persistent organic pollutants (Funke et al., 2003).Other forms of dietary stress, such as intensive fasting periods, have been linked to immune and inflammatory expression pathways.Examples include differences in the expression of heat shock proteins in pinniped species such as the northern elephant seal (Mirounga angustirostris) and the grey seal (Halichoerus grypus; Bennett et al., 2014;Khudyakov et al., 2017) during calorically intensive moulting and breeding periods.However, the link between a stressor and observed gene expression is not always clear.In terrestrial animals, the expression of both inflammatory transcripts (including tumour necrosis factor-alpha and interleukin proteins) and anti-inflammatory transcripts (e.g.heat shock proteins) have been linked to similar stressors, leading to uncertainty about which cytokine is responding to stress versus regulating this response (Jacquier-Sarlin et al., 1994;Meng & Harken, 2002).Similar expression profiles may represent differential responses depending on the species and the stressor.Thus, it is critical to identify the species-specific response that a transcript may provide before they can be used as monitoring and diagnostic tools for a particular species.
For humpback whales (Megaptera novaeangliae), there is a significant knowledge gap in our understanding of health-related gene transcripts, and how their life-history strategy influences baseline gene expression.However, health monitoring of Southern Ocean humpback whales is desirable for humpback conservation and the conservation of broader Antarctic ecosystems (Bengtson Nash et al., 2018).As long-lived capital breeders, humpback whales amass energy during hyperphagic summers at high latitudes by feeding on Antarctic krill (Euphausia superba).They then expend this energy migrating to equatorial breeding grounds during the winter.This dramatic, but naturally occurring, fasting event has the potential to influence gene expression similar to the seasonal expression changes that occur in other species (e.g.Thompson et al., 2012).
Furthermore, the expression of certain genes is known to vary intraspecifically based on features such as age, sex and reproductive status (Camolotto et al., 2010;Grath & Parsch, 2016;Gui et al., 2004).
For humpback whales, differences in migratory demands between males and females may also be a driver for significant differences in baseline gene expression.For example, lactating females face some of the highest energetic demands, due to calorically intensive milk production (Braithwaite et al., 2015).This could affect what a 'healthy' gene expression profile looks like for these females compared to other reproductive classes.To effectively use gene expression as a health monitoring tool, these natural sources of variation must be quantified at baseline levels.
Eastern Australian humpback whales (EAHW) are an ideal population for such baseline studies.The EAHW has shown some of the strongest growth rates of any population recovering in the post-whaling era (Noad et al., 2007;Paterson et al., 1994).Thus, this population can be considered a 'healthy' source of information to compare to other less-successful whale populations.Additionally, the EAHW is expected to reach carrying capacity sometime between 2021 and 2026 (Noad et al., 2019).This means that a future decline in the abundance of EAHW can be expected.However, baseline expression measurements collected before carrying capacity was reached can be used to monitor whale health during this decline and rapidly identify impacts that may lead to a population collapse.
In this study, we quantified the baseline expression of a panel of health-related genes in migrating humpback whales.To do this, we developed an assay that is capable of concurrently measuring several transcripts (i.e.all primer sets can be run in a single RT-qPCR reaction).This study aimed to (1) select a range of adipokine genes, the transcripts of which are considered to be candidate health indices for humpback whales based on the available information for terrestrial and marine mammals; (2) apply RT-qPCR to measure expression levels of these genes in blubber samples from free-ranging, migrating humpback whales; and (3) assess differences in the expression of the selected adipokines due to sex (male or female), lactation status (lactating and non-lactating) and migratory stage (northbound whales migrating from the feeding grounds to breeding grounds, and southbound whales returning to the feeding grounds).
By quantifying these natural sources of intraspecific variation, the health-related variation in gene expression can be identified and validated against other known indices (e.g.steroid hormones, body condition).This will allow the assay developed in this study to be utilized as a powerful conservation tool that is capable of simultaneously measuring multiple indices of health.With adequate sample sizes, the baseline information provided in this study can be compared to future studies to detect changes in the health of the EAHW population as well as of other humpback whale populations.With minor modifications, this assay may also provide health information for other data-deficient mysticete species.

| Sample collection
Blubber tissue biopsies were collected from free-ranging humpback whales off the coast of North Stradbroke Island (Minjerribah) and Peregian Beach, Queensland, Australia.Samples were collected during the 2020-2021 northbound and southbound annual migrations (June to November) to and from their breeding grounds in the Great Barrier Reef respectively.A total of 93 blubber biopsies were analysed including 57 male and 36 female samples.Of the females, 13 were observed to be accompanied by a calf and were assumed to be lactating, 42 were collected from northbound individuals, and 51 were collected from southbound individuals.
Biopsies were collected using the Paxarm® MK24C remote biopsy system.Typically, a small (~40 mm depth by 7 mm diameter) sample of skin and blubber was taken from the flank of each whale.
The sample was collected from the ocean surface as quickly as possible (about 1-3 min from sampling) and then a small disk (approx.0.03 g) was cut from the deepest point of the sample that remained inside of the cutting head of the biopsy dart when floating.This sample was stored in 1 mL of Qiagen RNAprotect (product number: All samples were accompanied by notes on the biopsied animal's behaviour (e.g.swimming direction) and calf presence/absence (as an additional method to identify females), and by identification photographs of the dorsal fin and underside of the fluke, when possible, to prevent multiple biopsies of the same animal.The sex of each individual was determined post-sampling using a blubber RNA-based X-chromosome inactivation method and cross-validated using ZFX/ SRY amplification from skin tissue (Jayasankar et al., 2008;Linsky et al., 2022).

| RNA extraction and cDNA synthesis
RNA and complementary DNA (cDNA) were obtained using the optimized protocol presented in Linsky et al. (2022).Briefly, the Qiagen RNeasy Universal Mini kit (product number: 73404) and the manufacturer's protocol with some minor modifications were used to extract RNA from blubber tissue.Modifications included increased tissue usage (≥50 mg) ground in liquid nitrogen with a mortar and pestle for homogenization, the addition of a 2% non-ionic surfactant solution (Triton X-100) in the lysis buffer to help break down the high lipid content of the blubber tissue, and an additional DNA digestion using the Qiagen RNase-free DNase set (product number: 79254) according to manufacturer protocol.After extraction, RNA was collected with a single 50 μL H 2 O elution and RNA yields (>9 ng/ μL) were measured using a NanoDrop 1000 spectrophotometer (Thermo Fisher Scientific).
cDNA synthesis was performed using the Bioline Tetro cDNA synthesis kit and random hexamers (product number: BIO-65042) according to manufacturer protocols.Successful cDNA synthesis was confirmed (approximately 200 ng/μL) by DNA yield measurement using a NanoDrop 1000 spectrophotometer.

| Gene selection and primer design
To compare gene expression between individuals, PCR primers were selected for control transcripts (constitutively expressed 'housekeeping' or 'reference' genes) to account for differences in RNA quantity due to variation in the sample collection, RNA extraction and cDNA synthesis.Transcripts for frequently used controls (Ywhaz, Beta-Actin, RPL-4, COX-1 and GAPDH) were identified by aligning known target transcripts in other species with the humpback whale megNov1 genomic assembly (GCA_004329385.1)using the Basic Local Alignment Search Tool (BLAST; Altschul et al., 1990).Primers for the identified regions were designed using the GenBank primer designing tool, Primer-BLAST (Ye et al., 2012) and supplied by Sigma Oligos (Sigma Aldrich).This process was repeated for transcripts considered to be candidate health indices for humpback whales.
Gene transcripts selected as candidates for energetic condition indices included (1) LEP, coding for leptin; (2) LEPR, coding for leptin receptor; and (3) ADIPOQ, coding for adiponectin.Additionally, gene transcript (4) AHR, coding for aryl hydrocarbon receptor was selected as an index of stress related to xenobiotic metabolism.Gene transcripts (5) TNF-alpha, coding for tumour necrosis factor-alpha, and (6) IL6, coding for interleukin 6, was selected for their immunoinflammatory function.Finally, gene transcript (7) HSP-70, coding for a 70-kDa heat shock protein, was also selected as an indicator of physiological stress, particularly due to its role in counterbalancing inflammation caused by the prior two transcripts.An initial specificity check for each primer set was carried out using the UCSC in Silico PCR search tool against the taxonomically closest available genome, the northern minke whale (Balaenoptera acutorostrata, balAcu1).The final primer set (Table 1) was selected based on reliable amplification and similar melt and annealing temperatures, which is critical for each gene to be concurrently amplified in single RT-qPCR assay platform.
Using a subset of 10 samples, primers for each gene transcript (including multiple sets for Beta-Actin, RPL-4 and LEP) were tested in a reaction containing Thermo Fisher Powerup™ Syber® Green Master Mix (product number: 4368577), 10 −1 cDNA and 0.25 μM of primer (for each forward and reverse per reaction).Each RT-qPCR thermocycling profile consisted of 2 min at 50°C and 2 min at 95°C, followed by 43 cycles of 15 s at 95°C, 15 s at 64°C and 60 s at 72°C.
After initial amplification, a melt curve step to measure the relative fluorescence units with 0.5°C intervals followed the final PCR cycle.Following RT-qPCR amplification, high-resolution melt analysis was used to confirm that PCR products of multiple samples and replicates produced similar melt profiles for amplicon sequences, indicating high specificity.PCR product specificity was additionally confirmed by Sanger sequencing performed by the Australian Genome Research Facility (AGRF).Sequences were aligned with known humpback whale genomic sequences in the BLAST database to confirm target specificity.

| RT-qPCR analysis
Samples were amplified in triplicate using a 10 −1 dilution of the cDNA template per reaction.All primers had a final concentration of 0.25 μM per reaction.Thermo Fisher Powerup™ Syber® Green Master Mix (product number: 4368577), forward and reverse primers, and diluted template were added to a 96-well plate and run using a Bio-Rad CFX96 Touch thermocycler.Primers were run together on a plate such that each row contained one of each primer set while each column was loaded with a sample, allowing for a control gene and seven experimental genes to be run concurrently on a single plate.Run in triplicate, this allowed four samples per plate.Primer efficiency was determined by running a dilution series consisting of six 10x dilution steps starting with the undiluted template (in triplicate).For some primer sets, higher dilutions went beyond the range of detection, and these reactions were excluded from the analysis.
The thermocycling profile of all PCR runs consisted of 2 min at 50°C and 2 min at 95°C, followed by 43 cycles of 15 s at 95°C, 15 s at 64°C and 60 s at 72°C.A melt curve step to measure the relative fluorescence units with 0.5°C intervals followed the final PCR cycle.
High-resolution melt analysis was performed using the Bio-Rad Precision Melt Analysis™ software (product number: 1845015) to identify reactions that produced non-specific products.These reactions were excluded from subsequent analysis.

| Model controls
The Bayesian model used in this study included all the gene transcripts simultaneously.This design is intended to normalize expression at the level of the individual, effectively replacing the need for a control with the normalization of several target transcripts (Matz et al., 2013).However, to test the assumption that controls were not required in a novel application, three models were constructed with different control parameters.The 'naïve' model contained all transcripts (seven experimental and RPL-4 as a control), using a diffuse normal prior with a large variance (10 8 ) for each gene.In this model, the control was treated as another experimental gene in the normalization process, meaning it varied freely with equal weight to the other seven genes.A 'control' model was also created, which allowed for a maximum variance of 1.2-fold on average for the control gene (RPL-4).Unlike the typically used normalized relative expression (Delta Delta Ct method), where the control gene's expression is assumed to be fixed, this allowed some flexibility in control gene expression, but the experimental genes were still considered in relation to a relatively stable control gene.Finally, a 'no control' model was constructed which contained only the experimental genes.This model had identical parameters to the naïve model but did not include data for RPL-4.If the normalization technique was robust regardless of control gene presence, each model should produce similar results.However, differences between the models may suggest global expression changes relative to the control, or changes in the control regulation itself.Each model was constructed according to the formulas in the following section.

| Differential expression modelling
Statistical analysis of RT-qPCR data was carried out in RStudio version 2022.07.0 using R version 4.2.1 (R Core Team, 2022).First, raw RT-qPCR data (i.e.quantification cycle detected, Cq) were converted to transcript counts representing the number of molecules of each transcript present in each RT-qPCR reaction.To calculate transcript count, primer efficiencies were determined by fitting a linear model to the dilution series of each primer set (as described in Pfaffl, 2001).
The slope of this model was then used to determine the primer efficiency.
To estimate the Cq value at which a single transcript would be detected (Cq1) for each primer set, the formula developed for the TA B L E 1 Primer sets used for RT-qPCR analysis.Note that source accession indicates a previously identified transcript that was aligned to the megNov1 assembly.All primers were designed from the closest aligned sequences in the humpback whale genome.

| Assay validation
Among the candidate control gene transcript primers (GAPDH, YWHAZ, Beta-Actin, RPL-4 and COX-1), consistent amplification was observed for Beta-Actin, GAPDH, RPL-4 and COX-1 primer sets.The specificity of amplified products for each primer set was confirmed by aligning product sequences to known targets using BLAST.From this set of controls, RPL-4 was selected as the most stable transcript using the BestKeeper housekeeping selection tool (Pfaffl et al., 2004).Among the experimental genes, six transcripts showed consistent amplification and specific product sequences (LEP, LEPR, ADIQ, AHR, TNF-alpha, HSP-70).An additional transcript, IL-6, amplified inconsistently but was also confirmed to produce specific product sequences.SERPINA-12 amplified inconsistently and showed high variability in product sequence (evidenced using high-resolution melt analysis).It was therefore not included in the final assay design.The final assay design contained primer sets for seven experimental transcripts (excluding SERPINA-12) and one control transcript (RPL-4, Table 1).IL-6 was measured in the assay but produced a binary expression profile with a high proportion of non-detect reactions.The IL6 results were not suitable for the GLMM methodologies used, and so IL6 expression measurements were excluded from subsequent analysis.

| Model controls
To assess appropriate control techniques in adipokine analysis, three

| Differential expression modelling
The greatest changes in gene expression were observed between northbound and southbound individuals (Figure 1a).During the southbound migration, stress and inflammation-regulating transcript HSP-70 expression were significantly upregulated (p = .05),while the expression of TNF-alpha, which was included due to its immune functions of regulating inflammation and fighting infections and cancer via apoptosis, was significantly downregulated (p = .03).Likely expression increases (>87.5% of MCMC simulations showed an increase in expression) were also observed during the southbound migration for stress and xenobiotic metabolism transcript AhR and energy regulation and appetite transcript LEP expression (p = .15and p = .22respectively).
An individual's sex did not affect the expression of these experimental adipokines (Figure 1b, model results in supplementary materials).However, the lactation status of females showed some likely patterns in gene expression.For lactating females, AhR and HSP-70 were likely downregulated (p = .17and p = .22respectively) compared to non-lactating individuals.All but one of the lactating females included in this study were sampled during the southbound migration.This means that the expression of these genes for the 'average lactating female' may be distinct from the typical upregulation in these transcripts that are observed in southbound males and non-lactating females (Figure 1a vs.This study describes a RT-qPCR assay with the potential to measure health from minimally invasive outer blubber biopsy tissue. The assay measures the expression of a control gene transcript (RPL-4) in addition to seven gene transcripts (LEP, LEPR, ADIPOQ, AHR, TNF-alpha, IL-6, HSP-70) that may play roles in humpback whale energetics, stress, immune regulation and migration.The expression of these genes may also vary due to a range of naturally occurring intraspecific variation; thus, this study aimed to parse out significant differences in gene expression due to nonhealth-related individual features (i.e. and lactating/ non-lactating females) as well as changes due to migration stage (northbound/southbound), which likely act as a proxy of nutritional state.Overall, the results of this study indicate that the expression of the selected genes was minimally influenced by sex.However, likely variation in gene expression was observed due to lactation status and significant variation was observed between migration stages.These findings emphasize the importance of considering an individual's migratory stage as a crucial piece of information in health-related studies.Conservative approaches should also consider lactation status as a potential source of gene expression variation.By identifying these sources of variation, future studies can account for these effects on 'baseline' gene expression.This their southward migration (returning to the feeding grounds).
One likely explanation for this observation was that it was a response to caloric restriction, a form of dietary stress.During the migratory season, humpback whales are in an extreme caloric deficit.Previous studies have identified upregulated HSP-70 during fasting periods in humans and mice (Ehrenfried et al., 1996;Shahabi, 2011), and in metabolically demanding fasting periods for moulting pinnipeds that bear similarities to the humpback migratory period (Bennett et al., 2014;Boyd et al., 1993).Inflammatory cytokines, including TNF-alpha, are also known to be downregulated with excessive weight loss in medical studies (Moschen et al., 2010).Thus, the HSP-70 and TNF-alpha expression patterns observed for humpback whales seem to be parallel to caloric restriction effects in other aquatic and terrestrial species.A caloric restriction hypothesis may also explain the non-significant but likely trend of increasing adipokine expression during the southbound migration (see Figure 1a).This global increase may be due to an 'outward burn' with outer blubber tissue becoming more metabolically active as inner blubber reserves are depleted throughout the migration.For health studies, this would mean that the relationships between experimental genes may be of more interest than the mere fold-change values across conditions (i.e. if most transcripts are being upregulated between migratory stages, is one transcript being upregulated more than another?).
An alternate, and less likely, explanation for this pattern is the presence of acute stressors.In other marine mammal studies, increased HSP-70 expression (along with AhR which showed a non-significant upregulation trend) has been linked to an increased presence of environmental contaminants (Fossi et al., 2014;Jia et al., 2015;Lehnert et al., 2017).In this study, the proportion of cohorts sampled (juveniles, males, non-lactating females and lactating females) is known to have shifted over the course of the migratory season.These cohorts may have differences in blubber pollutant load due to individual foraging behaviour and metabolism (Remili et al., 2020), indicating that increased AhR may signal a shift towards a more affected cohort.However, little foraging occurs during the migratory season, meaning that these expression differences would have had to persist up to a couple of months after feeding activities, suggesting the metabolism of pollutants stored in blubber or an effect of increased pollutant concentrations as lipids are metabolized and pollutants are left behind.The primary reason that the acute stressor explanation seems unlikely is that other gene expression changes that should be expected in an acute response were not observed.In an inflammatory response, HSP-70 is thought to promote the production of TNF-alpha (Tsan & Gao, 2004) although there is also evidence for the converse, with upregulated HSP-70 functioning to protect cells from the apoptotic effects of TNF-alpha (Jacquier-Sarlin et al., 1994;Meng & Harken, 2002).Regardless of the direction of the gene expression pathway, if the upregulated HSP-70 was related to an acute stressor, we should also expect to see an upregulation rather than the observed downregulation of TNF-alpha.A priority for future research can be to confirm that the opposing regulation of HSP-70 and TNF-alpha observed here is tied to declining energy stores (e.g.body condition or adiposity) and not other known biomarkers of acute stress (e.g.glucocorticoids).
Unlike in migratory terrestrial animals, no significant change was observed in the expression of common appetite regulation genes, LEP, ADIPOQ and LEPR between the migratory stages in humpback whales.Adaptations to the function of these genes in long-migrating mysticetes have been previously suggested (Ball et al., 2013).This study provides some evidence that leptin and adiponectin do not serve to regulate the southbound (return to feeding grounds) migration by conveying condition information as they do in some terrestrial species (or at least these changes are not detectible in outer blubber).However, the data in this study do not provide insight into what drives the shift between the summer feeding period and the migration towards the breeding grounds.
Perhaps, a molecular approach to body condition measurement of cetaceans requires a seasonal shift in measured genes.For example, appetite hormones may signal adequate energy reserves to begin the northbound migration while other stress and cytokine pathways could signal sufficient depletion of energy stores to initiate the southbound departure from breeding to feeding grounds.
Ideally, expression-based data should be collected at both ends of the migratory corridor and compared to body morphometry to validate these theories.
In this study, 12 of the 13 lactating females were sampled during the southbound migration.Thus, we should expect to see a southbound bias with a neutral to upwards trend in HSP-70 and AhR expression compared to that of general non-lactating individuals.Instead, a likely downward trend in expression was observed for these genes.The possible attenuation of these genes fits with what may be expected in lactating individuals.During the breeding period and southbound migration, lactating females undergo the greatest shift in energy reserves as they fast and produce calorically intensive milk for the growth of their calf (Braithwaite et al., 2015;Chittleborough, 1957;Christiansen et al., 2016).At the collection site of this study, these lactating females are at the beginning of the southbound migration and in the middle of this period of intensive milk production.Thus, it is not surprising that while other reproductive classes (males and non-lactating females) may be beginning to prioritize energy conservation for the southbound migration, lactating females are still prioritizing the metabolism of blubber for milk production.This pattern is consistent with measurements of attenuated HSP-70 expression in milk somatic cells of calorically restricted dairy cows (Eitam et al., 2009), suggesting that lactating whales naturally adapted similar pathways as artificially selected cattle to prioritize milk production even during calorically restricted periods.
Future studies may benefit from the consideration of lactation status both as a natural source of variation (i.e.typical expression) and as the prioritization of milk production potentially puts lactating females at a higher risk should there be declines in prey availability or stressors that increase energy expenditure during the migratory period.
This study relied on a modelling-based approach to qPCR analysis, namely a GLMM using MCMC sampled distributions.This offers advantages in predictive power over traditional approaches, such as relative expression and 2^-ddCT.One advantage of the MCMC GLMM method is the ability to normalize experimental genes instead of using a control.In this study, the informed and no control models produced similar predictions and supported normalization being an appropriate technique in these instances.
However, the naïve model deviated from the other two models showing that a subset of stable transcripts among a global change in the other transcripts can lead the model to incorrectly fit opposing regulation (i.e. the stable transcript is assumed to be regulated in the opposite direction of the other transcripts).In this study, the stability of RPL-4 was known as it had already been identified as a stable reference gene in other cetacean studies (e.g.Martinez-Levasseur et al., 2013), and these results were independently confirmed in this study using the BestKeeper housekeeping selection tool (Pfaffl et al., 2004).This suggests that there is a caveat to the idea that normalization is an effective replacement for the use of control genes.If a small number of stable transcripts are unidentified in a model with an otherwise global increase, the stable transcripts will appear to be downregulated.This may also significantly change the results and interpretation of experimental transcripts, as is the case here with TNF-alpha being considered upregulated compared to RPL-4 but downregulated compared to the global adipokine trend.By constructing both informed and naïve models, global expression trends can be identified from the naïve model and more detailed patterns in gene expression within this trend can also be examined using the informed model.Ideally, studies using this methodology should include pre-validated controls and continue using a multimodel (informed, naïve and no control) approach to ensure that stable experimental genes are not accidentally misclassified as the known control RPL-4 was in the naïve version.
In summary, this study presents a health assay capable of concurrently measuring a range of health-related gene transcripts.
Changes in gene expression due to potential natural sources of variation were quantified (migratory stage, sex and female lactation status), so these sources of variation can be effectively considered when using this assay for health assessments.Results showed that migratory stage has a significant impact on HSP-70 and TNF-alpha that is consistent with fasting effects in other species.However, common energetic hormones (LEP, ADIQ, LEPR) showed the opposite effect from terrestrial species (upregulated instead of downregulated during fasting), likely due to increased metabolic activity in outer blubber as an individual's energy stores are depleted.Unlike migratory stage, sex did not change the expression of any target genes, meaning that sex is unlikely to be a confounding variable in health studies using this method.Likely differences in gene expression were also observed for lactating and non-lactating females, suggesting that lactation should also be considered as a potential source of variation.Overall, the results indicate that this assay has the potential to quantify pathways relating to energetics and stress.Thus, the expression of the genes measured in this study may be a powerful diagnostic tool to quantify individual health, and with adequate sample sizes, to scale this health information to identify population-level trends and disturbances.
76104) and stored on ice while on board the vessel.Samples were immediately transferred to −80°C onshore the same day (1-8 h from sampling).Storage continued at −80°C until RNA extraction, which occurred between 3 and 10 months after sample collection.
instrument in Matz et al. (2013) was generalized to the Bio-Rad CFX96 Touch thermocycler used in this study.The use of this equation was validated when low transcript variance was observed in one of the primer set dilution series reactions (HSP-70), allowing us to confirm that the observed Cq1 approximately matched the value calculated using this formula.Because the Cq1 value has little influence over relative quantification analysis (Matzet al., 2013), it was assumed that the remaining Cq1 estimates using this formula were well within allowable error (Cq1 ± 3).With this information, the number of target molecules present in each reaction before qPCR (i.e. the transcript count) was determined using the cq2counts function in the MCMC.qpcrpackage(Matz et al., 2013).Next, a generalized linear mixed model (GLMM) with Poissonlognormal errors was fitted using a Markov Chain Monte Carlo method (MCMC) using the MCMC.qpcrfunction.This Bayesian approach was selected due to its ability to handle variation in the control transcript, bias and applicability to qPCR data that includes reactions with high Cq values and non-amplification.The count of each transcript was modelled using fixed terms representing the target transcript (RPL4, LEP, LEPR, ADIPOQ, AHR, TNF-alpha,, and fixed interactions between the target transcript and sex, migratory stage and lactation status and a interaction between target transcript, sex and migratory stage were used to account for sex-specific differences in gene expression within each migratory stage.Random effect terms accounted for variance in cDNA preparations and variance of sample replicates for each gene from each individual.To avoid high variability in early MCMC simulations, a burn-in (the number of samples discarded from the beginning of the MCMC chain) of 4500 samples was selected by fitting several model iterations and observing the approximate simulated sample size where the model reached stability.The minimum effective sample size for each parameter was calculated to be 2124 for a 95% confidence interval with 90% tolerance (i.e. the minimum effective samples for the mean of the posterior distribution to be within 10% of one standard deviation from the true mean 95% of the time) using the methods outlined byVats et al. (2019).To consistently meet the minimum effective samples required, each MCMC was constructed from 50,000 iterations.Fitted MCMC GLMMs were assessed for violations of the assumptions of residual linearity, homoscedasticity and normality.As Bayesian models used simulated data sets, the results reflect the likelihood of repeating the finding given a new set of observations (as opposed to the rejection of a null hypothesis).Significance was defined as p < .05using two-tailed p-values calculated in the MCMC.qpcrpackage such that an observed change in expression (either upregulation or downregulation) that occurs in 97.5% of MCMC iterations produces a p-value of 0.05.Observed changes in expression that did not reach significance but had >87.5% likelihood of occurrence (p < .25)were discussed as non-significant but 'likely' so that future studies can opt to take a conservative approach in considering potential sources of baseline variation.A full description of the equations used to calculate model parameters and p-values can be found in the supplementary materials.
Chain Monte Carlo generalized linear mixed models were constructed.The first represented the 'naïve' assumption that all transcripts were equally weighted.The second represented an 'informed' assumption, where the control gene RPL-4 was assumed to have a maximum variance of 1.2-fold (i.e. the control is mostly stable) among all individuals.The third was a 'no control' model, where experimental genes were equally weighted, and RPL-4 was absent from the analysis.Results show agreement between the three models in almost all the expression change estimates of each transcript (LEP, LEPR, ADIQ, AHR, TNF-alpha, HSP-70, RPL-4) for each effect (migratory stage, sex and female lactation status).However, there were differences in the modelled expression of TNF-alpha and RPL-4 between the northbound and southbound migrations (Figure1a, see * in panels).For the control gene RPL-4, a slight upregulation for southbound individuals was found in the informed model, and a substantial (about −1fold change) downregulation was observed in the naïve model.For TNF-alpha, southbound individuals showed a substantial downregulation in the informed and no control models, with an upregulation in the naïve model.If the global upregulation from the naïve model were observed for all transcripts (including RPL-4) that may indicate technical issues, such as differences in RNA integrity between migratory stages.However, the naïve model's southbound decline in RPL-4 expression indicates that the observed global increase (see ADIQ, AhR, HSP-70, LEP; Figure 1a) is an effect of the migratory stage.In the informed model, all the southbound expression increases were attenuated compared to the naïve model.However, the extreme decline in TNF-α represents a true downregulation compared to the global adipokine expression trend.This is confirmed by the results of the no-control model, where only adipokines are present (a).Thus, the naïve model informs us that there is a global upregulation of the experimental adipokine transcripts, while the informed model represents the most precise information on experimental transcript expression, accounting for this global change.For the remainder of these results, we refer to the informed model as the detailed proportional trends in adipokine expression are of most interest in health applications.
Figure 1c, model results in supplementary materials.This distinction also explains the likely changes in expression due to sex/migratory stage interactions (AhR, p = .21and HSP-70, p = .24).F I G U R E 1 Expression changes for the informed, naïve and no control models between (a) southbound individuals compared to northbound, (b) males compared to females and (c) pregnant versus non-pregnant females.A '0' value represents the average expression of the latter class, while dots represent the average fold difference of the former class.The bars represent the 95% credible interval of the fold change.Migration direction yielded the only significantly different shifts in the expression of HSP and TNF.* indicates disagreement between model results as discussed in 'Analysis of Differential Adipokine Expression' section.
will allow gene expression to be compared to other known indices (e.g.steroid hormones and body condition) to accurately measure individual health.The most pronounced source of natural variation identified in this study was migratory stage.Within this variation, the most dramatic gene expression changes were an upregulation of HSP-70 and an opposing downregulation of TNF-alpha, as individuals shifted from their northward (towards the breeding grounds) to