Neuroprotectants attenuate hypobaric hypoxia-induced brain injuries in cynomolgus monkeys.

Hypobaric hypoxia (HH) exposure can cause serious brain injury as well as life-threatening cerebral edema in severe cases. Previous studies on the mechanisms of HH-induced brain injury have been conducted primarily using non-primate animal models that are genetically distant to humans, thus hindering the development of disease treatment. Here, we report that cynomolgus monkeys ( Macaca fascicularis) exposed to acute HH developed human-like HH syndrome involving severe brain injury and abnormal behavior. Transcriptome profiling of white blood cells and brain tissue from monkeys exposed to increasing altitude revealed the central role of the HIF-1 and other novel signaling pathways, such as the vitamin D receptor (VDR) signaling pathway, in co-regulating HH-induced inflammation processes. We also observed profound transcriptomic alterations in brains after exposure to acute HH, including the activation of angiogenesis and impairment of aerobic respiration and protein folding processes, which likely underlie the pathological effects of HH-induced brain injury. Administration of progesterone (PROG) and steroid neuroprotectant 5α-androst-3β,5,6β-triol (TRIOL) significantly attenuated brain injuries and rescued the transcriptomic changes induced by acute HH. Functional investigation of the affected genes suggested that these two neuroprotectants protect the brain by targeting different pathways, with PROG enhancing erythropoiesis and TRIOL suppressing glutamate-induced excitotoxicity. Thus, this study advances our understanding of the pathology induced by acute HH and provides potential compounds for the development of neuroprotectant drugs for therapeutic treatment.

INTRODUCTION

High altitude, which represents one of the most extreme environments on Earth, creates hypobaric hypoxia (HH) conditions to which only a small proportion of people can adapt. As the brain is the most hypoxia-intolerant organ, most people who ascend too rapidly to high altitudes cannot acclimatize to the HH environment and frequently suffer from a range of symptoms caused by acute hypoxia (Wilson et al., 2009). In severe cases, patients can develop serious high-altitude cerebral edema (HACE). With over 35 million people around the world traveling to high-altitude regions annually for tourism, sport, or work (Martin & Windsor, 2008), brain damage caused by HH has become a health hazard for lowland people who move to high altitudes. So far, the main mechanisms underlying brain damage include metabolic disturbance of neural cells, increased permeability of brain microvasculature, and oxidative stress (Basnyat & Murdoch, 2003; Wilson et al., 2009). However, details on the cascades and gene regulatory networks through which brain damage is induced by HH have not yet been determined.

Several drugs, such as acetazolamide and dexamethasone, are proposed to treat or even prevent HH-induced brain damage and edema based on increased oxygen-carrying capacity or anti-inflammatory effects (Imray et al., 2010). However, their efficacy is still controversial and numerous side effects have been reported in relation to their clinical use (Imray et al., 2010; Rabinstein, 2006). Progesterone (PROG), a well-known endogenous neuroprotective sex hormone, is reported to be associated with acute HH acclimatization (Stein, 2008). In addition, serum levels of PROG are found to increase in men who travel to high-altitude areas (Basu et al., 1997), suggesting the possibility that PROG may ameliorate HH-induced brain damage. Furthermore, as first reported by our group (Chen et al., 2013a), 5α-androst-3β,5,6β-triol (TRIOL) is a novel neuroprotectant and analogue of the endogenous neuroprotective steroid cholestane-3β,5α,6β-triol (Hu et al., 2014), thus its effects on HH-induced brain damage are also worth investigation.

Previous pathological and molecular studies of HH-induced brain injuries have focused primarily on small mammals such as mice and rats (Imray et al., 2010). However, these rodents display distinct hypoxia tolerance and physiological responses to HH in comparison to humans. For example, both acute HH conditions and exhaustive exercise are necessary to induce HACE-like symptoms in rats (Guo et al., 2013), and very rapid ascent speed (velocity of 50 m/s within 5 min to 6 000 m) is required to induce HACE-like symptoms in mice (Huang et al., 2015). In addition, the potential drug targets in rodents may differ from those in humans due to the large genetic differences between rodents and primates after ~96 million years of evolutionary divergence (Nei et al., 2001). In particular, the central nervous system (CNS) of primates is much more complex than that of rodents (Lissa et al., 2013). Thus, given the limitations of rodent models, the establishment of non-human primate models for HH study is urgently needed and should provide valuable insights into the molecular cascades and gene regulatory networks underlying HH-induced brain diseases and further our understanding of the effects and molecular mechanisms of neuroprotectants as potential drugs.

In the current study, we elucidated the spatial and temporal influence of acute HH on gene expression and examined the possible effects of PROG and TRIOL on HH-induced brain damage in vivo. Specifically, we established a non-human primate model of HH-induced brain damage and profiled the transcriptomes of white blood cells (WBCs) and brain tissue from cynomolgus monkeys before and after acute HH exposure (i.e., increasing altitude). Extensive gene regulatory analyses revealed a dynamic change in the WBC transcriptome response to HH as well as the gene regulatory network of newly identified transcriptomic hub genes. Moreover, the application of the two neuroprotectants effectively protected the brains of cynomolgus monkeys from HH-induced injury. Finally, PROG and TRIOL exerted their effects via different pathways, the former through erythropoiesis and the latter by suppression of glutamate-induced excitotoxity.

MATERIALS AND METHODS

Experimental animals

All experiments were conducted in accordance with the Chinese Laws for the Protection of Animals. The experimental protocols were approved by the Ethics Committee of Zhongshan School of Medicine, Sun Yat-Sen University according to the ARRIVE (Animal Research: Reporting of In Vivo Experiments) guidelines (Kilkenny et al., 2010). All animal-based procedures were performed in strict adherence to the National Standards of Treating Experimental Animals (2006 version). Efforts were taken to minimize suffering and to ensure the welfare of monkeys during experimentation.

Twenty-four male cynomolgus monkeys (6.0–6.5 years old, 6.8–7.5 kg) were obtained from Gaoyao Kangda Laboratory Animals Science & Technology Co., Ltd, Zhaoqing, Guangdong Province, China. The animals were transported to the Third Military Medical University, Chongqing, China. Monkeys were singly housed in cages in a controlled environment at a temperature of 25±1 °C, relative humidity of 60%, and circadian 12 h light/dark cycle and were fed routinely.

For experimentation, monkeys were randomly divided into four groups, each containing six individuals. The first group was maintained under normobaric normoxia (NN) conditions (at an altitude of 320 m) as a control. The other three groups were subjected to hypobaric hypoxia (HH), as detailed below. Groups 3 and 4 were treated with progesterone (HH+PROG) and 5α-androst-3β,5α,6β-triol (HH+TRIOL), respectively.

After 48 h at a simulated altitude of 7 500 m, the HH monkeys were sacrificed. The NN monkeys were also sacrificed after 3 d in the NN environment. All monkeys were euthanized by bloodletting from the carotid under anesthesia using a mixture of injectable ketamine hydrochloride (0.06 mg/kg) and xylazine hydrochloride (0.02 mg/kg).

Pharmacological treatments

Intravenous (10 mg/mL) and extended-release intramuscular injections (50 mg/mL) of TRIOL were provided by Guangzhou Cellprotek Pharmaceutical Company Co., Ltd. (China). The TRIOL (Batch No. 101124) was dissolved in vehicle containing 0.9% sodium chloride and 20% hydroxypropyl-β-cyclodextrin (HP-β-CD) for intravenous injection, and in vehicle containing 12% glycerin, 20% HP-β-CD, and 0.19% CMC-Na for extended-release intramuscular injection. Injectable progesterone (20 mg/mL) was purchased from Zhejiang Xianju Pharmaceutical Co., Ltd. (China), with a CFDA ratification No. of Guo Yao Zhun Zi-H33020828. The HP-β-CD (Batch No. 100309) was purchased from Xi'an Deli Biology & Chemical Industry Co., Ltd. (China). Injectable ketamine hydrochloride (150 mg/mL) was purchased from Shenyang Veterinary Drugs Co., Ltd. (China) with an approval No. of (2011)060022668 by the Ministry of Agriculture of the People's Republic of China (PRC). Injectable xylazine hydrochloride (100 mg/mL) was purchased from Dunhua Shengda Veterinary Drugs Co., Ltd. (China) with an approval No. of (2010)070031581 by the Ministry of Agriculture of the PRC.

Acute hypobaric hypoxia and neuroprotective steroid agent treatments

For HH treatment, the simulated altitude in the hypobaric chamber was increased from 320 m to 6 000 m at a velocity of 3 m/s, and from 6 000 m to 7 500 m at a velocity of 2 m/s. Altitude simulation was started at 320 m (stage A), then suspended at 3 000 m (stage B), 4 500 m (stage C), and 6 000 m (stage D) for 50 min each, and finally maintained at 7 500 m for an initial 24 h (stage E) and then another 24 h (stage F) for a total of 48 h. At every stage, monkeys were maintained under HH conditions for 30 min at the corresponding altitudes, and were then subjected to blood sample collection and drug administration within 20 min. Because of safety concerns for the experimenters, the altitude was decreased from 7 500 m to 6 000 m for blood collection and drug administration at stage E within 30 min. Throughout the HH experimental exposure, the temperature, relative humidity, and air-flow velocity of the chamber were maintained at 22 °C, 60%, and 5 L/min, respectively.

For the HH+PROG group, animals were intramuscularly injected with 15 mg/kg of PROG three times, first at 12 h before HH treatment, then during stages D and E. For the HH+ TRIOL group, animals were intravenously injected with 10 mg/kg of TRIOL 5 min before HH treatment, and at stages B, C, D, and E. Owing to TRIOL's short half-life (t1/2) of 0.5 h, the extended-release intramuscular injection of TRIOL (30 mg/kg) was also administered at stages D and E.

Assessment of skeletal muscle coordination behavior

Behavioral assessment was performed by trained experimenters blind to group information. Monkey behavior was monitored 2 h after arrival at the simulated altitude of 7 500 m. We observed their activities over the next 30 min and employed a modified scale method to assess skeletal muscle coordination, with higher scores indicating more severely impaired behavioral states. The scoring system for monkey behavior assessment is presented in Supplementary Table S1.

Measurement of brain water content

Left brain hemispheres were used to measure brain water content with a precision electronic scale (Sartorius, BSA224S, Germany). Briefly, whole left hemispheres were immediately removed and measured as 'wet weight'. The left hemispheres were then placed in an electric thermostatic oven at 60 °C and weighed repeatedly every day until reaching a constant value (dry weight). The percentage of brain water content was calculated as: water content (%)=(wet weight-dry weight)/wet weight×100%.

Transmission electron microscopy (TEM)

The frontal cortices of brains were divided into 1 mm×1 mm×1 mm pieces and fixed in 2.5% glutaraldehyde overnight at 4 °C. The fixed brain tissues were dehydrated through an ethanol series, embedded in epoxy resin, and cut into 60 nm ultrathin sections with a Leica EM UC6 Ultramicrocut (Germany). The sections were mounted on copper grids and then stained in uranyl acetate and citric acid. Images were captured using a transmission electron microscope (Tecnai G2 Spirit Twin, FEI, USA). The experimental procedures for TEM were completed in a double-blind manner.

Nissl staining

The frontal cortices of brains were cut into 5 mm-thick pieces immediately after removal and then fixed in 4% paraformaldehyde. All samples were dehydrated in an increasing series of ethanol, cleared in xylene, and embedded in paraffin. Brain sections of 5 μm-thickness were cut using a rotary microtome (Shandon Finesse 325, Thermofisher Scientific, USA) and then deparaffinized with xylene and hydrated, followed by Nissl staining (0.5% cresyl violet) for histopathological assessment of neuronal injury. For Nissl staining, injured neurons were characterized by dark staining, condensed nuclei, shrunken cell bodies, or weak staining with irregular shapes. Images were acquired with a bright-field microscope (Nikon ECLIPSE Ti-U, Japan).

WBC isolation

WBCs were isolated from whole blood samples using Red Cell Lysis Buffer (Tiangen, RT122, China). In brief, red blood cells (RBCs) were lysed by gently mixing the blood sample with Red Cell Lysis Buffer (1:3), followed by incubation on ice for 5 min. The precipitate was then harvested by centrifugation at 2 500 g for 5 min at 4 °C. The WBC samples were ready for total RNA extraction after the isolation procedure was repeated three times.

RNA extraction, library construction, and sequencing

Total RNA was extracted and purified from WBCs and frontal cortices using Trizol reagent (Invitrogen, 15596, USA) according to the manufacturer's standard protocols. RNA quality and concentration were assessed using an Agilent Bioanalyzer 2100 and RNA 6000 Nano Kit (Agilent Technologies, USA). One of the TRIOL-treated monkeys was removed from RNA-seq analysis because the blood and brain RNA samples were poorly preserved due to depletion of dry ice during delivery to the sequencing company. One of the stage-A WBC samples from the HH+PROG group was also removed due to severe degradation of RNA.

For library construction, oligo-dT-coupled beads were first used to enrich poly-A+RNA molecules. First-strand cDNA synthesis was performed using random hexamers and Superscript II reverse transcriptase (Invitrogen, USA). Second-strand cDNA synthesis was performed using E. coli DNA PolI (Invitrogen, USA). A Qiaquick PCR purification kit (Qiagen, Germantown, MD, Germany) was used to purify the double-strand cDNA. cDNA was sheared with a nebulizer (Invitrogen, USA) into 100–500 bp fragments. The fragments were ligated to a Lumina PE adapter oligo mix after end repair and the addition of a 3' dA overhang. The 310–350 bp fragments were then collected by gel purification. After 15 cycles of PCR amplification, the libraries were subjected to paired-end sequencing using an HiSeq2500 (Illumina, USA) platform. Finally, 357.58 Gb of transcriptome sequences were generated from the 91 samples for transcriptome analysis.

Measurement of intraneuronal calcium concentration ([Ca2+]i)

The intracellular [Ca2+]i of cultured primary rat cortical neurons was measured as per Minta et al. (1989). Briefly, after 8 d of in vitro culture on glass cover slips, primary cortical neurons were washed once with Hank's balanced saline solution (HBSS) loading buffer. We added 100 μL of HBSS loading buffer containing 5 μmol/L Fluo-3/AM (Molecular Probes, F1241, USA) dropwise onto the cover slips, which were then incubated with neurons for 45 min at room temperature before being washed twice with HBSS testing buffer to remove residual fluorescent dye. Neurons were then transferred to the perfusion chamber of an OctaFlow Perfusion System (ALA Scientific Instruments, Farmingdale, NY, USA), which was placed on the stage of a TCS SP2 laser scanning confocal microscope (Leica Microsystems, Mannheim, Germany). Using an excitation wavelength of 488 nm and emission wavelength of 526 nm, at least 10 randomly selected neurons on every slide were scanned continuously at 6 s intervals to determine real-time fluorescence intensity. Dynamic changes in the fluorescence intensity of every neuron were recorded automatically by the Leica confocal system. Changes in [Ca2+]i (Δ[Ca2+]i) were represented by changes in fluorescence intensity before and after drug administration and were expressed as: Δ[Ca2+]i=(F (peak fluorescence intensity after administration) – F0 (baseline fluorescence intensity before administration))/F0×100%.

Statistical analysis of behavioral and neurological experiments

Statistical analyses were performed using SPSS 19.0 and graphs were drawn using GraphPad Prism 6.0. Most data, including brain water content and number of injured neurons, were analyzed using Student's t-test or one-way analysis of variance (ANOVA) with Dunnett's post-hoc test for multiple comparisons. Behavioral impairment scores were analyzed using the Kruskal-Wallis test. All data are represented as means±SD unless otherwise noted. P-values of less than 0.05 were considered statistically significant.

Source of genome reference and gene annotation

The reference genome and gene annotation of the cynomolgus monkey used in this study were originally generated by Yan et al. (2011) and downloaded from GigaDB ( http://gigadb.org/dataset/100003). Specifically, the following files were downloaded:

Genome assembly: ftp://parrot.genomics.cn/gigadb/pub/10.5524/100001_101000/100003/CE.cns.all.fa.gz

Gene annotation (GFF): ftp://parrot.genomics.cn/gigadb/pub/10.5524/100001_101000/100003/CE.gff.gz

Gene annotation (CDS): ftp://parrot.genomics.cn/gigadb/pub/10.5524/100001_101000/100003/CE.cds.fa.gz

Gene annotation (PEP): ftp://parrot.genomics.cn/gigadb/pub/10.5524/100001_101000/100003/CE.pep.fa.gz

Transcriptome alignment and expression level measurement

The Hisat2 (v2.0.4) (Kim et al., 2015) package was employed to align transcriptome reads to the reference genome of cynomolgus monkey with default parameters. Gene expression was measured as RPKM (reads per kilobase of transcript per million mapped reads) using an in-house script (supplementary file "count_reads_num. pl"), and library sizes were normalized with the geometric method in DESeq2 (v1.10.1) (Love et al., 2014) (Supplementary Data). Genes expressed robustly (i.e., mean RPKM>1) in at least one sample group were used for principal component analysis (PCA).

Identification of differentially expressed genes (DEGs)

We applied a strict protocol to identify DEGs between groups of samples. Candidate DEGs were first identified using DESeq2 (v1.10.1) (Love et al., 2014), edgeR (v3.12.1) (Robinson et al., 2010), and Cuffdiff (v2.2.1) (Trapnell et al., 2013) separately. For DESeq2, the default method was used to normalize library sizes, the fitting type of dispersions to mean intensity was set as "parametric", and the method to test DEG significance was set as "nbinomWaldTest". For edgeR, the trimmed mean of M-values (TMM) normalization method was used to normalize library sizes, the quantile-adjusted conditional maximum likelihood (qCML) method was used to estimate dispersions, and the genewise exact test ("exactTest") was used to calculate the significance of DEGs. For Cuffdiff, the method for library size normalization was set as "--geometric-norm" and other parameters were set as default. P-values were adjusted to allow for multiple testing by the Benjamini-Hochberg False Discovery Rate (FDR) (Benjamini & Hochberg, 1995). A DEG was required to meet the following criteria: (1) at least two of the three DEG-detection software reported a FDR<0.05; (2) the lower quartile of expression level in the up-regulated group was greater than the upper quartile of expression level in the down-regulated group; and, (3) the mean RPKM was not less than 5 in the up-regulated group. Of note, for convenience, the FDR values shown for the genes mentioned in the Results section were all from edgeR.

To verify the robustness of our DEG-detection strategy and to calculate the false-positive rate, samples were randomly divided into pseudo groups. For WBC samples, three pseudo groups were generated, each consisting of two samples randomly chosen from stage A, stage D, and stage F, respectively. For brain samples, two pseudo groups were generated, each consisting of three samples randomly chosen from the HH and NN groups, respectively. The procedures for DEG detection described above were performed, and the numbers of DEGs between pseudo groups of the same tissue were considered as false positives. In contrast to the thousands of significant DEGs identified between real groups, no DEGs were identified between pseudo groups for WBCs and brain tissue, confirming the robustness of our strategy and reliability of our results.

Gene ontology annotation, Kyoto Encyclopedia of Genes and Genomes (KEGG) annotation, and enrichment analysis

Gene ontology (GO) terms (Ashburner et al., 2000) for the cynomolgus monkeys were annotated according to their orthologous relationships with rhesus macaques and humans (Supplementary Table S5). KEGG pathways (Kanehisa & Goto, 2000) were annotated employing online KAAS (KEGG Automatic Annotation Server) tools ( www.genome.jp/tools/kaas/) (Moriya et al., 2007) using all protein sequences of the cynomolgus monkeys as inputs with a BBH (bi-directional best hit) method (Supplementary Table S5).

Fisher's exact tests were employed to identify whether a list of genes (foreground genes) was enriched in specific GO terms or KEGG pathways, with comparisons of the number of foreground genes annotated to the specific GO/KEGG, number of foreground genes not annotated to the specific GO/KEGG, number of background genes (excluding foreground genes) annotated to the specific GO/KEGG, and number of background genes (excluding foreground genes) not annotated to the specific GO/KEGG. P-values were adjusted to allow for multiple testing by means of the Benjamini-Hochberg False Discovery Rate (Benjamini & Hochberg, 1995).

Co-expression module and hub gene identification

The WGCNA (v1.51) (Langfelder & Horvath, 2008) package in R (v3.2.3) was employed for co-expression network analyses using the percentage of transformed RPKM of all robustly expressed genes. Co-expression modules were identified using standard methods, with correlation between genes calculated by the Pearson method, and the distance for hierarchical clustering calculated by the Ward method (Ward Jr, 1963). Similar modules were merged by WGCNA with the parameter MEDissThres=0.25. The eigengene (vector presenting expression dynamics of a module between samples) and kME (correlation of the gene with corresponding module eigengene) in every module were also calculated by WGCNA. For every module, genes with kME>0.9 were considered hub genes.

Analysis of impact of drugs on transcriptome dynamics

For WBCs, the sequencing of the HH and drug-treated groups was not performed in the same batch. Hence, laboratory conditions, personnel differences, and reagent lots could have caused variations between the groups, i. e., batch effects (Leek et al., 2010). It should be noted that stage A samples from the PROG and TRIOL groups were collected prior to drug injection; therefore, the stage A DEGs between the HH and drug-treated groups were likely a reflection of batch effects rather than drug-specific effects, and were thus excluded from subsequent analyses on DEGs between HH and drug-treated groups in stages D and F.

For brains, we identified the DEGs between the NN, HH, and drug-treated groups in pairwise comparisons. Firstly, we focused on the impact of drug treatment on HH-responding genes, which were DEGs between the NN and HH groups. If an HH-responding gene was a DEG between the HH and drug-treated groups but not a DEG between the NN and drug-treated groups (i. e., the expression level of the gene in the drug-treated group was not significantly different from that in the NN group but was significantly different from that in the HH group), the gene was considered to be a drug strongly responding gene (SRG). If an HH-responding gene was a DEG between the NN and drug-treated groups but not a DEG between the HH and drug-treated groups (i.e., the expression level of the gene in the drug-treated group was not significantly different from that in the HH group but was significantly different from that in the NN group), the gene was considered to be a drug non-responding gene (NRG). The other HH-responding genes with expression levels that fell between the drug-treated, NN, and HH groups were considered drug weakly responding genes (WRGs). If a gene was not significantly differentially expressed by HH in contrast to NN (i.e., not a DEG between the HH and NN groups) but was significantly differentially regulated by a drug in contrast to both NN and HH conditions, it was considered to be a drug-induced DEG. Even if a gene was significantly up- or down-regulated by HH, if that gene was further significantly up- or down-regulated but not restored by the drug, it was also considered to be a drug-induced DEG, not an NRG.

RESULTS

Monkeys develop severe brain injuries under acute HH conditions that can be significantly alleviated by PROG and TRIOL

Twenty-four healthy, adult male cynomolgus monkeys were classified into four groups (n=6): i. e., normobaric normoxia (NN), HH, HH+PROG, and HH+TRIOL. To characterize the behavioral and pathological changes caused by acute HH and to assess the drug effects, animals from the HH, HH+PROG, and HH+TRIOL groups were placed in a hypobaric chamber, which mimicked ascending altitudes from 320 to 7 500 m at a velocity of 2–3 m/s (HH group). Altitude simulation for each monkey started at 320 m (stage A) and was then sequentially suspended at 3 000 m (stage B), 4 500 m (stage C), and 6 000 m (stage D) for 50 min each. Altitude was then maintained at 7 500 m for 24 h (stage E) and then a further 24 h (stage F) for a total of 48 h (Figure 1A). During HH treatment, monkeys in the HH+PROG and HH+TRIOL groups were given either PROG or TRIOL, respectively (Figure 1B, C; see Methods). In parallel, monkeys in the NN group were placed in a normobaric normoxia environment at an altitude of 320 m to serve as controls.

Figure 1

Acute HH-induced behavioral and cerebral impairments were attenuated by PROG and TRIOL treatments

A: Experimental procedure for HH treatment of cynomolgus monkeys (HH group, n=6). Vertical bars along y-axis indicate medical definitions of high, very high, and extreme altitudes. Horizontal bars represent duration spent at each altitude. Blood drop and brain icons indicate altitudes and

The impact of acute HH exposure on behavior was assessed 2 h after ascending to 7 500 m. In sharp contrast to the NN monkeys, which exhibited high activity levels with frequent walking, climbing, and eating, the HH monkeys lost their balance and lay prostrate with limited body movements, or even remained in lateral or dorsal recumbency during the entire behavioral assessment. All six monkeys in the HH group displayed anorexia, vomiting, motor deficits, and ataxia. In contrast, monkeys in the HH+PROG and HH+TRIOL groups were able to maintain their balance and remain seated, with frequent head and forelimb movements. Several PROG or TRIOL-treated animals exhibited normal feeding behavior or short-distance movements (a few steps), although vomiting still occurred in some monkeys (two PROG-treated monkeys and four TRIOL-treated monkeys). We then employed behavioral deficit scoring to obtain a quantitative measure of skeletal muscle coordination by recording the movements of each monkey over 30 min (Supplementary Table S1) (Kito et al., 2001). All NN monkeys scored 0 (i. e., walked normally), whereas the HH monkeys scored between 14 (i. e., could not stand or sit for long) and 18 (i. e., no movement). PROG and TRIOL-treated HH monkeys scored between 10 (i. e., able to stand spontaneously) and 14; this was significantly lower than the scores of untreated HH monkeys (ANOVA P<0.05,Figure 1D), indicating that PROG and TRIOL were able to partially restore motor coordination in HH monkeys.

Cerebral edema is normally manifested as the increase of the water content in the brain tissue. We found that brain water content was significantly higher in the HH group (76.71%±0.30%) than in the NN group (76.16%±0.26%) (ANOVA P<0.01,Figure 1E). The proportion of brain water content was significantly reduced in both the HH+PROG and HH+TRIOL groups to 76.12%±0.46% and 76.28%±0.23%, respectively, similar to levels in the NN group (Figure 1E). Additionally, the pericapillary space was markedly wider in untreated HH brains than in other groups (Figure 1F), suggesting that blood-brain barrier (BBB) disruption and vasogenic edema, which are considered to cause cerebral edema in patients with HACE (Bärtsch & Swenson, 2013), had been induced by HH and reversed by drug administration. Nissl staining of the frontal cortex revealed many injured neurons in the HH group only. These cells were characterized by dark staining, condensed nuclei, shrunken cell bodies, or weak staining with irregular shapes (Figure 1G). The proportion of injured neurons in the HH group was 21.1%, significantly higher than that in the NN group (1.4%) (ANOVA P<10–3). This ratio was restored to 7.5% in the HH+PROG group (ANOVA P<0.01) and 6.7% in the HH+TRIOL group (ANOVAP<0.01,Figure 1H; see Methods).

Taken together, these findings suggest that cynomolgus monkeys develop severe motor impairments and brain damage after acute exposure to HH environments, and this can be significantly mitigated by either PROG or TRIOL.

WBC dynamic transcriptome analysis reveals three major regulatory modules responding to acute HH

To systematically investigate the genetic regulatory networks that participate in the HH response, we sequenced and analyzed the transcriptome of WBCs obtained at each of the six HH stages, A–F (n=6, 36 WBC samples in total). Based on PCA of 12 198 robustly expressed genes, the degree of HH exposure (i. e., a function of altitude plus exposure time), rather than individual differences, appeared to be the main factor affecting the genome-wide WBC gene expression patterns (Figure 2A). We then identified DEGs based on comparisons between any two of the six HH stages (see Methods) and obtained a total of 5 174 DEGs across all comparisons (ranging from 0–3 864 for different comparisons). As stage A (simulated altitude of 320 m) represents the NN condition, the number of DEGs between this and other stages reflects the degree of gene expression change in response to varying levels of HH exposure. The number of DEGs increased incrementally from stages B to D, reaching its highest at stage D. During stages E and F, the number decreased (Figure 2B), implying that the increased gene expression changes associated with HH were reversed in later stages, probably because of organismal acclimatization to hypoxia after 24 h of exposure to HH (Berglund, 1992).

Figure 2

Transcriptomic dynamics of WBCs in response to acute HH.

A: Principal component analysis (PCA) using genes robustly expressed in at least one of six HH stages (A–F). Colors denote stages of HH at which WBC samples were collected, as described in Figure 1A. B: Numbers of differentially expressed genes (DEGs) between stage A and each of the other five stages. C: Expression patterns of DEGs in three HH-responding modules, shown as log2-transformed reads per kilobase of transcript per million mapped reads (RPKM) fold changes of DEGs in each stage compared to stage A. Red and light blue lines indicate expression patterns of individual DEGs up- and down-regulated by HH, respectively. Bold lines show mean expression levels of all up- or down-regulated DEGs, with error bars indicating SD. UP and DOWN indicate exact number of up- or down-regulated DEGs by HH in each module. D: Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enriched in three HH-responding modules. Only significantly enriched pathways with a false discovery rate (FDR)-adjusted P-value of <0.05 were plotted. Dot color denotes P-value after FDR correction, and dot size denotes ratio of DEGs versus all expressed genes in each pathway. Different background colors denote classification of pathways. E: Interaction of pathways and symptoms according to reported articles. Blue color denotes pathways of impulse module; red color denotes pathways of late-responding module. Solid arrow line denotes reported direct interaction; dashed arrow line denotes reported associated relationship. F: Local co-expression network revealing genes in vitamin D receptor (VDR) complex as hub genes. Dark blue circles denote DEGs in impulse module; red circles denote genes in VDR complex and HIF1A as hub genes.

Co-expression analysis was next used to identify co-expression modules associated with HH changes (see Methods). Eight modules were identified, ranging in size from 64 to 5 607 genes (Supplementary Figure S1A, B). Three modules were significantly enriched with DEGs (Fisher's exact test P<10–15) and showed strong correlation between the expression pattern and progression of HH (ANOVA P<10–4; Supplementary Table S2). Hence, these were considered the primary HH-responding modules (HH modules). The first HH module contained 3 114 genes, 2 123 (68%) of which were DEGs (Supplementary Table S2). This module (termed the 'impulse module') showed an 'impulse pattern', whereby the expression levels were monotonically increased or decreased in the early stages of HH (stages A–D), but recovered to a certain extent during the later stages (stages E–F) (Figure 2C). The second HH-responding module (termed the 'sustained module') contained 571 genes, 479 (84%) of which were DEGs (Supplementary Table S2). In this module, changes in expression levels of genes that responded to early-stage changes in altitude were sustained in the later stages of HH (Figure 2C). The third module (termed the 'late-responding module') contained 552 genes, 398 (72%) of which were DEGs (Supplementary Table S2). Genes in this module showed no differential expression until the later stages of HH (Figure 2C).

Impulse responses of hypoxia/HIF-1 mediated innate immunity and inflammatory processes to acute HH

Genes in the impulse module were significantly implicated in several classical pathways critical to hypoxia response and inflammation (Figure 2D–E, and Supplementary Figure S1C). As noted in previous small mammal and human cell line studies, we confirmed that the HIF-1 signaling pathway also played a regulatory role in adaptive responses promoting optimal oxygen utilization (Semenza, 2009) and in inflammatory responses (Palazon et al., 2014) to hypoxia in in vivo primate experiments. Other innate immune and inflammatory pathways such as the NF-κB, Toll-like receptor, and FoxO signaling pathways, were also activated, probably forming a positive feedback loop with HIF-1 (Han et al., 2016). The inflammatory factor IL1B and chemokine CXCL1, which constitute links with HIF-1 and NF-κB signaling (Hatfield et al., 2010; Jung et al., 2003), were up-regulated by HH in the impulse module and are known to contribute to vascular impairment and vasogenic edema (Stamatovic et al., 2006), as observed in the brains of the HH-treated group (Figure 1F). Collectively, HIF-1/hypoxia pathways mediating the inflammatory response may contribute to pathogenesis in monkeys and were expressed in the impulse module during HH progression (Figure 2E). In contrast, the sustained module included genes mostly involved in adaptive immune responses (Figure 2D and Supplementary Figure S1C). Specifically, most key T cell markers, T cell receptors, and their specific kinases were consistently down-regulated by HH, suggesting overall suppression of the adaptive immune system in the HH environment.

Emergence of RBC-associated inflammatory processes in late stages

The late-responding module mainly contained genes involved in erythropoiesis and platelet activation (Figure 2D and Supplementary Figure S1C). Several genes encoding key factors in erythropoiesis, including erythroid transcription factors (e.g., GATA1 and KLF1), erythropoietin receptor EPOR, alpha and beta-globin genes (e.g., HBA1, HBA2, and HBB), hemoglobin stabilizing protein AHSP, key synthases of heme (e.g., HMBS, ALAD, and UROS), tropomodulin-1 TMOD1, and E3 ubiquitin-protein ligase TRIM58, were significantly up-regulated during stages E and F (Supplementary Figure S1D). This suggests that the enhanced expression of erythropoiesis genes was induced as an adaptive response after 28–52 h (stages E–F, Figure 1A) of HH exposure, which can regulate erythropoiesis and improve oxygen supply, thereby attenuating the impact of hypoxia in later stages. Key genes of blood coagulation and inflammation were also significantly up-regulated in late stages, including GP1BA, F2R, and other important procoagulant/antifibrinolytic factors (e.g., VWF, SERPINE1, and F13A1) that are associated with impairments of RBCs and vascular inflammation (Bergmeier et al., 2006; Gemmati et al., 2016; Sparkenbaugh & Pawlinski, 2013). In addition, the expression of the haptoglobin gene HP, a major scavenger of free hemoglobin, was also significantly up-regulated by HH (Supplementary Figure S1D), reflecting a homeostatic response to the increased release of free hemoglobin from impaired RBCs.

Intriguingly, in the late-responding module, we also found genes encoding critical markers of neutrophil degranulation, including neutrophil elastase (ELANE), myeloperoxidase (MPO), and cathepsin G (CTSG), the expression levels of which were up-regulated by HH (Supplementary Figure S1D). These activated neutrophils are tightly connected with BBB disruption, vascular inflammation, and vasogenic edema in injured brains (Ikegame et al., 2010; Liu et al., 2018), which could contribute to the development of the HACE-like symptoms observed in monkeys.

Vitamin D receptor (VDR) signaling is a novel key regulator involved in HH response

We next searched for intra-modular hub genes (i. e., genes with the highest degree of connectivity with other genes in the same module) in the three HH-responding modules to identify core regulators of the regulatory networks (see Methods). A total of 326 hub genes were identified (Supplementary Table S3), including those encoding key transcription factors, kinases, and receptors from the HIF-1 and FoxO signaling pathways (e.g., STAT3, FOXO4, and PTEN in the impulse module), T cell receptor signaling (e.g., LCK and PTPN4 in the sustained module), and erythrocyte development and maintenance (e.g., TMOD1, TRIM58, and KLF1 in the late-responding module).

Interestingly, our hub gene analysis also identified the vitamin D receptor (VDR) signaling pathway as an important but previously unreported pathway involved in the HH response. Specifically, genes encoding two key transcription factors, VDR and retinoid X receptor alpha (RXRA), and two co-activators, nuclear receptor coactivator 1 (NCOA1) and tripartite motif containing 24 (TRIM24), were identified as hub genes in the impulse module and were up-regulated at an early stage (Figure 2F). VDR and RXRA encode proteins that together form a heterodimeric complex that serves as a nuclear receptor for calcitriol, the hormonally active metabolite of vitamin D (Kato, 2000). Previous studies have revealed that VDR can inhibit the transcriptional activity of NF-κB (Chen et al., 2013b) and HIF-1 (Chung et al., 2009). VDR may also be up-regulated by Toll-like receptors and participate in the regulation of both the innate and adaptive immune systems (Adams & Hewison, 2008). These observations raise the possibility that the VDR pathway may function as a hub regulatory network that balances HIF-1 signaling and its inflammatory responses. Consistent with this, the vitamin D-binding protein in humans is associated with HH adaption in high-altitude dwelling native peoples (Ahmad et al., 2013), whereas VDR in mouse endothelial cells protects brains from hypoxia/reoxygenation-induced BBB disruption (Won et al., 2015).

Diverse effects of PROG and TRIOL on HH-induced WBC transcriptome dynamics

To evaluate the effects of two neuroprotective steroid agents on the blood cell transcriptome, we collected WBC samples for RNA-seq analysis from the HH+PROG and HH+TRIOL groups at stages D and F (Figure 1B), because stage D exhibited the most dynamic response to HH (Figure 2B) and stage F was the final stage of our analysis. WBC samples were also collected from the same animals at stage A, before drug injection, to serve as controls (Figure 1B). PCA revealed that neither PROG nor TRIOL treatment of HH animals had a significant impact on their WBC transcriptome dynamics, and the degree of HH exposure remained the primary differentiating factor between samples (Figure 3A).

Figure 3

Regulation of WBC transcriptomic dynamics by PROG and TRIOL

A: Principal component analysis (PCA) using genes robustly expressed in at least one of nine sample groups (i.e., no-drug treated group, PROG-treated group, and TRIOL-treated group in stages A, D, and F, respectively). Colors denote HH stages at which WBC samples were collected during HH+drug experiments, as described in Figure 1B, whereas shapes denote treatments. B: Number of differentially expressed genes (DEGs) between drug-treated and untreated groups in each HH stage studied. Stage A samples from both PROG and TRIOL groups were collected prior to drug injection; thus, DEGs of A/PROG*-A and A/TRIOL*-A reflect random variations between two groups of experimental monkeys rather than drug-specific effects. C: Erythrocyte-associated DEGs up-regulated by PROG. Red color and arrows denote genes up-regulated by PROG. D: Expression patterns and functions of TRIOL-induced DEGs in three HH-responding modules, shown as log2-transformed reads per kilobase of transcript per million mapped reads (RPKM) fold changes of DEGs in each sample group compared to stage A samples without TRIOL treatment. Light red and blue lines indicate expression patterns of individual DEGs in TRIOL-treated and untreated groups, respectively. Bold lines show mean expression levels of all TRIOL up- or down-regulated DEGs in each group, with error bars representing SD. UP and DOWN indicate exact number of up- or down-regulated DEGs by TRIOL in each module. Text below each module outlines functions associated with TRIOL-induced DEGs, with up and down arrows denoting gene-associated functions up- or down-regulated by TRIOL, respectively.

To further quantify the potential effects of the two neuroprotective agents on specific genes and pathways, we identified agent-induced DEGs by comparing the gene expression between samples from agent-treated (HH+PROG or HH+TRIOL) and untreated HH groups at stages D and F. We identified 134 DEGs induced by PROG in stages D (70 DEGs) and F (72 DEGs) and found that these DEGs overlapped significantly with those from the late-responding module (Figure 3B and Supplementary Figure S2A, Fisher's exact test: P<10-12). In particular, key erythrocyte-associated genes were among the DEGs with the largest expression changes (top DEGs) in stage F and significantly up-regulated by PROG, including the alpha and beta-globin genes (HBA1: FC=3.03, FDR=0.000 5; HBA2: FC=2.99, FDR=0.000 8; HBB: FC=3.40, FDR=0.000 1), heme synthases (SLC25A39: FC=2.10, FDR=0.008 3; ALAS2: FC=3.43, FDR=0.000 3), and bisphosphoglycerate mutase (BPGM: FC=2.33, FDR=0.026 6) (Figure 3C). Thus, although PROG appeared to have a relatively subtle effect on the WBC transcriptomes, we speculate that this steroid may function specifically to enhance erythropoiesis and elevate oxygen transport in the blood, thereby helping to relieve the symptoms caused by HH exposure.

Compared with PROG, TRIOL had significantly broader effects on the WBC transcriptome (Figure 3B), inducing 1 676 DEGs during stages D (195 DEGs) and F (1 605 DEGs). Interestingly, these TRIOL-induced DEGs overlapped significantly with HH-induced DEGs in all three HH-responding modules (Supplementary Figure S2A). Functional enrichment analysis revealed that TRIOL significantly up-regulated HIF1/hypoxia and RBC-associated inflammation, and down-regulated T cell receptor signaling during stage F (Figure 3D and Supplementary Figure S2B). Thus, in contrast to PROG, which had a specific function in erythropoiesis, TRIOL treatment altered the expression of a broader variety of genes in WBCs in terms of HH-related functions.

Profound transcriptomic alterations underlie pathological effects of HH-induced brain injuries

To investigate the brain gene regulatory network underlying the acute HH response, and the effects of PROG and TRIOL, we performed RNA-seq on frontal cortex samples (brains) from all monkeys after HH treatment (Figure 1A, B). Based on PCA of 14 200 robustly expressed genes, the brain samples could be clearly separated into the three groups: NN, HH, and drug treatment. The HH and NN groups could be separated at both the PC2 and PC3 level, as was performed for the WBC samples, indicating that HH was the main factor affecting the global brain gene expression pattern for monkeys without drug treatment (Figure 4A). Moreover, both agent-treated groups clustered tightly together, but separately from the HH group at the PC2 level (34.8% of total variance) and from the NN group at the PC3 level (23.8% of total variance) (Figure 4A), suggesting marked effects of both neuroprotective steroid agents on the brain transcriptome.

Figure 4

Transcriptomic changes in frontal cortex in response to acute HH, PROG, and TRIOL

A: Principal component analysis (PCA) using genes robustly expressed in at least one of four sample groups (normobaric normoxia [NN], HH, HH+ PROG, and HH+TRIOL). Colors denote sample group. B: Expression changes of HH-induced differentially expressed genes (DEGs) after PROG or TRIOL treatment, shown as log2-transformed reads per kilobase of transcript per million mapped reads (RPKM) fold changes of DEGs in each sample group compared to NN group. HH-induced DEGs are categorized into strongly responsive genes (SRGs), weakly responsive genes (WRGs), and non-responsive genes (NRGs) according to their degree of expression level recovery after drug treatment. Red and light blue lines denote expression patterns of individual DEGs up- and down-regulated by HH, respectively. Bold lines show mean expression levels of all HH up- or down-regulated DEGs in each group, with error bars representing SD. Percentages of SRGs and WRGs versus all HH-induced DEGs are presented above each plot. C: Expression changes of drug DEGs after PROG or TRIOL treatment, shown as log2-transformed RPKM fold changes of DEGs in each sample group compared to NN group. Red and light blue lines denote expression patterns of individual DEGs up- or down-regulated by drugs, respectively. Bold lines show mean expression levels of all DEGs up- or down-regulated by drugs in each group, with error bars representing SD. D: Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enriched in HH-induced or drug-recovered DEGs. Only significantly enriched pathways with a false discovery rate (FDR)-adjusted P-value of <0.05 are plotted. In HH DEGs, SRGs, WRGs, and NRGs, red and blue dots denote pathways enriched in DEGs up- and down-regulated by HH, respectively. In drug DEGs (PROG DEGs and TRIOL DEGs), red and blue dots denote pathways enriched in DEGs up- and down-regulated by drugs, respectively. Color intensity of dot indicates level of significance, and size of dot denotes ratio of DEGs versus all expressed genes in each pathway. E: Postulated regulation of excitatory glutamate signaling by TRIOL. Rectangles represent proteins of TRIOL DEGs involved in glutamate signaling pathway. Blue colors denote down-regulation by HH. Gene symbols of DEGs associated with each protein in the pathway are listed below. F: Time-course showing intraneuronal calcium [Ca 2+]i in primary cultured rat cortical neurons stimulated by glutamate with different doses of TRIOL. Glutamate stimulation was added at the fiftieth second in the experiment.

We identified 2 992 DEGs between the NN and HH groups (see Methods), with 1 446 up-regulated and 1 546 down-regulated by HH (Supplementary Table S4), reflecting considerable regulatory and functional rewiring in brainexpressed genes in response to acute HH exposure. DEGs showing the largest expression fold changes (top 1% of all DEGs) comprised the hemoglobins (carrier of oxygen) HBA1 (FC=10.69, FDR=0.000 0) and HBA2 (FC=4.12, FDR=0.000 2), and several angiogenesis-related genes, including adrenomedullin ADM (FC=13.58, FDR=0.000 0), vascular endothelial growth factor VEGFA (FC=5.48, FDR=0.000 0), and apelin APLN (FC=3.29, FDR=0.000 0), which were all up-regulated in HH brains. It is worth noting that many tightly associated DEGs were also significantly up-regulated by HH, including VEGFB (F=1.58, FDR=0.000 0), receptor of vascular endothelial growth factor FLT1 (FC=2.62, FDR=0.000 0), and receptor of adrenomedullin RAMP2 (FC=1.45, FDR=0.006 2). This strongly suggests that an adaptive response to HH, by increasing local tissue oxygenation and increasing oxygen delivery, occurred in monkey brains after 52 h of HH exposure. However, VEGF expression can also promote BBB leakage (Schoch et al., 2002), as observed in all monkeys exposed to HH (Figure 1F). Of note, CARTPT, a gene that encodes neuropeptides controlling feeding (Lambert et al., 1998), locomotor activity (Kimmel et al., 2002), and stress response (Rogge et al., 2008), also emerged as a top DEG (FC=6.88, FDR=0.000 3) between NN and HH groups, and its up-regulation in HH brains likely played a role in inducing the symptoms of anorexia and motor deficits in HH monkeys.

As with WBC transcriptome analysis, functional enrichment analyses revealed that HIF-1 signaling was significantly enriched by DEGs up-regulated in HH brains (Figure 4D), including EPAS1 (FC=1.37, FDR=0.028 5) and EGLN1(FC=1.48, FDR=0.004 7), two genes strongly associated with high-altitude adaptation in highland populations (Simonson et al., 2010; Yi et al., 2010), and HIF3A (FC=4.36, FDR=0.000 0), an important regulator of HIF1A and EPAS1 involved in adaptive responses to hypoxia (Ravenna et al., 2015). Moreover, given that the top DEGs (e.g., ADM, VEGF, APLN, HBA1, and HBA2) are regulated by HIF-1 signaling (Chen et al., 2012; Eyries et al., 2008; Ramakrishnan et al., 2014; Saha et al., 2014), this suggests that HIF-1 signaling also plays a main regulatory role in adaptive responses to acute HH in brains. Additionally, genes involved in the complement and coagulation cascades, and in leukocyte transendothelial migration, were up-regulated in HH-exposed brains (Figure 4D and Supplementary Figure S3A) and WBCs (Figure 2D), implying a link between peripheral inflammation and neuroinflammation in the CNS.

As a consequence of the decrease in oxygen supply under HH conditions, several pathways involved in aerobic respiration were found to be down-regulated in HH brains, including oxidative phosphorylation, respiratory electron transport chain, tricarboxylic acid (TCA) cycle, and pyruvate metabolism. On the other hand, pathways involved in anaerobic respiration were up-regulated, including glycolysis and mTOR signaling (Figure 4D and Supplementary Figure S3A). These findings suggest that, regardless of the adaptive response by stimulating angiogenesis and increasing oxygen delivery, the monkey brains suffered an energy deficiency, which may have induced autophagy and neuronal injuries through mTOR signaling (Dunlop & Tee, 2014).

Protein misfolding and endoplasmic reticulum (ER) stress are involved in many types of brain injury and neurodegenerative diseases (Giffard et al., 2004; Soto & Estrada, 2008). Indeed, we found that many key genes involved in the protein folding process were down-regulated by HH (Supplementary Figure S3A), including those encoding members of the HSP40 family, prefoldins, and FK506 binding proteins (Supplementary Table S4). Members of the HSP40 family are crucial partners for HSP70 chaperones (Qiu et al., 2006), and are important for hypoxic tolerance (Jain et al., 2014) and injury resistance (Hasegawa et al., 2018) in the CNS. Thus, their down-regulation may exacerbate HH-induced brain injury.

Taken together, our results revealed profound neural gene expression alterations after exposure to acute HH, and rewiring of genetic pathways, which may contribute to an adaptive response to acute HH and underlie the pathological effects of the associated brain injuries.

Brain transcriptomic changes induced by HH are reversed by treatment with PROG and TRIOL

We then assessed how PROG and TRIOL affect the expression of HH-altered genes (see Methods). Of all HH-altered genes, the expression levels were restored to normal (defined by reference to the NN group) in 50.67% and 48.09% of cases by PROG and TRIOL treatment, respectively, i. e., 'strongly responsive genes' (SRGs) (Figure 4B). Furthermore, 31.52% and 29.45% of HH-altered genes were partially restored by PROG and TRIOL treatment, respectively, i. e., 'weakly responsive genes' (WRGs) (Figure 4B). In addition, the expression levels of 17.15% and 20.82% of HH-altered genes were unaffected by PROG and TRIOL treatment, respectively, i.e., 'non-responsive genes' (NRGs) (Figure 4B). These results indicate that the abnormal expression levels of most HH-induced DEGs were completely or partially restored by PROG (82.19%) or TRIOL (77.54%) treatment. It should be noted that the SRGs and WRGs associated with the two drug treatments largely overlapped with each other (Supplementary Figure S3B), suggesting that PROG and TRIOL have broadly similar effects on reversing HH-induced transcriptomic changes in the brain.

PROG and TRIOL SRGs were particularly enriched in components of the key pathways of aerobic respiration (Figure 4D and Supplementary Figure S3A), suggesting an improvement in aerobic energy supply in response to PROG and TRIOL treatment. Specifically, of the 44 genes involved in oxidative phosphorylation and down-regulated by HH, 25 (56.8%) were PROG SRGs and 39 (88.6%) were TRIOL SRGs (Fisher's exact test, P=0.002). This implies that TRIOL may be more effective at improving the aerobic energy supply than PROG. The abnormal HH-induced expression of genes with roles in protein folding was also restored by both PROG and TRIOL treatment (Figure 4D and Supplementary Figure S3A), which may relieve neuronal injuries induced by protein folding and ER stress. In addition, the up-regulation of the behavior-related gene CARTPT by HH was also attenuated by both drugs to a certain degree (Supplementary Table S4), consistent with the alleviation of anorexia and motor deficits in drug-treated monkeys.

Of note, the up-regulations of several top DEGs involved in angiogenesis by HH were not recovered by PROG or TRIOL, including ADM, VWF, and VEGFA, suggesting that angiogenesis was also active in drug-treated groups (Supplementary Table S4). Intriguingly, PROG but not TRIOL further up-regulated several key genes involved in angiogenesis and oxygen delivery, including VEGFA (HH-vs-PROG: FC=1.23, FDR=0.032 3), its receptor FLT1 (HH-vs-PROG: FC=1.23, FDR=0.032 3), and hemoglobins HBB (HH-vs-PROG: FC=2.33, FDR=0.002 4) and HBA2 (HH-vs-PROG: FC=1.97, FDR=0.000 0) (Supplementary Table S4). These results further emphasize the potential role of PROG in relieving HH-induced symptoms by elevating oxygen transport in both WBCs and brains.

Gene regulatory mechanisms by which TRIOL suppresses glutamate-mediated excitotoxicity

In the brain transcriptome, we detected many genes that were specifically rewired by either PROG or TRIOL, even though their expression levels were not altered by HH treatment (see Methods, Figure 4C). These included 380 and 743 genes whose expression levels were induced by PROG and TRIOL, respectively. Intriguingly, among them, we found that TRIOL, but not PROG, significantly induced the enrichment of DEGs in many key excitatory neuronal signaling pathways (Figure 4D and Supplementary Figure S3A), including the glutamatergic synapse, L-glutamate import, glutamate receptor, and calcium signaling pathways. This suggests that TRIOL may have a specific effect on the regulation of glutamate-induced excitotoxicity. Glutamate induces excitotoxicity by promoting the excessive influx of intracellular calcium [Ca2+] through ionotropic glutamate receptors, particularly through the N-methyl-d-aspartate (NMDA) receptor. This is the main cause of hypoxia-induced neuronal death in many hypoxia-induced brain injuries (Wroge et al., 2012). Many key genes involved in this process were down-regulated by TRIOL, including NMDA receptors GRIN2A (FC=1.79, FDR=0.000 0) and GRIN2B (FC=2.23, FDR=0.000 3), vesicular glutamate transporter (VGLUT) SLC17A7 (FC=1.63, FDR=0.000 0), excitatory amino acid transporter (EAAT) SLC1A2 (FC=1.60, FDR=0.018 4), group I metabotropic glutamate receptor GRM5 (FC=1.48, FDR=0.030 4), and its downstream adenylate cyclase family members ADCY1 (FC=1.52, FDR=0.029 1), ADCY3 (FC=1.39, FDR 0.009 8), ADCY5 (FC=1.55, FDR 0.003 8), and ADCY9 (FC=1.59, FDR=0.000 2) (Figure 4E and Supplementary Table S4).

Glutamate-induced intracellular calcium ([Ca2+]i) overload is part of the signature indicating excitotoxicity in neurons (Scannevin & Huganir, 2000). To further examine whether TRIOL can indeed reduce neuronal excitotoxicity by directly suppressing the intracellular accumulation of calcium, we measured glutamate-induced [Ca2+]i in primary rat cortical neurons with and without TRIOL treatment. We observed that stimulation with 100 μmol/L glutamate alone dramatically increased [Ca2+]i in neurons within 50 s, whereas supplementation with TRIOL decreased these glutamate-induced effects in a dose-dependent fashion (Figure 4F). Taken together, our results suggest that TRIOL may suppress excitotoxicity via the indirect transcriptional regulation of key glutamate signaling-related genes through direct modulation of membrane permeability for calcium.

DISCUSSION

In this study, we demonstrated that cynomolgus monkeys exposed to acute HH exhibited various symptoms, including anorexia, vomiting, motor deficits, and ataxia, as observed in human patients suffering acute high-altitude diseases and which are mostly attributed to pathological changes in the brain (Wilson et al., 2009). The pathophysiological diagnoses in the brain suggest that cynomolgus monkeys developed severe cerebral edema, similar to that found in human HACE patients (Hackett & Roach, 2004). Hence, these data suggest cynomolgus monkeys can be used as non-human primate models for studying the disorders of oxygen homeostasis. Moreover, the genetic regulatory networks shown in this study to be involved in the development of hypoxia-induced brain malfunctions provide valuable insights into the relationship between the mechanism of oxygen homeostasis regulation and related diseases. Our study confirmed the fundamental role of several classical pathways, such as HIF-1 signaling, in the organism level response to hypoxia. Moreover, we also revealed that several molecules, such as VDR, serve as central connectors of different acute HH response pathways; these have been largely ignored in previous studies using human cell lines or mice as experimental models.

The dynamic transcriptome of WBCs from non-human primates may help to improve our understanding of the progression of HACE. First of all, most genes were categorized as belonging to the impulse module, including those involved in HIF-1 signaling; the expression of these genes was stimulated during early HH stages but recovered in late stages, suggesting that hypoxia is indeed the primary cause of injury under acute HH conditions. The hub genes in the impulse module may account for the subsequent expression changes in genes in the late stages of HH and may be responsible for the fundamental pathophysiological changes involved in acclimatization or injury. For example, RBC-associated inflammation is a novel pathway in late stages and may interact with the progress of HACE. Our transcriptomic analyses combining WBCs and brains also highlighted a complex organismal response to acute HH, with active crosstalk between the acclimatization that maintains oxygen homeostasis and adverse pathological feedback caused by inflammation. Specifically, in addition to the activation of the HIF-1 signaling pathway, consistent with previous reports (Semenza, 2009; Wilson et al., 2009), we also found that HH may induce a hierarchy of transcriptional programs of both innate and adaptive immunity, while progressively promoting inflammation in both the blood and brain.

The development of treatments for brain injury associated with HH, such as HACE, remains a formidable challenge. Our pharmacogenomic studies reported here provide preliminary but nevertheless valuable insights into the feasibility of a generalized neuroprotection strategy. Further, they present possible drug candidates, PROG and TRIOL, which both alleviated acute HH-induced brain injuries in our primate model. PROG is reported to be neuroprotective by stimulating breathing, increasing oxygen-carrying capacity, and modulating the GABAa receptor to suppress excitotoxicity (Singh & Su, 2013), whereas our results suggested a role in stimulating erythropoiesis and increasing oxygen delivery. Medroxyprogesterone (MPA), a synthetic analogue of PROG without neuroprotectant properties (Singh & Su, 2013), is used to treat chronic mountain sickness by stimulating breathing and increasing blood oxygen levels (Wright et al., 2004). However, MPA has no significant prophylactic effect for acute mountain sickness (AMS) (Wright et al., 2004), which further emphasizes the importance of neuroprotection in treating HACE. It is intriguing that the two neuroprotectants, PROG and TRIOL, appear to have different modes of action and act through distinct functional mechanisms to rewire the expression of genes in the blood and brain. Our previous study reported that TRIOL can protect primary cortical neurons from hypoxia/reoxygenation injury by maintaining mitochondrial functions (Chen et al., 2013a). Here, we further revealed that the neuroprotective effect of TRIOL may be mediated by attenuation of glutamate-induced excitotoxicity.

In conclusion, our research extends the transcriptomic analysis of acute HH responses at an integrated spatial and temporal level in non-human primates. We presented a kinetic transcriptomic analysis of WBCs and the transcriptomic profile of brains after exposure to acute HH. Our data revealed thousands of genes regulated by acute HH and categorized to different co-expression modules depending on their spatial and temporal dimensions. Our study also revealed that neuroprotective steroids such as PROG and TRIOL have significant recovery effects on acute HH-induced behavioral deficits and brain damage, suggesting the potential of using novel neuroprotectants as therapeutic drugs for acute high-altitude diseases. However, we also acknowledge that there are ~25 million years of evolutionary divergence between cynomolgus monkeys and humans (Yan et al., 2011), and some routine physiological data such as arterial blood pressure, arterial oxygenation, and ventilation were not collected in this study due to the limitation of experimental conditions. Thus, we cannot rule out that the pathological and gene regulatory changes observed in the cynomolgus monkeys may differ from that in humans suffering HH-induced disorders. Future studies with relevant physiological measurements could help better understand how animals respond to severe hypoxia, and how neuroprotective steroids function to protect HH-induced brain damage.

DATA AVAILABILITY

The underlying sequence data used in this study were deposited with the National Center for Biotechnology Information (NCBI) under accession No. PRJNA431946 and in the CNSA (https://db.cngb.org/cnsa/) of CNGBdb with accession code CNP0000293. The DEGs are shown in Supplementary Tables S3 and S4.

SUPPLEMENTARY DATA

Supplementary data to this article can be found online.

COMPETING INTERESTS

The authors declare that they have no competing interests.

AUTHORS' CONTRIBUTIONS

G.M.Y., G.J.Z., and Y.Q.G. designed the overall project; P.Z., Q.Y.L., W.Y., X. Y.Z., B.Z.L., L.X.S., and B.T.Q. carried out the transcriptome analyses; W.Y., Y.X.G., J.S.C., B.Z.L., W.B.Z., Z.B.Y., G.X., F.Y.L., and Z.W.T replicated animal models using the hypobaric chamber; B.Z.L., J.S.C., W.B.Z., M.Y., and Y.Q.W. prepared the samples for transcriptomic and pharmacodynamic testing; Q. Y. L., P. Z., and C. X. G. managed library construction and sequencing; B.Z.L., W.Y., J.S.C., and Y.L. performed the ultrastructural and histological assays; S. Z. L., J. X. Z., and Y. J. H. provided the TRIOL and pharmaceutical preparations; W.Y., J.S.C., and F.Y.L. handled the shipment and management of animals; J. W. and H. M. Y. contributed reagents/materials/analysis tools; W.Y., P.Z., Q.Y.L., B.Z.L., and L.X.S. drafted the manuscript; G.J.Z. and D.N.C. revised the manuscript. All authors read and approved the final version of the manuscript.

ACKNOWLEDGEMENTS

We thank China National Genebank for providing computing resources.