Natural disaster and immunological aging in a nonhuman primate.

Weather-related disasters are increasing in frequency and severity, leaving survivors to cope with ensuing mental, financial, and physical hardships. This adversity can exacerbate existing morbidities, trigger new ones, and increase the risk of mortality-features that are also characteristic of advanced age-inviting the hypothesis that extreme weather events may accelerate aging. To test this idea, we examined the impact of Hurricane Maria and its aftermath on immune cell gene expression in large, age-matched, cross-sectional samples from free-ranging rhesus macaques (Macaca mulatta) living on an isolated island. A cross section of macaques was sampled 1 to 4 y before (n = 435) and 1 y after (n = 108) the hurricane. Hurricane Maria was significantly associated with differential expression of 4% of immune-cell-expressed genes, and these effects were correlated with age-associated alterations in gene expression. We further found that individuals exposed to the hurricane had a gene expression profile that was, on average, 1.96 y older than individuals that were not-roughly equivalent to an increase in 7 to 8 y of a human life. Living through an intense hurricane and its aftermath was associated with expression of key immune genes, dysregulated proteostasis networks, and greater expression of inflammatory immune cell-specific marker genes. Together, our findings illuminate potential mechanisms through which the adversity unleashed by extreme weather and potentially other natural disasters might become biologically embedded, accelerate age-related molecular immune phenotypes, and ultimately contribute to earlier onset of disease and death.

Survivors of extreme adverse events, such as natural disasters, have increased incidence of cardiovascular disease and chronic low-grade inflammation (1 –6). At the molecular level, severe hardship can also alter immune cell gene expression (7 –9) and advance hallmarks of aging such as prematurely aging T cell populations (10) and accelerating epigenetic aging (11 –13). While everyone ages, not everyone ages at the same rate, and individuals of the same chronological age can have dramatic differences in the onset or severity of age-related diseases (14 –16). Characteristic changes in age-associated phenotypes are known as biological aging. The difference between biological and chronological age (i.e., age acceleration or deceleration) is associated with morbidity and mortality in humans (17, 18), can be manipulated through experimental interventions in nonhuman primates and mice (19, 20), and can change throughout the life course due to environmental exposures (e.g., smoking) (21, 22). Age acceleration can be measured using biological clocks, which can quantify the difference between biological and chronological age (17, 18). The prevalence of shared morbidities as a result of extreme adversity and aging invites the hypothesis that exposures to adversity may accelerate aging (11, 12, 23, 24). Indeed, adversity associated with traumatic stress (e.g., war) or environmental disasters, like wildfire smoke, can increase inflammatory markers, indicating that premature immune aging may be a particularly salient mechanism by which disasters translate into disease (7, 12, 25 –27). At the molecular level, environmental conditions can induce persistent gene regulatory changes that affect gene transcription for many years, a mechanism through which adversity and other experiences may become biologically embedded and mechanistically explain lasting immune changes (28 –33). Despite these tantalizing links, prior studies of adversity and aging acceleration have typically lacked biological measures in the same population prior to the adverse event and almost exclusively investigated violence- and deprivation-related adversity (e.g., war, abuse). Thus, whether aged molecular phenotypes extend to natural disasters like extreme weather events (e.g., hurricanes, floods, tornados) remains unknown.

Hurricane Maria, the most destructive hurricane recorded in the history of Puerto Rico, made landfall in Puerto Rico on September 20, 2017, as a category 4 hurricane. One kilometer off the southeastern coast of Puerto Rico lies Cayo Santiago, a 15.2-ha island home to a population of 1,800 free-ranging rhesus macaques that have been studied for decades. Cayo Santiago was the first part of Puerto Rico hit and thus bore the full force of Hurricane Maria. The storm severely damaged homes and infrastructure across the Puerto Rican mainland and destroyed most of the vegetation, freshwater cisterns, and research structures on Cayo Santiago (Fig. 1 ) (34). The Cayo Santiago Field Station is an ideal system for studying the immune consequences of aging and environmentally induced adversity for several reasons. First, rhesus macaques share many behavioral and biological features with humans, including aging phenotypes. However, the rhesus macaque lifespan is approximately one-quarter of the human lifespan (35), permitting effective sampling of a significant portion of the aging process in just a few years. Third, macaques on Cayo Santiago are provisioned with food and water—even in the immediate aftermath of the hurricane—thus controlling for nutritional disruptions that might additionally affect the aging process (36 –38).

Fig. 1.

Hurricane Maria caused severe destruction to Cayo Santiago and was associated with altered immune cell gene regulation. (A) Hurricane Maria struck the island of Cayo Santiago on September 20, 2017, causing severe devastation to the vegetation and substantial remodeling of the environment. Aerial photographs show Cayo Santiago in August 2008 (Left) and January 2020 (Right). Photos reused with permission of Joyce Cohen, WOM Productions, and Michelle Skrabut La Pierre. (B) Vegetation, as measured by NDVI, decreased by 63% due to the storm (t test, P = 3.7 × 10−25). (C) Temperature in areas that lost vegetation due to Hurricane Maria (i.e., de-vegetated) were significantly warmer than areas still vegetated after the storm (t test, P < 1 × 10−10). (D) Peripheral whole blood was collected cross-sectionally from the rhesus macaque population on Cayo Santiago island approximately annually 1 to 4 y before (n = 435) and 1 y after (n = 108) Hurricane Maria struck the island. Animals sampled posthurricane were an average of 0.89 y younger than animals sampled prehurricane (median age prehurricane = 6.94, median age posthurricane = 6.05, t test of means P value = 1.16 × 10−3). (E) PC3 of global gene expression captures the shared effects between increasing age (β = −0.19, P = 8.96 × 10−8) and (F) exposure to Hurricane Maria (β = −1.79, P = 9.12 × 10−5) on immune gene expression. PC3 explains 3.4% of the overall variance in gene expression.

To test how Hurricane Maria influenced immune cell gene regulation and aging, we leveraged our collection of peripheral blood samples and detailed demographic data ( Dataset S1) from age-matched cross-sectionally sampled subsets of the Cayo Santiago rhesus macaque population annually in the 4 y prior to Hurricane Maria (n = 435) and 1 y after (n = 108) Hurricane Maria (Fig. 1). We hypothesized that exposure to the hurricane would recapitulate molecular changes associated with the natural process of aging.

Results

Age-Associated Immune Cell Gene Expression Patterns.

We first quantified how chronological age was associated with immune cell gene expression. As expected, chronological age was strongly associated with immune cell gene regulation. Age was significantly associated with the first and third principal components of gene expression, which explained 55.2% and 3.4% of the variance, respectively (PC1: age β = −0.39, P = 0.005, ; PC3: age β = −0.19, P = 8.96 × 10−8, Fig. 1). At the level of individual genes, we found that 16% of the 7,009 detectably expressed genes were significantly associated with chronological age (n genes = 1,131, false discovery rate [FDR] < 10%, Dataset S2, in a model controlling for Hurricane Maria exposure, sex, and RNA quality [RQN], n = 543). Genes more highly expressed in older animals were enriched for biological processes associated with inflammation (Benjamini-Hochberg corrected pBH = 5.72 × 10−5) and innate immune system activity (e.g., negative regulation of lymphocyte mediated immunity, pBH = 0.005) ( Dataset S3). Transcription factor (TF) binding site analysis identified that these genes were putatively under control of TFs in the bZip protein family (e.g., AP-1, Fosl2, JunB) as well as enriched for components of TFs NR5A2 and NFkB (FDR < 5%), which are both implicated in cytokine production and inflammatory diseases (39, 40) ( Dataset S4). By contrast, genes more highly expressed in younger animals were enriched for biological processes associated with translation (pBH = 1.86 × 10−23) and immunoglobulin production (pBH = 0.108) ( Dataset S5). We did not find significant evidence for any sex differences in aging. We also examined nonlinear age-related changes in gene expression using autoregressive integrated moving average (ARIMA) models (41) and detected 1,246 genes with nonzero trends that clustered into five aging trajectories ( ). Clusters 1, 2, and 3 changed linearly across aging (i.e., increasing, decreasing) and 62.9% (n = 656) of genes in these clusters overlapped with genes detected by our linear modeling approach. Notably, genes in clusters 4 and 5 exhibited nonlinear trajectories with increasing age. Genes in cluster 4 followed an inverted U-shape trajectory across the lifespan (increase in midlife and decrease in late life), while cluster 5 genes exhibited a U-shaped trajectory (decrease in midlife and increase in late life). Cluster 5, in particular, included several important immune genes such as CCL5 and CD4, which encode for the T cell chemoattractant CCL5, expressed by activated T cells, and protein CD4, which is critical in the antigen detection process for regulatory and helper T cell subsets, suggesting that expression of these key immune genes reaches a nadir in midlife but rapidly increases in late life.

Exposure to a Natural Disaster Is Associated with Altered Immune Cell Gene Expression.

We next evaluated whether experiencing Hurricane Maria was associated with differences in immune cell gene expression. Exposure to Hurricane Maria was significantly associated with PC3 (β = −1.79, P = 9.12 × 10−5) (Fig. 1) and, at the gene level, differential expression in 4% of genes (n = 260, FDR < 10%). We found no sex differences in the gene expression response to Hurricane Maria. Genes more highly expressed in animals sampled after Hurricane Maria (n = 54) were implicated in inflammation (e.g., positive regulation of IL-8 production, pBH = 0.052) ( Dataset S6). Genes with reduced expression in animals sampled after Hurricane Maria (n = 206) were involved in translation (pBH = 2.58 × 10−6), chaperone cofactor-dependent protein refolding (pBH = 5.75 × 10−2), and processes suggestive of a decrease in adaptive immune activity (e.g., regulation of T cell differentiation, pBH = 0.096) ( Dataset S7). Genes with reduced expression in animals sampled after Hurricane Maria were putatively controlled by HSF1 (FDR < 5%) ( Dataset S8), the central TF of the heat shock response (HSR). Together, down-regulation of translation and protein folding activity suggest a disruption to protein folding networks that promote proteostasis—the loss of which is a hallmark of aging (13, 42). Specifically, HSC70/HSPA8, a constitutively expressed component of the HSP70 family, exhibited the strongest association with exposure to Hurricane Maria. Expression of HSC70 was 2× lower after Hurricane Maria (β = −1.03, se β = 0.09, FDR = 1.02 × 10−24). Decreased HSC70 activity may exacerbate cardiovascular diseases (43, 44) and Alzheimer’s disease (45, 46), which are both age-related diseases that are more prevalent in survivors of severe adversity (1, 6). Further, among the 260 genes differentially expressed between animals sampled before and after Hurricane Maria, gene-gene coexpression networks were disproportionately disrupted for genes associated with the HSR. Seven gene pairs were differentially correlated before versus after Hurricane Maria (FDR < 10%). These pairs were composed of 11 unique genes, meaning that several genes were represented in multiple pairs. Five of these 11 genes were involved in the HSR, an overrepresentation compared to the background rate among hurricane-associated genes, underscoring that the HSR is particularly disrupted in individuals that experienced Hurricane Maria (binomial test, P = 6.97 × 10−4).

Exposure to a Natural Disaster Broadly Recapitulates Gene Expression Differences Associated with Aging.

We quantified the extent to which the effects of age and exposure to the hurricane on immune cell gene expression were similar in two ways. First, 40% of the 260 hurricane-associated genes were also significantly associated with age (n = 104), which is 2.5× more than expected by chance (Fisher’s exact test, odds ratio [OR] = 3.7, P = 4.06 × 10−21) (Fig. 2 and Dataset S2). Moreover, these 104 genes were significantly more likely to exhibit alterations in expression that were consistent with the prediction that exposure to the hurricane was akin to increases in chronological age. In other words, we found strong and significant concordance in the hurricane and aging effects in these 104 genes: 84.6% of the genes (n = 88) significantly associated with age and the hurricane had higher expression with both aging and exposure to Hurricane Maria or lower expression with both aging and hurricane exposure (binomial test, P = 3.18 × 10−13) (Fig. 2). Second, we quantified if the expression of genes was similarly associated with hurricane exposure and aging in both magnitude and direction. Across all genes, the effects of aging and exposure to Hurricane Maria were significantly positively correlated (r = 0.23, P = 1.33 × 10−84) (Fig. 2). The strength of this correlation increased more than twofold when we limited to genes that were more strongly associated with aging and the hurricane (genes with FDR < 10% for both age and Hurricane Maria: r = 0.58, P = 1.05 × 10−10) ( ). Importantly, effects of the hurricane and aging persisted when we controlled for potential covariates of body condition, wound presence, and sampling condition (see section ).

Fig. 2.

Similar effects of Hurricane Maria and aging on immune cell gene expression. The standardized effect sizes of experiencing Hurricane Maria on immune gene expression are positively correlated with standardized effect sizes of aging on immune gene expression (r = 0.23, P = 1.33 × 10−84). Genes that are significantly associated (FDR < 10%) with aging and Hurricane Maria correlate even more strongly (r = 0.58, P = 1.05 × 10−10). Genes significantly associated with both effects are colored purple. Inset plots show schematics of the gene expression differences represented in each quadrant.

Genes with lower expression levels in samples from animals that experienced the hurricane and in older individuals suggested a disruption to protein folding processes (e.g., proteostasis, pBH = 5.94 × 10−5) and telomerase activity (pBH = 0.02) ( Dataset S9). Genes with reduced expression in samples after Hurricane Maria and with older age were enriched for the TF HSF1 binding site (FDR = 0.005) ( Dataset S10). Further implicating the HSR pathway, we found that genes coding for molecular chaperones that aid in protein folding and protein degradation (i.e., ubiquitination) overwhelmingly had lower expression in samples from older individuals and those that experienced Hurricane Maria (Fig. 3). In fact, genes involved in molecular chaperoning and ubiquitination were among the most perturbed (Fig. 3 and ), comprising 9 of the top 10 most affected by the hurricane ( Dataset S2). Decreased molecular chaperone activity leads to protein misfolding and aggregation and decreased clearing of disabled proteins from the cell via ubiquitination and is implicated in many age-related diseases (Fig. 3 and ) (47 –49). Together, these data suggest that exposure to Hurricane Maria may increase protein misfolding and contribute to the onset and progression of age-associated diseases.

Fig. 3.

Hurricane Maria disrupted genes involved in protein homeostasis and was associated with accelerated gene expression aging. (A) Aging and experiencing Hurricane Maria were associated with reduced expression of key molecular chaperone proteins. Genes shown were differentially expressed with aging and/or hurricane exposure (FDR < 10%). (B) Molecular chaperone proteins in the protein homeostasis pathway that were significantly associated with the hurricane and aging. These proteins promote proper folding of client proteins and clear terminally misfolded proteins from the cell. (C) Predicted biological ages based on gene expression profiles for prehurricane (green) and posthurricane (brown) samples. Posthurricane samples have a predicted biological age that is, on average, 1.96 y greater than their chronological age (P = 0.0003).

Only 16 genes were associated with the hurricane and aging in opposing directions (meaning that they were up-regulated with aging and down-regulated following Hurricane Maria; Fig. 2). Genes in this group included IRF1, CCL5, PRF1, and GSTP1, which are integral to the aging process, yet may represent a signature of an environmental stress-specific response that differs in direction from age-related processes.

Lastly, we tested whether animals exposed to Hurricane Maria had immune cell expression profiles consistent with signatures of accelerated aging. To do so, we used human transcriptional age predictors that have been successfully applied by our group to predict age from gene expression in rhesus macaques (35, 50). We hypothesized that samples after the hurricane would have older biological age predictions, which might be indicative of accelerated biological aging in these animals. Reproductively mature adult animals that experienced Hurricane Maria had gene expression profiles that were, on average, 1.96 y biologically older than adults sampled before Hurricane Maria (in a linear model of predicted age controlling for chronological age and hurricane exposure, P = 3.0 × 10−4; Fig. 3). When scaled to the human lifespan, immune gene expression of adults who experienced the hurricane was equivalent to an age acceleration of 7 to 8 y.

Aging and Exposure to Natural Disaster Are Both Associated with Higher Expression of Innate Immune Cell Marker Genes.

Immune cell composition changes across the lifespan and could therefore be reflected in the transcriptome. To examine age or hurricane-associated differences in immune cell composition, we used immune cell marker genes identified from single-cell RNA sequencing ( Dataset S11). Overall, cell-specific markers of canonical anti-inflammatory immune cells had lower expression in older individuals and those that experienced the hurricane, while markers of myeloid-derived cells associated with the inflammatory response had higher expression (Fig. 4 ). Our findings are consistent with known changes in cell populations across aging and also demonstrate that experiencing Hurricane Maria recapitulated age-associated differences at the level of immune cell composition. Specifically, granulocyte, classical (CD14+) monocyte, nonclassical monocyte (CD14-), NK cell, and cytotoxic T cell marker genes were more highly expressed with aging, while helper T cell and B cell marker genes were more lowly expressed (pBH < 0.05). Animals exposed to Hurricane Maria had greater expression of classical monocyte marker genes and lower expression of helper T cell marker genes than individuals who did not experience the hurricane (pBH < 0.05) (Fig. 4).

Fig. 4.

Aging and experiencing Hurricane Maria were associated with broadly similar immune cell marker gene expression. Shaded points represent marker genes that were significantly differentially expressed (FDR < 10%) across aging (A) and Hurricane Maria (B). Open points represent marker genes that were not significantly associated with the effect of aging or hurricane exposure. Asterisks denote cell types for which differentially expressed marker genes were either positively or negatively enriched (pBH < 0.05).

Discussion

Here, we leveraged a rare natural experiment on a cross section of animals sampled either before or after experiencing a major hurricane. Our findings suggest that differences in immune cell gene expression in individuals exposed to an extreme natural disaster were in many ways similar to the effects of the natural aging process. Further, we observed evidence for accelerated biological aging in samples collected after animals experienced Hurricane Maria. This adds yet another factor to the growing list of adverse experiences that can alter immune gene expression, disrupt protein homeostasis networks, and appear to accelerate markers of biological aging. Importantly, we identify a critical mechanism—immune cell gene regulation—that may explain how adversity, specifically in the context of natural disasters, may ultimately “get under the skin” to drive age-associated disease onset and progression. Indeed, the leukocyte transcriptome was, on average, 1.96 y older in animals that experienced Hurricane Maria compared to those that did not, corresponding to predicted biological age acceleration of ∼7 to 8 human years. These alterations were also not short-lived: we detected the hurricane-associated differences in immune cell gene expression 1 y after the event. Together, our findings support the idea that the peripheral immune system is altered following an extreme adverse event and may help to explain how natural disasters increase morbidity and mortality rates.

The most salient result associated with aging and in animals that experienced Hurricane Maria were lower levels of expression of genes associated with the HSR and molecular chaperone activity. Canonically, the HSR is activated in response to heat stress, ultraviolet radiation, infection, and other acute environmental stressors, yet we observed substantial down-regulation in HSR genes despite the extensive loss of shade and temperature increase across Cayo Santiago. Activation of the HSR induces molecular chaperones to assist in protein folding and guide misfolded and aggregate proteins toward degradation pathways (51, 52). Components of the HSR involved in protein homeostasis, known as proteostasis networks, decrease throughout aging (51), which leads to greater incidence of the protein misfolding and aggregation implicated in many diseases (e.g., Huntington’s disease, Alzheimer’s disease, and amyotrophic lateral sclerosis). Our results indicate that aging and experiencing Hurricane Maria were both associated with lower levels of molecular chaperone gene expression. Thus, the hurricane-associated down-regulation of the HSR may indicate dysregulation of proteostasis networks and contribute to observed disease susceptibility in survivors of extreme adversity.

Notably, several components of the HSR are also directly involved in the glucocorticoid response (53 –55), which is activated in response to external stressors and well documented to be disrupted by adverse experiences (56). Specifically, greater expression of FKBP5 following Hurricane Maria may indicate increased FKBP51, which is implicated in aberrant stress responses and glucocorticoid resistance following adversity (53 –58). These findings add to the increasing evidence that the aberrant GC responses frequently observed in human survivors of extreme adversity may be mechanistically driven by disruptions to FKBP5 at the genomic or epigenomic level and extend to survivors of natural disasters (29, 59, 60).

Changes in gene expression at the bulk tissue level are also sensitive to changes in the relative cell proportions in the tissue. We leveraged this fact and found that adaptive immune-associated cell marker genes (e.g., helper T cells) were more lowly expressed, while inflammatory innate immune cell marker genes (e.g., classical monocytes) were more highly expressed in samples from animals exposed to Hurricane Maria. Additionally, we observed a general pattern of greater expression of markers of proinflammatory cells in animals sampled after Hurricane Maria (Fig. 4). The fact that hurricane exposure broadly recapitulates age-associated changes in key immune cell marker genes adds to recent studies in humans finding “aged” T cell composition (10). Chronic proinflammatory immune activity is characteristic of older populations, leading to the suggestion that up-regulation of proinflammatory immune cells, and their associated activity (e.g., inflammatory cytokines), may be a mechanism by which severe adversities contribute to morbidity and premature mortality (3, 8, 10, 61).

Our study centered on understanding the effects of a major natural disaster, which included both the effects of the hurricane itself as well as its aftermath. It is difficult, and potentially not necessary, to disentangle the initial event from the aftermath because the devastating repercussions of severe natural disasters contribute to their extreme nature and are enduring sources of adversity. Due to the opportunistic nature of our study, the available samples from 2018 were from a single social group. It is therefore possible that the animals sampled after the hurricane differed from the prehurricane samples in ways that may confound our results. We thus carefully examined all potential confounds we could measure (body condition, wound presence, and sampling conditions) and found that the hurricane and aging effects persisted despite any influence from these covariates (see section ).

Our study relied on opportunistic sampling surrounding an extremely unpredictable event, which is accompanied by some limitations in our study design. The cross-sectional, rather than longitudinal, nature of our study means that we were unable to measure aging rates within the same individuals either before or after the hurricane. This design therefore limits our ability to disentangle aging effects from cohort effects. However, the age-associated results we identified recapitulate those found in studies of humans and other closely related species (35), suggesting that any cohort effects are likely small. Importantly, because age and hurricane exposure are not correlated in our study, we are able to disentangle their effects on immune gene expression as animals before and after the hurricane were drawn from similar age distributions—with the only difference being that animals sampled in 2018 experienced Hurricane Maria. Nevertheless, future work would benefit from longitudinal studies that can more accurately quantify within-individual changes in age and in the face of a natural disaster.

Notably, we did not find evidence that the observed immune disruption in animals that experienced Hurricane Maria varied by sex, contrary to expectations from human studies (62 –64). However, it has been suggested that heightened morbidity and mortality observed in human female populations after natural disasters may be due to sociocultural constraints rather than biological differences (65, 66). Since we did not find that sex moderated immune responses to environmental devastation in rhesus macaques, our findings support the idea that sociocultural constraints placed on human female populations, rather than conserved biological differences, contribute to the observed health disparity human female populations face after natural disasters.

We observed heterogeneity in animals’ gene expression after Hurricane Maria that was not explained by any variables we tested. Lack of social support both following adversity and throughout aging is associated with increased disease risk (67 –69). Additionally, it is increasingly recognized that gene expression is highly responsive to social constraints and that social stress can drive broadly similar immune modifications to the stress response in humans and nonhuman primates (70 –72). The effects of sociality can be highly context-specific but, broadly, socially stressed individuals would be expected to have an exacerbated response to extreme adversity. However, in this system, animals that were most socially isolated before the storm adjusted their social milieu the most after the storm (34), suggesting social flexibility might buffer these individuals from the effects of the hurricane yet obscure predictions about which individuals might be most socially stressed. Future studies should investigate the extent to which social factors, such as social dominance, integration, and flexibility, influence immune gene regulation and cell composition following natural disaster.

In conclusion, we find support for the hypothesis that environmental adversity and the natural aging process were associated with broadly similar disruptions in peripheral immune regulation. While the biological costs of adversity are becoming increasingly clear, our findings encourage efforts to develop a thorough understanding of the modifiers of aging and adversity, which may provide rapidly realized mitigation of detrimental immune system effects for populations that have experienced extreme adversity.

Materials and Methods

Study Subjects and Experimental Design.

Data were collected from rhesus macaques (Macaca mulatta) on the island of Cayo Santiago, a long-term field site where animals have been continuously monitored between 1938 and 2018. Animals are individually identified by tattoos and ear notches. Whole blood samples, body weight measurements, and identification of the presence of wounds were collected during the annual trap-and-release period. Trapping of animals on Cayo Santiago occurs annually, ∼1 to 3 mo before the breeding season, between October and February, to minimize interference with the reproductive season. Specifically, before the annual trap-and-release period in years 2013, 2014, 2015, 2016, and 2018, age and sex-matched individuals were selected for biological sampling as part of other ongoing studies. When Hurricane Maria struck in 2017, we took advantage of the rare opportunity to sample part of the population in 2018 when conditions allowed. Thus, while different animals were sampled each year, they were broadly age and sex-matched across each year and across the hurricane (Fig. 1). Additionally, because natural disasters are difficult to predict, to have pre- and posthurricane samples, this project took advantage of samples collected for other ongoing projects.

Although Hurricane Maria was an extreme weather event, mortality in the population was relatively limited in the wake of the storm—only 2.75% of the population died in the immediate aftermath of the storm ( ). Notably, there was no differential survival based on sex or age in the year after the storm, as the average age of death within 1 y of Hurricane Maria did not differ from the average age of death in the 2 decades prior to the storm (P = 0.13). This supports recent findings that no overall increase in mortality has been found after the three previous major hurricanes to hit Cayo Santiago (Hurricanes Hugo in 1989, Georges in 1998, and Maria in 2017) (73).

Blood Sampling.

Whole blood was drawn from sedated rhesus macaques by veterinary staff. A 2.5 mL blood sample was collected in PAXgene Blood RNA Tube (PreAnalytiX GmbH) and stored at room temperature for 2 to 4 h before being transferred to 4 °C for 24 h and finally to −80 °C for long-term storage until RNA extraction. Before Hurricane Maria, we collected 435 blood samples from 357 animals from December 2013 to December 2016 (median age of 6.94 y; n females = 174, n males = 261) and sampled 108 animals after Hurricane Maria from October to December 2018 (median age = 6.05 y; n females = 57, n males = 51). Detailed metadata for each sample can be found in Dataset S1. As part of ongoing population management, a subset of the animals sampled in 2016 were captured on Cayo Santiago and transported to the Sabana Seca Field Station (n = 95), where the same procedure was performed. Due to the limited infrastructure on the island of Cayo Santiago and as part of the same population management, animals sampled in 2018 (n = 108) were also sampled at Sabana Seca.

RNA Extraction, Sequencing, and Data Preprocessing.

RNA for 499 samples was extracted using the MagMAX for Stabilized Blood Tubes RNA Isolation Kits (Thermo Fisher), and RNA from 48 samples was extracted using the PAXgene Blood RNA Kit IVD (Qiagen) following the standard protocols for each procedure. RQN was quantified using AATI Fragment Analyzer. RNA sequencing libraries were prepared using 50 ng of total RNA and following a recently developed 3′-biased protocol, TM3′SEq. (74). Libraries were amplified with 16 PCR cycles. All other procedures followed the published protocol or manufacturer recommendations. Libraries were combined in equimolar quantities and sequenced on an Illumina NovaSeq S2 flowcell (R1 was 25 bp and R2 was 80 bp to capture the complementary DNA [cDNA] transcript) to an average read depth of 3.5 million reads per sample. We mapped cDNA reads to the M. mulatta reference assembly Mmul_10 using kallisto (75) (average mapping rate = 71.1%).

Read Count Normalization.

Prior to normalization, we removed reads mapping to seven genes encoding hemoglobin and ribosomal RNA subunits, which together comprised 54% of all mapped reads. We aggregated expression levels at the level of the gene since 3′-seq cannot accurately assign reads to all alternative transcripts of a gene. We removed lowly expressed genes with fewer than 3.28 transcripts per million, the value that most evenly split the distributions of lowly expressed genes and highly expressed genes, resulting in 7,009 detectibly expressed genes for downstream analyses (https://github.com/mwatowich/Natural-disaster-and-immunological-aging-in-a-nonhuman-primate). Read counts were normalized using the voom function from the limma package in the R environment (76), which normalizes counts overall to account for between-sample variation by estimating the mean-variance relationship of the log-counts among all samples.

Vegetation and Temperature Analyses.

We measured how the vegetation on Cayo Santiago changed following Hurricane Maria by quantifying normalized difference vegetation index (NDVI) from Landsat 8 satellite images available on Sentinel-hub EO-Browser from ∼2 y before and 2 y after Hurricane Maria. NDVI measures the density of greenness in a given area and is the most commonly used remote sensing measure of greenness (77). To elucidate how temperatures across the island had been altered by the extensive vegetation loss due to Hurricane Maria, we installed remote temperature sensors 0 to 2 m from the ground in areas with 1) vegetation cover before and after Hurricane Maria (“vegetated”) and 2) vegetation cover before Hurricane Maria that lost vegetation cover after the storm (“de-vegetated”) (Fig. 1). Temperature sensors were installed after Hurricane Maria and were active from June to August 2018 and April 2019 to present. Detailed methods are described in Testard et al. (34).

Modeling the Effects of Age and Hurricane Maria on Gene Expression.

We carried out statistical analyses in the R version 4.0.2 environment unless otherwise specified. To evaluate global trends in gene expression variation, we performed a principal component analysis on the correlation matrix of normalized gene expression values. We used R package EMMREML (78) to employ a mixed modeling approach to quantify the effects of age, sex, and exposure to Hurricane Maria on gene expression while controlling for effects of RQN, which was quantified during RNA library preparation, and relatedness between individuals. Kinship was estimated using ANGSD and ngsRelate (79) from mapped 3′seq bam files. Data from all 543 samples were included in these models. Genes that passed an FDR threshold of 10% were considered to be differentially expressed. q values were calculated using the R qvalue package (80).

Modeling Interactions between Sex and Age and Hurricane Maria on Gene Expression.

As prior research has identified sex differences in how the immune system ages (81) and how different sexes respond to adversity (62, 65, 82), we tested whether the sexes differed in their responses to aging and Hurricane Maria by modeling an interaction between sex and chronological age and, in a separate model, sex and hurricane exposure. Only three genes passed our FDR of 10% for an interaction between sex and age, and no interaction effects of sex and hurricane exposure passed our FDR of 10%.

Temporal Trends in Age-Related Gene Expression.

To further investigate how gene expression changed with increased age, we used an ARIMA model adapted from Márquez et al. (81). This approach determines the best-fitting ARIMA model for each gene as a function of age and selects the models that exhibit gene expression trajectories that are significantly different from random fluctuations across age (41). We controlled for batch effects of RQN on gene expression by using a partial residual gene expression matrix calculated by subtracting the residuals of RQN (quantified by our linear modeling approach) from our normalized gene expression matrix. We limited our analysis to adults (>3 y, n = 505), as infant gene expression often strongly differed from adults, decreasing our ability to effectively group genes based on temporal trajectories occurring throughout the majority of animals’ lifespan. Finally, we grouped gene trajectories by calculating a distance matrix of the correlation of predicted values across aging.

Enrichment Analyses.

To identify biological processes that were enriched in genes differentially expressed with chronological age and exposure to Hurricane Maria, we conducted Gene Ontology (GO) enrichment analysis using the R package topGO (83). We used the weight01 algorithm to conduct Fisher’s exact tests. To investigate GO enrichment in genes associated with age, we compared genes that significantly increased in expression across age (from the results of our linear model) to a background of all other genes. We then repeated this analysis for differentially expressed genes that had significantly lower expression in older individuals. We repeated these analyses for genes significantly up- and down-regulated in samples collected before and after Hurricane Maria. To characterize enrichment of biological functions in genes associated with both age and Hurricane Maria, we compared differentially expressed genes that had i) greater expression in older and hurricane-exposed individuals, ii) lower expression in older and hurricane-exposed individuals ( Dataset S9), and iii) greater expression in older individuals but lower expression in hurricane-exposed individuals against a background of all other genes in the dataset. There were no genes that passed an FDR threshold of 10% that were more lowly expressed in older and more highly expressed in hurricane-exposed individuals; thus, we did not perform an enrichment test on this subset of genes.

To characterize putative gene regulatory mechanisms, we tested for enrichment of TF binding motifs within 2 kb upstream and downstream of TF start sites of genes significantly associated with aging and Hurricane Maria using the program HOMER (84). We use all other genes in our dataset as the background set of genes for this test. We searched for known vertebrate TF binding motifs and report all TF motifs that passed an FDR threshold of 10%. We tested the same subsets of genes differentially expressed with aging, Hurricane Maria, and shared effects as our GO enrichment analyses.

Transcriptomic Clock.

To test whether samples collected after Hurricane Maria were predicted to have accelerated biological ages, we predicted ages of the animals in our dataset from their normalized gene expression profiles using human transcriptomic age predictors (50). This approach predicts transcriptomic age via weighted scores for each gene multiplied by scaled normalized gene expression data. To calculate this metric, we used genes that overlapped between the Peters dataset and our dataset (matched using Ensembl homolog annotations from biomaRt) (85, 86). Peters et al. performed a meta-analysis using a discovery and replication dataset to produce more strongly supported markers (50). The number of overlapping genes between our dataset and the meta-analysis was limited, and we therefore used markers from their discovery dataset that overlapped with genes in our dataset that were highly associated with aging (FDR < 1%) for a total of 321 genes.

Using the transcriptomic age predictor and equation provided by Peters ​et al​. (50), we predicted the chronological ages of animals in our dataset from their gene expression data. We estimated the predictor (​Z)​ using Eq. on our normalized gene expression matrix.

Predictions of biological age did not control for sex, as we did not detect sex effects on gene expression. After predicting chronological age from normalized gene expression, we scaled the predictions to the age distribution in our sample population using Eq. .

There was a significant correlation between scaled predictions and chronological age among the prehurricane samples (r = 0.42, P = 7.12 × 10−19). To determine the extent to which experiencing Hurricane Maria moderated predicted age, we modeled the scaled predictions as a function of chronological age and exposure to Hurricane Maria. We first predicted transcriptomic ages of sexually mature adults (>1 y) for similar reasons as in our ARIMA modeling approach and found that biological age was 1.96 y older on average in samples after Hurricane Maria (P = 0.0003). This effect was slightly weaker but still significant when we predicted biological age with all genes in our dataset that were associated with age at an FDR of <10% (n = 838; β = 1.64, P = 0.003). To test whether results held when we included infants, we performed the analysis again including samples from all animals (<1 y) and found that biological age was an average of 1.1 y older in samples from after Hurricane Maria (P = 0.03).

Testing Other Sources of Transcriptional Variation.

We tested whether body condition, presence of wounds, or sampling condition moderated the differential gene expression we observed in samples after Hurricane Maria, as we hypothesized that each may be a potential confound in our analysis. These effects were modeled using the previously described linear model (gene expression as a function of age, exposure to Hurricane Maria, sex, RQN, and estimated kinship) and including the additional covariate. None of these variables significantly altered the effects of the hurricane on gene expression.

We hypothesized that animals after the hurricane might have poorer body condition due to vegetation loss and adversity related to social changes following the hurricane; however, we found no difference in body condition in individuals sampled before or after the hurricane (Wilcoxon rank sum test: P = 0.09; ). Body condition was calculated using a loess regression of body weight controlling for age, sex, and location where weight was recorded (i.e., sampling locations with different scales). At the gene level, only two genes passed an FDR threshold of 10% with body condition.

We expected that the immune-related differences in gene expression we observed following Hurricane Maria may have been driven by wounds that might be more prevalent following the hurricane or might have healed more slowly if animals had other undetected infections. We tested whether the absence or presence of a wound (e.g., open cuts, broken fingers, etc.) affected gene expression and found that no genes passed a liberal false discovery rate of 20% (n animals with wound data = 203).

A subset of animals in our study were trapped on Cayo Santiago and transported to an inland facility where the biological sampling took place; this included 95 animals in 2016 that were part of an ongoing population management and all animals in 2018 because research infrastructure on Cayo Santiago was still destroyed at the time. Because of this potential confound with the hurricane, we tested whether effects of the hurricane were similar in the entire sample (n = 543) and in only off-island sampled animals from 2016 and 2018 (n = 203) using the same modeling framework controlling for age, sex, RQN, and kinship. The effect of the hurricane on gene expression was highly correlated between the smaller sample (off-island blood drawn, n = 203) and the main model (n = 543) among all genes (Pearson’s r = 0.54, P < 1 × 10−16) and those that passed an FDR of 15% for both models (n = 45, r = 0.75, P = 3.71 × 10−9). Further, hurricane and age-associated effects on gene expression were also still correlated in the same direction in the models with the reduced sample size (r = 0.17, P = 7.4 × 10−49) among all genes and those that passed an FDR of 15% for both effects (n = 57, r = 0.59, P = 1.52 × 10−6).

Characterizing Gene–Gene Relationships.

To characterize whether relationships between gene pairs were altered by Hurricane Maria, we quantified correlations between pairs of genes differentially expressed with Hurricane Maria using “correlation by individual level product” (87). We used gene expression values that controlled for residuals of RQN, sex, and aging. To test whether a given pair of genes was differentially correlated in samples before Hurricane Maria versus samples from after the hurricane, we first scaled gene expression residuals within each group (prehurricane, posthurricane) (87). Next, we multiplied normalized gene expression values for each gene in the pair and used the EMMREML package to fit this vector of products as a function of exposure to Hurricane Maria while controlling for kinship. To test whether genes that were in differentially correlated pairs across Hurricane Maria were more likely to be involved in the HSR than by chance, we tested for enrichment using a binomial test with a background rate of all genes with a GO term name involving “heat shock,” genes encoding for HSP subunits, and genes encoding for DNAJ subunits (17 of 230 hurricane-associated genes).

Marker Gene Enrichment Tests.

We characterized whether and how immune cell marker genes were associated with aging and Hurricane Maria exposure. We identified immune cell-specific marker genes from single-cell sequencing data of rhesus macaque peripheral blood mononuclear cells (PBMCs). Single-cell RNA sequencing data were generated from two adult female rhesus macaques from the large breeding colonies at Yerkes National Primate Research Center. For each individual, PBMCs were isolated with density gradient centrifugation and cells were captured on a 10× Genomics Chromium controller. Sequencing libraries were prepared using Chromium Single-Cell 3′ v2 chemistry according to the manufacturer’s specifications. Libraries were sequenced on an S2 flow cell on an Illumina NovaSeq 6000. Resulting reads were mapped to the MacaM v7 reference genome (88) (including mitochondrial and intronic sequences) with –count function in CellRanger v.3.0.2 (89). Individuals were demultiplexed based on standing genetic variation and doublets were removed with demuxlet (90) using the following parameters: –tag-group CB –tag-UMI UB –field GT –geno-error 0.05 –doublet-prior 0.076 –min-uniq 30 –min-snp 5. Count matrices were filtered to include only cells with at least 200 genes expressed and only genes expressed in at least 3 cells with Seurat v3 (91). The resulting count matrices included ∼14k genes and ∼2.3k cells per individual. Marker genes were identified using the Seurat (91) and monocle3 (92 –95) packages in the R environment. One hundred and forty-nine genes passed a specificity threshold of 0.1 and sensitivity threshold of 0.1 ( Dataset S11).

Because our bulk RNA sequencing data were generated from peripheral whole blood, which likely contains pertinent immune cells types beyond those in PBMCs, we also characterized granulocyte marker genes published in Palmer et al. (96), which identified canonical granulocyte marker genes from human peripheral blood samples. We expect that the majority of granulocyte marker genes in our dataset are specific to neutrophils because they comprise the vast majority of granulocytes in primates. We used a subsampling approach to test whether marker genes for each cell type changed significantly with increased aging or hurricane exposure. Specifically, we constructed a null distribution of median standardized effect sizes by resampling all standardized effect sizes and considered marker genes significantly perturbed if the marker genes’ median standardized beta fell outside of Bonferroni corrected CIs of 0.18% and 99.82%.