The multi-tissue gene expression and physiological responses of water deprived Peromyscus eremicus

Abstract

The harsh and dry conditions of desert environments have resulted in genomic adaptations, allowing for desert organisms to withstand prolonged drought, extreme temperatures, and limited food resources. Here, we present a comprehensive exploration of gene expression across five tissues (kidney, liver, lung, gastrointestinal tract, and hypothalamus) and 19 phenotypic measurements to explore the whole-organism physiological and genomic response to water deprivation in the desert-adapted cactus mouse (Peromyscus eremicus). The findings encompass the identification of differentially expressed genes and correlative analysis between phenotypes and gene expression patterns across multiple tissues. Specifically, we found robust activation of the vasopressin renin-angiotensin-aldosterone system (RAAS) pathways, whose primary function is to manage water and solute balance. Animals reduced food intake during water deprivation, and upregulation of PCK1 highlights the adaptive response to reduced oral intake via its actions aimed at maintained serum glucose levels. Even with such responses to maintain water balance, hemoconcentration still occurred, prompting a protective downregulation of genes responsible for the production of clotting factors while simultaneously enhancing angiogenesis which is thought to maintain tissue perfusion. In this study, we elucidate the complex mechanisms involved in water balance in the desert-adapted cactus mouse, P. eremicus. By prioritizing a comprehensive analysis of whole-organism physiology and multi-tissue gene expression in a simulated desert environment, we describe the complex response of regulatory processes.

Introduction

Genomic adaptations play a pivotal role in enabling life to persist in the harsh and dynamic conditions of desert environments. Evolutionary processes have shaped the genomes of these organisms to enhance their capacity to withstand prolonged drought, extreme temperatures, and limited food resources [1,2,3,4,5]. Understanding these genetic underpinnings of desert adaptation not only contributes to our comprehension of evolutionary biology, but it also holds promise for insights into how these adaptations may be leveraged to address challenges posed by water scarcity and climate change in other organisms, including humans, and in other ecosystems. Studies of desert mammals have provided evidence of positive selection on genes related to food storage [4, 6], water reabsorption [5,6,7,8], osmoregulation [1, 9,10,11], fat metabolism [1, 3, 12,13,14], thyroid-induced metabolism [15], and salt regulation [16]. These genetic insights suggest the molecular basis of observed phenotypes, including enhanced metabolic water production [17,18,19], reduced water loss [17, 20, 21], tolerance to high-salt diets [6, 22], and coping with starvation and dehydration [9, 20, 23,24,25], which are all common in desert-dwelling mammals. However, it remains unclear how water deprivation affects gene expression in individual organs, and how this in turn contributes to whole-organism phenotypes.

Osmoregulation, or the process through which animals manage water and solute balance, is critical for desert animals. It involves the maintenance of internal fluid homeostasis, with water intake being dependent on factors such as drinking, dietary sources, and metabolic water, while water output is regulated through processes including waste removal (i.e., urine and feces), respiration, perspiration, and reduced food intake [26,27,28]. Failure to maintain water and solute homeostasis can result in impaired renal, reproductive, and cardiovascular function, affect an animal’s ability to regulate its body temperature, and ultimately lead to death [27]. During water deprivation, a decrease in extracellular water volume results in heightened plasma osmolality due to an elevated concentration of solutes, primarily sodium, which is detected by osmoreceptors [29,30,31]. Osmotic balance is then intricately managed through two distinct multisystem mechanisms, the renin-angiotensin-aldosterone system (RAAS) and by vasopressin [26, 29, 32, 33].

In response to changes in osmotic pressure, osmoreceptors in the hypothalamus become activated, aldosterone is released from the adrenal glands, and vasopressin is produced in the hypothalamus and released from the posterior pituitary gland [30, 31, 34]. This promotes water reabsorption in the kidneys by enhancing the water permeability of epithelial cells lining renal collecting ducts [35]. Additionally, the kidneys retain sodium while also excreting solutes via urination [36]. Proteins related to the transport of water are translocated from within the cell to the cell surface, forming water channels, resulting in increased water reabsorption from the tubule system of the nephrons back into the bloodstream [37,38,39] and the retention of sodium in the distal tubules [40]. Renin is released from the kidney, and acts on angiotensinogen that is released from the liver [29, 41]. This triggers the formation of angiotensin I, which is then converted into angiotensin II, a hormone involved in regulating blood pressure and peripheral circulation, by the release of angiotensin-converting enzyme in the lungs [29, 42, 43]. As a result, RAAS impacts water and solute retention by promoting vasoconstriction, stimulating the release of aldosterone, and enhancing the reabsorption of sodium and water in the kidneys. This orchestrated response helps regulate blood pressure and maintain fluid balance in the body.

If water homeostasis is not achieved, blood volume continues to decrease, resulting in hemoconcentration. As dehydration leads to a higher concentration of erythrocytes, blood viscosity increases, enhancing the risk of spontaneous clot formation and hindering blood flow through vessels. This elevated blood viscosity might impede the efficient delivery of oxygen and nutrients to various tissues and organs, triggering alterations in vascular dynamics in nutrient-deprived, hypoxic environments [44,45,46]. Studies suggest that impaired flow-dynamics typical of dehydration may induce reversible changes in angiogenesis, altering the growth and development of blood vessels to regulate blood flow and distribution [47, 48]. The modulation of angiogenesis during dehydration reflects the body’s dynamic response in adjusting vascular networks to manage the metabolic demands of tissues.

Vasopressin receptors expressed in the kidneys, lungs, liver, hypothalamus etc., further affect vasoconstriction, glycogenolysis, water reabsorption, thermoregulation, and food intake [34]. Water deprived animals have been shown to limit food intake [49] as an adaptive mechanism allowing for osmotically sequestered water in the gastrointestinal (GI) tract to be reabsorbed into the systemic vasculature [50,51,52] which thereby reduces solute load [53]. Decreased food intake may be driven by the expression of vasopressin in the hypothalamus [34], suggesting another link between eating and drinking. To survive reduced food intake during water deprivation (i.e., dehydration anorexia) while maintaining blood glucose concentrations, previous rodent studies have shown increased glycogenolysis, lipolysis, and/ or gluconeogenesis [28, 54, 55].

In this this study, we have performed a comprehensive analysis of gene expression across five tissues relevant to the response to dehydration (kidney, liver, lung, gastrointestinal tract, and hypothalamus) and used 19 phenotypic measurements to assess the whole-organism physiological and genomic response to water deprivation in a hot and dry environment in the desert-adapted cactus mouse (Peromyscus eremicus). The results of this study include: (1) a robust activation of RAAS, as seen by upregulation of AGT across all five tissues, (2) upregulation of PCK1, reflecting an adaptive response to maintain blood glucose levels during decreased oral intake, (3) a broad decrease in genes related to coagulation, possibly in response to hemoconcentration and (4) a clear signal of vascular remodeling. Overall, the lung experienced the largest number of changes in gene expression, followed by tissues involved in RAAS and then the hypothalamus.

Methods

Animal Care, RNA extraction, and sequencing

Captive born, sexually mature, non-reproductive healthy male and female P. eremicus were reared in an environmental chamber designed to simulate the Sonoran desert [20, 24, 56, 57]. All mice were subjected to standard animal care procedures before the experiment which included a health assessment conducted by licensed veterinary staff following animal care procedures guidelines established by the American Society of Mammologists [58] and approved by the University of New Hampshire Institutional Animal Care and Use Committee under protocol number 210,602. Mice were provided a standard diet and fed ad libitum (LabDiet^® 5015*, 26.101% fat, 19.752% protein, 54.148% carbohydrates, energy 15.02 kJ/g, food quotient [FQ] 0.89). Animals were randomly selected and assigned to the two water treatment groups (n = 9 of each group, female mice with water, female mice without water, male mice with water, and male mice without water, total n = 36). Prior to the start of the experiment, a temperature-sensing passive integrated transponder (PIT) tag (BioThermo13, accuracy ± 0.5 °C, BioMark^®, Boise, ID, USA) was implanted subdermally. At the start of the experiment (day 0, time 0 h, 10:00), mice were weighed (rounded to the nearest tenth of a gram) on a digital scale and water was removed from chambers corresponding to those animals in the dehydration group. Mice were metabolically phenotyped for the duration of the experiment [20] using a pull flow-through respirometry system from Sable Systems International (SSI). Rates of CO₂ production, O₂ consumption, and water loss were calculated using Eq. 10.6, 10.5, and 10.9, respectively, from Lighton [59]. Respiratory quotient (RQ, the ratio of VCO₂ to VO₂) and energy expenditure (EE) kJ hr^− 1 were calculated as in Lighton [59], [Eq. 9.15]. For downstream analysis, we calculated the mean of the last hour of water loss, EE, and RQ for each mouse.

At the conclusion of the experiment (day 3, time 72 h, 12:00) as described in Blumstein and MacManes [20], body temperature was recorded via a Biomark^® HPR Plus reader, mice were weighed, animals were euthanized with an overdose of isoflurane, and 120 µl of trunk blood was collected for serum electrolyte measurement and analyzed with an Abaxis i-STAT^® Alinity machine using i-STAT CHEM8 + cartridges (Abbott Park, IL, USA, Abbott Point of Care Inc). We measured the concentration of sodium (Na, mmol/L), potassium (K, mmol/L), blood urea nitrogen (BUN, mmol/L), hematocrit (Hct, % PCV), ionized calcium (iCa, mmol/L), glucose (Glu, mmol/L), osmolality (mmol/L), hemoglobin (Hb, g/dl), chlorine (Cl, mEq/L), total CO₂ (TCO₂, mmol/L), and Anion gap (AnGap, mEq/L). Using Na, Glu, and BUN, we calculated serum osmolality. To test for statistically significant (p < 0.05) differences, we used a student’s two-tailed t-test (stats::t.test) between the sexes for each experimental group in R v 4.0.3 [60].

The lung, liver, kidney, a section of the large intestine (referred to as GI throughout), and hypothalamus were collected and stored in RNAlater (Ambion) at 4 °C for 12 h before being frozen at − 80 °C for long-term storage. Prior to RNA extraction, the tissues were removed from the RNAlater and a small section was dissected off. Care was taken to retain an anatomically similar region of tissue from each animal. Tissues were mechanically lysed using a Bead Beater, and RNA was then extracted using a standardized Trizol protocol. RNA libraries were prepared using standard poly-A tail purification, prepared using Illumina primers, and individually dual-barcoded using a New England Biolabs Ultra II Directional kit (NEB #E7765). Individually barcoded samples were pooled and sequenced paired end and 150 bp in length on two lanes of a Novaseq at the University of New Hampshire Hubbard Center of Genome Studies.

Genome alignment and Differential Gene expression

All the code used to analyze the data is located at the GitHub repository (https://github.com/DaniBlumstein/dehy_rnaseq). The P. eremicus genome version 2.0.1 from the DNA Zoo Consortium (dnazoo.org) was indexed, and reads from each individual were aligned to the genome using STAR version 2.7.10b [61] allowing a 10 base mismatches, a maximum of 20 multiple alignments per read, and discarding reads that mapped at <30% of the read length. Aligned reads were counted using HTSEQ-COUNT version 2.0.2 [62] with the flags ‘-s no -t exon’ which directs the software to use only uniquely mapped reads. All other flags for HTSEQ-COUNT were set to the default parameters.

Counts from HTSEQ-COUNT were exported as csv files, and all downstream statistical analyses were conducted in R v4.0.3 [60]. Tissue samples were removed based on low mapping rates and PCA outliers as these two criteria suggest poor sample quality, sample contamination, or mis labelling leaving 163 samples (Females without water: 9 GI, 9 hypothalamus, 8 kidney, 8 liver, 8 lung, 42 total. Males without water: 7 GI, 8 hypothalamus, 7 kidney, 8 liver, 8 lung, total 38. Females with water: 7 GI, 8 hypothalamus, 8 kidney, 8 liver, 8 lung, 39 total. Males with water: 9 GI, 9 hypothalamus, 9 kidney, 9 liver, 8 lung, 44 total). Counts were aggregated into a gene-level count by combining all counts that mapped to the same gene. Low expression genes (defined as having 10 or less counts in 8 or more individuals) were removed from downstream analyses. Differential gene expression analysis was conducted in R using DESEQ2 [63]. For the dataset as a whole, we performed three models to test for the effects of sex, water access, and tissue type. For each tissue, we performed two models, testing the effect of and identifying genes specific to sex and water access with a Wald test. Results were visualized using GGPLOT2 [64]. We then conducted all downstream analyses (except for CCA, see below) on each tissue independently. Count and sample data were filtered to include only the tissue of interest and low expression genes were removed.

Weighted Gene Correlation Network Analysis

To identify the regulation of gene expression associated with responses to water access, we performed a weighted gene correlation network analysis (WGCNA), a network-based statistical approach that identifies clusters of genes with highly correlated expression profiles (modules), [65] for each tissue independently. This approach allows us to relate gene expression with physiological phenotypes (mean EE, water loss, RQ, total weight loss, proportional weight loss, sex, body temperature, water access, and the panel of electrolytes). Prior to WGCNA, read counts were normalized within tissues using DESEQ2 [63]. Module detection was done using WGCNA::blockwiseModules with networkType set to “signed” but otherwise default parameters were used. We estimated a soft threshold power (β) for each tissue dataset by plotting this value against mean connectivity to determine the minimum value at which mean connectivity asymptotes, which represents scale-free topology (liver = 15, kidney = 21, GI = 14, lung = 20, hypothalamus = 14).

Canonical Correlation Analysis

We used a Canonical Correlation Analysis (CCA) implemented in the R package vegan [66] to investigate multivariate correlation of gene expression, by tissue, water access, and sex, with metabolic variables (mean EE, mean RQ, mean water loss, body temperature, and proportional weight loss) and display the three levels of information in a triplot. We used an ANOVA to identify what response variables were significant. Significant response variables were graphed as vectors and allowed us to identify their correlative nature; vectors pointing in the same direction are positively correlated, while vectors pointing in opposite direction are negatively correlated. To identify genes of interest, we selected genes that graphed two standard deviations away from the mean for CCA1 and CCA2.

Gene Ontology

To examine gene ontology of DE genes and WGCNA modules, we cross-referenced our gene IDs with Homo sapiens gene IDs via Ensembl before running Gene Ontology (GO) analyses. Each analysis above resulted in a list or lists of genes that were used as input for the GO analysis using the R package gprofiler [67]. From there, we identified the topmost 20 significant GO terms based on g: SCS corrected p-values [68] for each up and downregulated list for each tissue and for each significant module from the individual WGCNA analysis. We used REVIGO [69] to reduce GO term redundancy.

Consensus gene list and KEGG pathway analysis

We generated a high confidence consensus list of genes from the results of the three orthogonal analyses; DE, WGCNA, and genes located two standard deviations away from the origin in the CCA. We then selected three KEGG pathways (Renin-Angiotensin system KEGG pathway [hsa04614], vasopressin-regulated water reabsorption pathway [hsa04962], and insulin resistance KEGG pathway [hsa04931], [69, 70] based on genes in our consensus gene set and cross referenced the genes in those pathways with significantly differently expressed genes in our five tissue datasets.

Results

No health issues were detected by veterinary staff, no animals were removed prior to the end of the experiment, and all mice were active at the end of the experiment.

Genomic data

We obtained an average of 21.44 million reads (+- 11.6 million SD) per sample (PRJNA1048512). On average, 78.33% of reads were uniquely mapped per sample (+- 2.12% SD). Data on the number of reads and mapping rate per sample are located in Supplemental File 1, raw read files are archived at NCBI SRA BioProject: PRJNA1048512, and all gene expression count data and code used to analyze the data are located at the GitHub repository (https://github.com/DaniBlumstein/dehy_rnaseq).

Electrolytes and physiological phenotypes

The same mice used to generate the electrolyte and physiology data more fully described in Blumstein and MacManes [20] are also used in the study described herein for RNAseq analysis. When comparing males and females separately, the following electrolytes were significantly different and elevated in water deprived mice: Na (male and female Na p = 0.0016 and p = 0.0026 respectively), BUN (p = 0.001/0.003), Hct (p = 0.002/0.001), osmolality (p = 8.2e-05/0.0001), Cl (p = 0.02/0.007), Hb (p = 0.017/0.009), and TCO2 (female p = 0.017) (Table 1). When comparing males to females within each water treatment (with or without access to water), no significant differences were found in the electrolyte levels (Table 1). Both males and females experienced significant weight loss (p = 0.001, 0.005) and proportional weight loss (p = 2.2e-16, 2.659e-09) at the end of the experiment [20]. Body temperature was significantly lower for female mice without access to water, but not for males (p = 0.0003).

To relate whole-organism physiology data to gene expression data, we calculated the means for each mouse from the last hour of data collected in Blumstein and MacManes [20], [data located at: https://github.com/DaniBlumstein/dehy_phys] for WLR, EE, and RQ for the same adult females and adult males used in this study. Within sex, WLR significantly differed between water groups (male and female, p = 0.001 and 0.002), RQ was significantly different between water groups for males (p = 0.0003), and EE was not significantly different for either males or females.

Table 1 Mean measurements for serum electrolyte measurements (hct = hematocrit (% PCV), hb = hemoglobin (mmol/L), na = sodium (mmol/L), K = potassium (mmol/L), cl = chloride (mmol/L), TCO₂ = total Carbon Dioxide (mmol/L), BUN = blood Urea Nitrogen (mmol/L), cr = creatinine (µmol/L), Glu = glucose (mmol/L), iCa = ionized calcium (mmol/L), AnGap = anion gap (mmol/L), and osmolality (mmol/L), change in weight (g), and body temperature (°C) for female (n = 18), male (n = 18), and all Peromyscus eremicus (n = 36) with and without access to water. Data were collected and are further described in Blumstein and MacManes [20]

Differential gene expression

We cross-referenced our gene IDs with Homo sapiens gene IDs. Patterns of gene expression data are largely driven by tissue type (PC1: 42% variance and PC2: 26% variance, Supplemental File 2).

For all tissues except the hypothalamus, we found many differentially expressed genes (p < 0.05) between water treatments (Fig. 1, Supplemental File 3) and few differentially expressed genes between sex (Supplemental File 3). This resulted in 12,083 genes in the kidney (1056 downregulated and 793 upregulated), 12,622 genes in the GI (703 downregulated and 956 upregulated), 13,389 genes in the hypothalamus (8 downregulated and 1 upregulated), 11,125 genes in the liver (1122 downregulated and 933 upregulated), and 12,942 genes in the lung (2183 downregulated and 2190 upregulated).

Fig. 1

Log2FC of all genes across the lung (lu), liver (liv), gastrointestinal tract (gi), hypothalamus (hyp), and kidney (kid) of Peromyscus eremicus with water vs. without water. Blue and green colored dots indicate p < 0.05, whereas grey dots indicate p > = 0.05. The number of differentially expressed genes is annotated above the point cloud for that tissue

Weighted gene correlation network analysis

A total of 12,083 genes in the kidney were successfully assigned into 12 modules (Table 2) with the number of genes per module ranging from 31 to 7008. A full list of gene assignments is available in Supplemental File 4. We identified 24 individual modules using 12,622 genes for the GI (Table 2). The modules contained 33-3530 genes each (Supplemental File 4). In the lung, 12,942 genes were assigned to 13 modules (Table 2). Each module contained 23-7116 genes (Supplemental File 4). A total of 13,389 genes were assigned to 18 different modules in the hypothalamus (Table 2). Modules contained 30-2319 genes (Supplemental File 4). Finally, 11,125 genes were assigned to 18 modules in the liver (Table 2), with the number of genes per module ranging from 28 to 3116 (Supplemental File 4).

Table 2 The number of WGCNA modules for each phenotypic measurement for the lung, liver, gastrointestinal tract, hypothalamus, and kidney of Peromyscus eremicus

Gene ontology

149 total GO terms were identified in the upregulated and down regulated differentially expressed genes across all tissue with 14 terms identified in two different tissues (Supplemental File 5). Specifically, hydrolase activity and renal system development were uniquely identified in the GI. Lipid metabolic process and detoxification were uniquely identified in the liver. Immune responses were uniquely identified in the lung. Lastly, starvation response was uniquely identified in the hypothalamus (Supplemental File 5).

For each tissue WGCNA analysis, we identified GO terms for each significant module phenotype combination. After filtering each combination for the top 20 GO terms, we identified a total of 379 unique GO terms. 199 GO terms were identified in the hypothalamus with 168 terms only appearing in one module, 25 identified in two modules, five in three modules, and one in four modules (Fig. 2). Within the kidney 82 unique GO terms were identified. Specifically, 78 terms were in one module and four terms were in two modules. There were 122 unique GO terms identified in the liver with 100 of the GO terms appearing in one module, nine terms in two modules, and three terms in three different modules (Fig. 2). Lastly, the lung had 52 unique GO terms identified, the fewest number any tissue, and there was no GO term overlap for any of the modules (Fig. 2). When comparing GO terms between tissues, 53 terms were in two different tissues, five terms were in three of the tissues, and one term were in four tissues (Fig. 2).

We found several overlapping GO terms related to vascular development across all five tissues (vascular development [GO:0001944], circulatory system development [GO:0072359], angiogenesis [GO:0001525], blood vessel development [GO:0001568], blood vessel morphogenesis [GO:0048514], and regulation of vasculature development [GO:1,901,342]). These GO terms were identified in a myriad of gene modules (water treatment, EE, Hb, TCO2, WLR, and body temperature, Fig. 2). Additionally, we found several downregulated GO terms related to hemoconcentration (regulation of coagulation [GO:0050818], blood microparticle [GO:0072562], hemopoiesis [GO:0030097], serine hydrolase activity [GO:0017171], and regulation of hemostasis [GO:1,900,046]) were identified in the hypothalamus and liver (Fig. 2) as well as. Finally, GO terms related to metabolic processes were identified in multiple significant WGCNA modules in the kidney, but not other tissues.

Fig. 2

Visualization of gene ontology (GO) terms to show common WGCNA modules within and between the lung (lu), liver (liv), gastrointestinal tract (gi), hypothalamus (hyp), and kidney (kid) of Peromyscus eremicus. Visualized are selections of the top 20 significant GO terms for each phenotype module combination based on adjusted p-value. The number of genes in the GO term are indicated by size of the dots

Canonical correlation analysis

We examined the relationship between gene expression and water access, physiological variables (EE, RQ, WLR, proportional weight loss, body temperature), and tissue type using CCA [66]. CCA suggests a significant overall association between the physiological variables and gene expression across tissues (F = 105.45, p = 0.001; CCA1–38.02% and CCA2–21.14%, Fig. 3; Table 3). We identified 1233 genes two standard deviations from the origin and found overlap between the genes located two standard deviations in the CCA and genes assigned to WGCNA modules (kidney: 37, GI: 60, lung: 69, hypothalamus: 53, liver: 40). Here, the proportional weight loss was significantly correlated (p = 0.001, Table 3) with gene expression in the hypothalamus (Fig. 3). This suggests that proportional weight loss exerts an effect on gene expression in the hypothalamus.

Table 3 ANOVA results from the Canonical Correspondence Analysis (CCA) indicating correlations between normalized gene expression and physiological measurements for Peromyscus eremicus with and without access to water. Formula: gene expression ~ water access + proportional weight loss + RQ + EE + WLR + body temperature + sex + tissue)

Fig. 3

Canonical correspondence analysis (CCA) indicates correlations between normalized gene expression and physiological measurements for Peromyscus eremicus with and without access to water. The distribution of tissue samples in CCA space as a function of their gene expression values is shown (points colored by tissue type and water treatment). Outlier genes (defined as two standard deviations or more from the mean) are colored blue. Inlier genes (defined as less than two standard deviations from the mean) are colored grey. CCA reveals a significant relationship for proportion weight loss (F = 4.5153, P = 0.001) and, while not significant, a relationship for water loss rate (WLR) (F = 1.9564, P = 0.093)

Consensus gene set

We identified 41 genes (Table 4) assigned to at least one significant WGCNA module for all five tissues as well as located at least two standard deviations away from the origin in our CCA triplot (Fig. 3). All genes identified from these two orthogonal analyses were also DE in at least on tissue. Using manual annotation, the consensus gene lists were grouped into housekeeping, ion homeostasis, central nervous system, apoptosis, coagulation, blood pressure, angiogenesis, lipid, immune response, gluconeogenesis, and bitter taste (Table 4).

Table 4 Genes identified in all three analyses (i.e. significantly differentially expressed genes, assigned to a significant module in WGCNA, and are outliers in CCA) in all five tissues in the study (lung (lu), liver (liv), gastrointestinal tract (gi), hypothalamus (hyp), and kidney (kid) grouped by the general function of the gene

Discussion

Extensive research has been conducted on the genomic and physiological mechanisms that control water balance in mice [20, 71, 72]. This is particularly intriguing in the context of exploring adaptations to extreme conditions, given that the management of water has significant implications for survival. Studies have largely focused on characterizing the response in kidneys [25, 72,73,74], which, while important, represent only a fraction of the physiological and genomic response to dehydration. Indeed, the organismal response to an environmental stressor (such as dehydration), is likely to involve coordination of multiple organ systems and physiological process. However, there are key challenges to studying such responses at the organismal level. Although we understand that organ-systems operate in tandem with other organ systems, the identification of coordination at the level of gene expression is challenging, even when organismal physiology is well-characterized. In addition to complications related to biology, technical complications exist. Here, when looking at gene expression levels across tissues, results highlight differences between the tissues rather than to elucidate the ways in which processes in one depend on the actions of the other, thereby potentially obscuring the gene expression signal of coordination. Further, this work relies on one rodent species, P. eremicus. Future studies should leverage a comparative species approach, contrasting gene expression changes during water deprivation in desert adapted and non-desert adapted rodent species to identify gene expression differences associated with adaptation to desert environments.

The study of coordination at the level of physiology has been similarly difficult to elucidate, due to limitations in our ability to collect and analyze phenotypic data at a temporal scale that is relevant to the biological phenomenon under study but see [20, 56, 57, 74,75,76]. Despite the challenges, understanding the complex interplay of multi-tissue networked gene expression with whole-organism physiological response is critical to developing a more in-depth understanding of how organisms respond to environmental stressors. Our focus on physiological measurements in a simulated desert environment, as well as on collection of multi-tissue gene expression data in the context of a whole-organism response, represents a distinctive contribution to the field. The findings underscore the complexity of the genetic landscape governing physiological responses to water deprivation and emphasize the need for a broader organismal understanding of the physiological and genomic mechanisms orchestrating successful adaptations.

Global response to acute water deprivation

The inclusion of physiology and multi-tissue gene expression data in this study provides opportunities to understand how the tissues work together to compensate for water deprivation. During water deprivation, blood volume decreases, which if inadequately managed, may result in a decrease in organ perfusion and ultimately organ failure. The multisystem response to water deprivation involves a coordinated response, including renin-angiotensin-aldosterone system (RAAS) activation, reduced food intake, with compensatory activation of systems for preserving levels of serum glucose, a widespread reduction of genes responsible for clotting factors, and distinct indications of vascular restructuring to support perfusion.

The first response to water deprivation involves the upregulation of processes that are aimed at solute management. In part regulated by the RAAS. In support of the activation of the RAAS mechanism, multiple genes in the Renin-Angiotensin system KEGG pathway (hsa04614: [70, 77] were found to be differentially expressed. Among these genes were AGT (angiotensinogen, DE in the lung, kidney, liver, and GI), ACE2 (angiotensin-converting enzyme, DE in the kidney), REN1 (Renin, DE in the kidney), CMA1 (converts angiotensin I to angiotensin II, DE in the lung), and AGTR1 (Angiotensin II Receptor, DE in the kidney), which were all more highly expressed in water-deprived mice. Further, we found two complete pathways of genes to be significantly differentially expressed genes within hsa04614, enhanced vasoconstriction and coagulation cascade as well as vasoconstriction, inflammation, fibrosis, antinatriuresis, reactive oxygen species activation, and Na and water retention [77]. Interestingly, the candidate gene and essential component of the RAAS mechanism, AGT, was significantly upregulated in dehydrated mice, assigned to a significant WGCNA module, and defined as an outlier in the CCA analysis (Table 4). Together with the physiological data, the consistent genetic signature of RAAS activation, using orthogonal analytical methods, suggests a robust whole-body RAAS response which is critical for the cactus mouse to regulate water balance.

At the same time, independent of the RAAS pathway but stimulated by its products, vasopressin is released for the primary function of conserving water. Vasopressin binds to receptors, activating a signaling cascade in the kidneys which functions to retain water by inducing expression of water transport proteins in the late distal tubule and collecting duct, increasing the permeability of the membrane to water. The increased water permeability allows water to move through, from the collecting ducts back into the bloodstream, further aiding in regulating blood volume and fluid balance [42, 43]. Support for the genetic activation of the vasopressin pathway is shown with various genes within the vasopressin-regulated water reabsorption pathway [hsa04962], [69, 70] are significantly differentially expressed. Specifically, STX4 (DE in the kidney), RAB5C (DE in the kidney, liver, and lung), RAB11B (DE in the lung), AQP3 (DE in the lung), VAMP2 (DE in the lung), DYNLL2 (DE in the lung), GNAS (DE in the lung), and AVPI1 (DE in the lung and kidney) were upregulated in dehydrated mice. Interestingly, AQP2, a key gene in hsa04962, did not have any reads mapped to it.

It is worth noting here that the identification of differentially expressed transcripts on a KEGG pathway can be thought of in terms of a stoichiometric problem, albeit one where neither current tools nor current data allow us to satisfactorily solve. When attempting to use differential expression to support the upregulation of a given pathway, it may be that only one gene, the quantity of which is rate limiting, may be more highly expressed. Without understanding which transcripts are rate limiting, the interpretation of KEGG pathway mapping may suffer from either over- or under-valuation. Lastly, gene interaction maps (e.g., KEGG) that span organ systems do not currently exist, mostly because the models used for their development (e.g., yeast, flies) lack such complexity present in mammalian systems. Future studies could instead leverage analytical approaches such as GWAS or QTL to test for functional association. Additionally, CRISPR or spatial transcriptomics could also be used to test for functional association.

During water deprivation, animals are known to reduce the amount of solid food intake [28, 49, 54, 55, 78]. Known as dehydration associated anorexia, this secondary response reduces the amount of water required for digestion and facilitates water reabsorption from the kidneys and gastrointestinal tract back to systemic circulation [28, 53]. While the magnitude of the impact in cactus mice is unknown, diets high in fiber, like the diet consumed in this study, can result in high fecal volume, and has been shown to account for as much as 25% of total daily water loss in rats [53, 79]. Perhaps even more importantly, but unquantified, the processing of solid food requires the production of significant quantities of digestive enzymes, all of which are water rich. In humans, as much as 10 L of fluids is excreted daily [80], which represents a significant investment of water resources. While some of these enzymes continue to be produced during dehydration anorexia, their quantity is likely decreased. As a direct result, the reduction of oral food intake may result in significant water savings, which at least in the short term with the presence of sufficient glycogen stores (see below), should be related to a reduction in water use and electrolyte derangement and therefore enhanced survival.

With limited food intake, there is no external source of glucose, nevertheless, we observed that glucose levels were still maintained [20] which is critical for surviving dehydration. This can be achieved by enhancing glycogenolysis and gluconeogenesis [28, 54, 55], responses that are designed to maintain blood glucose levels during fasting or starvation. In mammals, there is a vasopressin receptor in the liver that when bound activates gluconeogenesis [81], with secondary contributions from the kidney [82]. This suggests a successful two-pronged response to vasopressin secretion, such as vasoconstriction and reabsorption of water from kidney as well as a role in feeding behavior and energy balance. In the current study, we identified the candidate genes PCK1 (DE in the GI, liver, and kidney), a main control point for the regulation of gluconeogenesis [83] as well as CYP2E1 (DE in the hypothalamus, GI, and lung), a gene induced by starvation with products involved in gluconeogenesis [84, 85], as upregulated in water-deprived animals. Further, these genes were assigned to a significant module in all five tissues and defined as an outlier in the CCA analysis (Table 4). Further, the insulin resistance KEGG pathway [hsa04931 [69, 70] involves various mechanisms contributing to altered glucose metabolism and insulin responsiveness. Notably, all genes but two genes in hsa04931 are significantly up regulated in dehydrated animals, all of which are upstream of PCK1 are in the pathway. Two of these important pathway genes (AGT and SLC27A2, Table 4) were found in our consensus gene set. This upregulation is likely a compensatory response to ensure the maintenance of blood glucose levels in the absence of dietary intake. Lastly, INMT (Table 4), a gene in the Tryptophan metabolism pathway [hsa00390], [69, 70] is upstream of the Glycolysis/gluconeogenesis pathway [map00010], [69, 70], suggesting complex precursor genes and gene pathways may be contributing to glucose homeostasis as well. The significant changes in expression observed in these genes underscore their role in modulating glucose metabolism, shedding light on how organisms adjust to sustain vital glucose levels during periods of reduced food intake and water deprivation.

During water deprivation, a key challenge is to cope with altered fluid balance and maintain effective blood circulation and nutrient delivery under water-deficient conditions. This process can involve the formation of new blood vessels or the remodeling of existing ones to optimize perfusion, addressing the challenges posed by reduced fluid availability and potential hemoconcentration. We identified two genes from our consensus gene set that are mediators of new blood vessel formation and morphogenesis (ANG, CASZ, Table 4). Prior research has indicated that persistent activation of neuronal systems could modify local blood circulation through angiogenesis. Specifically, rats reared in complex environments had increased capillary density in the visual cortex [45]. Long-term motor activity has been documented to induce the development of new blood vessels within the cerebellar cortex [86] and primary motor cortex [46], and hyperosmotic stimuli, similar to what experimental animals experience, has been shown to modify the vasculature and induced reversible angiogenesis throughout the hypothalamic nuclei [44]. Furthermore, Alim et al. [47] observed seasonal differences in vascularization, showing blood capillaries were thicker during the winter, suggesting less diffusion across the membrane, compared summer in the dromedary camel.

To maintain blood flow and prevent the formation of clots that could impede circulation, downregulation of coagulation factors during periods of reduced fluid intake could be a protective mechanism against the potential risks associated with increased blood viscosity due to hemoconcentration. As dehydration leads to reduced water content in the blood, resulting in thicker blood consistency there is a risk for blood clot formation. four genes in our consensus gene list involved in blood composition and concentration (ALB, HP, SERPINA6, SERPINF2, Table 4). This could indicate several things: (1) To mitigate the impact of hemoconcentration, genes responsible for clotting factors are downregulated to reduce the risk associated with increased blood viscosity as thicker blood is more prone to clotting. (2) to decrease clotting during body temperature dysregulation. Blumstein and MacManes (20) discuss heterothermy as a mechanism for substantial energy and water savings. In addition to this, relative hypothermia can result in inactivation of coagulation enzymes and/or platelet adhesion defect [87], suggesting further mechanisms for maintaining perfusion and limit clotting. However, it’s important to note that the timing of expression poses a potential concern regarding correlations. Decreased body temperature was measured in Blumstein and MacManes (20) during the dark phase of the experiment and tissues were collected during the light phase of the experiment. Using a single time point snapshot for each experimental condition might overlook certain interactions, given the lag between the expression of a gene (such as a transcription factor) and the expression of downstream effectors.

Tissue specific responses

When analyzing the RNAseq data across all the tissues, we observed that, as expected, samples of the same tissue type clustered together regardless of the water-access treatment, suggesting that the signature of tissues-specific gene expression overpowers the signature of experimental treatment. During individual tissue analysis, the degree of sample separation in PCA space (Supplemental File 2) and CCA space (Fig. 3) and number of differentially expressed genes suggests tissue specific responses. While it is well known that the hypothalamus is the central regulating unit in the brain for maintenance of energy homeostasis [88]. However, few genes were identified as differentially expressed in the hypothalamus (Fig. 1), while our other analyses, WGCNA and CCA, uncovered genes and GO terms that responded to water deprivation. This suggests that the hypothalamus may not be well suited for bulk RNAseq studies due to the heterogeneity of the tissue. Future studies, particularly those using single cell methods [89,90,91] may further clarify the role that the hypothalamus plays in the overall response to dehydration.

Conclusion

Here, we highlight the intricate mechanisms involved in regulating water balance in the desert cactus mouse, P. eremicus. Our emphasis on whole-organism physiological and multi-tissue (kidney, GI, hypothalamus, liver, and lung) gene expression analysis within a simulated desert environment allowed us to achieve an understanding of genomic mechanisms of water homeostasis. Previous genome scan studies [1, 3, 4, 13, 72] have identified many of the same processes we found (i.e., metabolic processes, renal system development, immune response, and starvation response) however, our study design allows us to further interpret function because we were able to both described the tissue specific location and link the processes to physiological measurements.

At a whole-organismal scale, we observed a robust response of the renin-angiotensin-aldosterone system (RAAS) in dehydrated cactus mice, with upregulation of AGT in all five tissues as well as upregulation of other pathway genes. Additionally, the compensatory action of PCK1 was activated in all tissues, further supporting reduced food intake during water deprivation and underscores the body’s adaptive response. However, despite efforts to maintain blood volume, hemoconcentration still occurs, but in response there was a downregulation of genes responsible for coagulation (e.g., SERPINF2) as a protective measure against blood clotting in all five tissues, a gene with the major role of regulating the blood clotting pathway. The consequential thickened blood consistency poses challenges to effective blood flow through vessels, compelling the body to initiate the construction of additional vessels to enhance blood movement, further supported by the upregulation of ANG in all five tissues, further illustrating the complex interplay of regulatory processes in response to fluid balance disturbances.