Deciphering regulatory patterns in a mouse model of hyperoxia-induced acute lung injury

Construction of the mouse model

All experimental procedures were performed and approved by Institutional Animal Care and Use Committee of Zunyi Medical University (IACUC; No. zyfy-an-2023-0047), ensuring compliance with ethical guidelines for animal research. A total of 12 C57BL/6J wild-type mice were purchased from Changsha Tianqin Biotechnology Co., Ltd. (Changsha, China), and randomly assigned into four groups: (i) Control Group (n = 6), exposed to zoom air in a regular environment; (ii) HALI group (n = 6), where mice were exposed to a hyperoxic environment (5.0 L/min pure oxygen for 72 h), with oxygen levels maintained at ≥90% saturation within an air-tight chamber, as previously described (Liu et al., 2019). This setup was closely monitored to maintain environmental conditions, including a temperature range of 25–27 °C and relative humidity between 50–70%. The chamber’s atmosphere was balanced by using soda lime to absorb the carbon dioxide exhaled by the mice, keeping CO₂ levels below 0.5%. The exposure protocol was rigorously scheduled for 23.5 h per day, with a 0.5-hour opening daily for essential maintenance tasks such as providing food, water, and clean bedding, ensuring consistency and replicability of the experimental conditions. No mice showed any signs of adverse weight loss in the two groups, and therefore no mice were excluded from the study. All experimental animals were anesthetized using isoflurane inhalation and subsequently euthanized, ensuring a painless end. Of note, euthanasia of animals strictly followed a protocol approved by our IACUC and the American Veterinary Medical Association. Animals needed for H&E staining were euthanized as described below (Tissue Collection and H&E staining procedure). To the further extent possible, our manuscript followed the ARRIVE guidelines (Checklist S1) (Perciedu Sert et al., 2020).

Tissue collection and H&E staining procedure

On the day of sacrifice, mice exposed to 90% oxygen for 72 h in a hyperoxic chamber were euthanized using an intraperitoneal injection of pentobarbital (50 mg/kg). After the completion of euthanasia, meticulous disinfection and preparatory measures were undertaken to ascertain the sterility of both the surgical field and the instruments employed. A representative lung tissue sample was typically procured post-thoracotomy, with the chosen site for this specific study being the lower lobe of the right lung. Lung tissues from the mice exposed to hyperoxia were then harvested and processed into sections on glass slides. These slides were deparaffinized in xylene and rehydrated through a series of graded alcohols (100%, 90%, and 70%). Subsequently, the sections underwent staining with hematoxylin for 3 min, followed by eosin staining for 2 min, and were then rinsed in distilled water. Finally, the sections were dehydrated through a second series of graded alcohols (70%, 90%, and 100%), cleared in xylene, and evaluated for lung injury using a light microscope.

RNA isolation and quality control

Lung tissues were collected from the mice and processed accordingly. Total RNA was extracted using TRIzol reagent (Magen, Guangzhou, China) following the manufacturer’s instructions. RNA quantity and purity were evaluated spectrophotometrically using a Nanodrop ND-2000 (Thermo Fisher Scientific, Waltham, MA, USA) based on the A260/A280 ratio. RNA integrity was further assessed using an Agilent Bioanalyzer 4,150 (Agilent Technologies, Santa Clara, CA, USA), and only RNA samples with an RNA integrity number (RIN) ≥8 was used for library construction.

Library preparation for RNA-seq

Library preparation was performed following the Abclonal mRNA-seq Library Preparation Kit (Abclonal, Wuhan, China) protocol. For each sample, 1 µg of total RNA was used to enrich and purify polyadenylated mRNAs using oligo (dT) beads. The purified mRNAs were then reverse transcribed into cDNA using reverse Rnase H and DNA polymerase I. The cDNA libraries were validated for appropriate size distribution using an Agilent Bioanalyzer 4150 (Agilent Technologies). Finally, the libraries were subjected to 150 bp paired-end sequencing on an Illumina NovaSeq 6,000 platform.

Quality control and read mapping

The mouse reference genome and gene annotation files (GRCm39, release M32) were downloaded from the GENCODE database (Frankish et al., 2021). Adapter sequences, low-quality reads (≥20% of bases with Phred score <15), and reads containing >5% ambiguous bases were trimmed using Trim Galore (v0.6.4) (Krueger, 2015). The quality of the trimmed reads was assessed using FastQC (v0.12.0) (Andrews, 2010). The cleaned reads were then aligned to the reference genome using STAR (v2.7.10b) (Dobin et al., 2013) with default parameters. Read counts were derived from the alignments and quantified at the gene level. In total, 12 samples were included and are evenly separated into 2 groups, with 6 biological replicates on each side. The statistical power of this experimental design, calculated in RNASeqPower (Hart et al., 2013) is 1.0.

Data deposit and acquisition

The datasets generated during this investigation, are meticulously deposited in the Gene Expression Omnibus (GEO) (Clough & Barrett, 2016). Processed gene counts from this study that are directly downloadable from the GEO website (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE237260).

The raw sequencing data, initially captured as fastq files, are accessible through the SRAToolkit (Sayers, O’Sullivan & Karsch-Mizrachi, 2022), enabling the conversion to full-length, paired-end sequence reads (fastq format).

Identification of differentially expressed genes

Gene expression for each sample was consolidated using ensemble IDs, and differential expression analysis was executed with DESeq2 (v1.36.0) (Love, Huber & Anders, 2014). Significant differentially expressed genes were identified based on the following criteria: log₂ fold change >1.5 or <−1.5 and P-values <0.05. Visualization of differentially expressed genes was accomplished using a volcano plot generated with the ggplot2 package (v3.3.6) (Wickham, 2011) and a heatmap created with the ComplexHeatmap package (v2.10.0) (Gu, Eils & Schlesner, 2016) in R 4.1.0 (Ihaka & Gentleman, 1996).

GO, KEGG and WikiPathways enrichment analyses

Gene Ontology (GO) (Ashburner et al., 2000), KEGG pathway (Kanehisa & Goto, 2000), and WikiPathways (Slenter et al., 2018) enrichment analyses were conducted to examine the biological functions of differentially expressed genes (DEGs). These genes were categorized into two groups according to their regulatory directions (up- or down-regulated), based on log₂ fold change values. Enrichment analysis was performed and visualized using the ClusterProfiler package (version 4.2.2) (Wu et al., 2021) in R software (version 4.1.0) (Ihaka & Gentleman, 1996). Terms with P-values less than 0.05 were considered significantly enriched.

Construction of WGCNA network

Differentially expressed genes were clustered into distinct groups using the WGCNA package (v1.70.3) (Langfelder & Horvath, 2008) in R, according to their expression patterns. Modules exhibiting a module significance score <0.05 were identified as recurrence-associated modules.

Similarly, distinct gene clusters were subjected to GO and KEGG pathway enrichment analysis using the ClusterProfiler package (v4.2.2) (Wu et al., 2021) to identify various aspects of biological functions. Terms with a P-value <0.05 were considered significantly enriched.

Identification of differential alternative splicing events

Splicing is an essential process that determines various cellular fates and directs diverse biological pathways. In this study, we utilized mapped read data from each sample (in BAM format) and combined it with mouse annotation files (in GTF format) as input. We employed rMATS (v3.3.0) (Shen et al., 2014) to quantify the splicing events for each sample and identify differential alternative splicing events

A statistical summary of alternative splicing for each sample was conducted using the summary.py script from rMATS (v4.1.2) (Shen et al., 2014). Differential alternative splicing events with a P-value <0.05 were considered significant. If a gene exhibited more than one type of alternative splicing, the most significant event was selected for visualization in a volcano plot using ggplot2 (v3.3.6) (Wickham, 2011) in R 4.1.0 (Ihaka & Gentleman, 1996).

Exploring the lncRNA and mRNA targets of miRNAs

There are lots of tools and databases available for predicting or collecting miRNA-target gene interactions. MultiMir (Ru et al., 2014) is one of the most comprehensive databases, incorporating both validated miRNA-target interaction databases (miRecords (Xiao et al., 2009, miRTarBase (Huang et al., 2020), TarBase Sethupathy, Corda & Hatzigeorgiou, 2006)) and predicted miRNA-target interaction databases (DIANA-MicroT-CDS (Paraskevopoulou et al., 2013), EIMMo, MicroCosm, miRanda, miRDB, PicTar, PITA, and TargetScan McGeary et al., 2019). Moreover, it provides an R package to enable user-friendly downloading of these results. In this study, we utilized MultiMir (version 1.22.0) (Ru et al., 2014) to query DEGs and differentially expressed long non-coding RNAs targeting miRNAs, and subsequently saved the obtained results.

Construction of ceRNA network

According to the ceRNA hypothesis, differentially expressed mRNA, differentially expressed lncRNA, and their corresponding targeting miRNAs were collected using the multiMiR (v1.22.0) (Ru et al., 2014) package in R (v4.1.0) (Ihaka & Gentleman, 1996). The lncRNA-miRNA-mRNA ceRNA network was constructed based on these results and visualized using Cytoscape (v3.10) (Shannon et al., 2003).

Generation of mouse model of HALI

A model of HALI was developed utilizing twelve C57BL/6J mice (Figs. 1A, 1B). Six mice were randomly selected and exposed to hyperoxia (≥ 90% oxygen) within an environmentally regulated chamber (25 °C–27 °C, 50–70% humidity) for 72 h (details in Methods). The remaining six mice were housed in an identical environment with standard oxygen concentration, serving as the control group. Hematoxylin and eosin (H&E) staining was conducted to verify the presence of HALI-associated pathological features. As depicted in Fig. 1B, exposure to hyperoxia resulted in considerable lung injury, characterized by alveolar capillary barrier disruption, pulmonary edema, and infiltration of inflammatory cells. These observations were identical to the clinical observation of the HALI (Kallet & Matthay, 2013), thereby confirming the successful establishment of the HALI murine model.

Identification of differentially expressed genes

The RNA-seq data from six mice with HALI and six control mice were analyzed for a comprehensive identification of genes and isoforms associated with HALI. Gene expression patterns of each sample were quantified with STAR (Dobin et al., 2013), and differential analysis was conducted utilizing DESeq2 (Love, Huber & Anders, 2014). A total of 727 differentially expressed genes (DEGs) were identified (Table S1), with genes exhibiting a log₂ fold change >1.5 or <−1.5 and the P-value <0.05 considered as significantly DEGs. To validate our analysis parameters, RNAseqPower was employed with a depth of 413, alpha of 0.1, effect size of 727/56953, and a sample size of 6, yielding an RNAseqPower of 1. This confirms the suitability of our parameters for DEGs identification in this dataset. Among these DEGs, 208 were up-regulated in the HALI group (including 43 lncRNAs) and 519 were down-regulated (including 205 lncRNAs) (Figs. 2A, 2B). Notably, genes associated with HALI were also detected differentially in this study (Table S1), such as Spink5, Gpx2, Lif, and Ear1, which are related to immune processes (Ahmed et al., 2022; Bonfill-Teixidor et al., 2024; Chen et al., 2024; O’Reilly et al., 2012) and, like Bax, may lead to cell death through the apoptosis process (Budinger et al., 2011). Gclc and Ftl1 are associated with ferroptosis (Luo et al., 2022; Wang et al., 2021).

Determination of the biological functions of DEGs

To address the biological functions of the DEGs, we conducted enrichment analyses of Gene Ontology (GO) (Ashburner et al., 2000), Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa & Goto, 2000), and WikiPathways (Slenter et al., 2018) via the ClusterProfiler (Wu et al., 2021) package in R (Table S2). The up-regulated DEGs in HALI were predominantly enriched in immune response-related functions, such as cellular response to interleukin-1, ERK1 and ERK2 cascade, and monocyte chemotaxis, within the biological process (BP) category of GO (Fig. 2C). Conversely, down-regulated DEGs (including Actn2, Agt, Casq2, Dhrs7c, Drd4, Grin1, Hamp, and so on) in HALI were connected to muscle cell development and myofibril assembly (Fig. 2C).

The results of the KEGG pathway enrichment analysis demonstrated that up-regulated DEGs were associated with cytokine-cytokine receptor interactions, tumor necrosis factor (TNF) signaling pathway, and interleukin-17 (IL-17) signaling pathway (Fig. 2D). Conversely, down-regulated DEGs were implicated in biological processes such as cardiac muscle contraction, neuroactive ligand–receptor interaction, and dilated cardiomyopathy, among others (Fig. 2D).

WikiPathways (Slenter et al., 2018) is an open, continuously updated, and curated pathway database that encompasses 202 pathways for Mus musculus, covering more than 4,500 genes (version 20230610). In this study, we utilized this database to augment the pathway resources of the KEGG database. The up-regulated genes were primarily linked to the chemokine signaling pathway, P53 signaling, and brain-derived neurotrophic factor (BDNF) pathway. In contrast, down-regulated genes were correlated with calcium regulation in cardiac cells, iron homeostasis, and so forth (Fig. 2E).

Given that the DEGs participated in multiple biological functions, it was crucial to discern the genes exhibiting distinct patterns to identify the complex biological processes involving these DEGs.

Construction of WGCNA network

Weighted gene co-expression network analysis (WGCNA) is a hierarchical clustering approach employed to categorize genes into discrete expression patterns. In our study, we utilized WGCNA to identify five distinct clusters of DEGs, employing a soft thresholding power of 14 and applying module filtering with specific cut-off scores (Fig. 3A, Table S3). Among these clusters, the G1 module encompassed merely four genes, and due to its score falling below the significance threshold, it was filtered out with subsequent analyses. The G5 module demonstrated a positive association with the HALI cohort while the G3, G4, and G2 modules displayed negative correlations (Fig. 3A). The G3, G5, and G4 modules comprised approximately 200 genes each. The G2 module, on the other hand, contained a minor fraction of DEGs (56 genes) (Fig. 3B).

The genes from various modules were subjected to GO and KEGG enrichment analyses. The enrichment results demonstrated that genes in distinct groups were involved in diverse biological processes (Figs. 3C, 3D, Table S4). Genes within the G5 module were correlated with immune response (e.g., response to interleukin-1, ERK1 and ERK2 cascade), which resembled the up-regulated genes’ enriched terms in DEG analysis. Genes within the G4 module were related to muscle system processes, muscle cell development, myofibril assembly, and so forth. Meanwhile, genes in the G2 module were associated with G1 to G0 transition, cardiac cell fate commitment, and other relevant processes. Genes in the G3 module were associated with alcoholism, systemic lupus erythematosus.

Upon examination of the enrichment analysis, it appears that genes exhibiting down-regulation in response to HALI events can be grouped into two distinct patterns. Conversely, up-regulated genes display a strikingly similar pattern across instances.

Identification of differential alternative splicing

Alternative splicing (AS) expands the number of proteins and is a key post-transcriptional regulatory mechanism. Therefore, investigating AS events is crucial for understanding the development of HALI. To this end, we employed rMATS to detect five types of AS events: skipped exon (SE), alternative 3′  splicing sites (A3SS), alternative 5′  splicing sites (A5SS), mutually exclusive exons (MXE), and retained intron (RI) (Shen et al., 2014). AS events with a P-value <0.05 were considered as differentially alternative splicing events (Table S5). For specific genes exhibiting multiple splicing events, the AS events with the minimum P-value and higher inclusion level difference were selected as the typical differential AS events (Table S6). As depicted in Fig. 4A, approximately 10% of AS events were differential. SE was the predominant AS type, followed by MXE, RI, A3SS, and A5SS, which had the fewest differential AS events. Most genes with differential AS exhibited only one type of AS. Notably, only three genes—Gas5, Smox, and Snhg17—demonstrated all five types of differential AS events (Fig. 4B). Gas5, a long non-coding RNA (lncRNA), functions as a competitive endogenous RNA (ceRNA) that sequesters miR429, thereby inhibiting Dusp1 in HALI (Li & Liu, 2020). Smox expression has been reported to increase under hyperoxia-induced neuronal damage in the retina (Narayanan et al., 2014).

Regarding genes not significantly expressed but exhibited differential alternative splicing (AS), we juxtaposed our DEGs with differentially AS events (Table S7). 159 genes showed differentially expressed but not AS, 1,744 genes showed AS but not differentially expressed, and only five genes showed both significance in AS events and DE genes (Asns, Lrrc74b, Nr4a1, Slc26a10, Snhg15) (Figs. 4C, 4D). Among those five genes, Snhg15 has been implicated in the suppression of cell apoptosis, expression of inflammatory cytokines, and oxidative stress response by up-regulating miR-362-3p expression in IPS-induced vascular endothelial cell (Liu et al., 2022). Overexpression of Nr4a1 could reduce the damage of Dex on hypoxia reoxygenation-induced in mouse pulmonary vascular endothelial cells (Dong et al., 2020).

Alternative splicing events involved in immune response and ferroptosis

Immune response plays a crucial role in the development of HALI (Chen et al., 2022; Soundararajan et al., 2022). Therefore, analyzing alternative splicing events of immune-related genes is of particular interest. To this end, we surveyed a list of immune genes from the Gene Ontology database (immune system process, GO:0002376), including 2,929 genes.

We discovered that 60 genes exhibited differential expressions during the HALI onset (Fig. 5A, Table S8). By intersecting these results with AS results from rMATS (Shen et al., 2014), we identified 1,237 events with significant alternative splicing when HALI developed, covering 422 genes. Our findings revealed those immune-related genes underwent alternative splicing during HALI (Fig. 5B). Notably, we observed a significant up-regulation of Bax during HALI (Fig. 5A), concurrent with significant AS happened to Tmbim6. Given in previous studies that Tmbim6 acts as a suppressor of Bax (Seitaj et al., 2018), the abnormal AS of Tmbim6 might account for the up-regulation of Bax.

Nearly half of the AS events were classified as skipped exons (SE), while only a small portion of AS events were categorized as alternative 3′  splicing sites (A3SS) and alternative 5′  splicing sites (A5SS) (Fig. 5C). Most genes harbored only one type of AS event; however, some genes exhibited multiple AS types. For example, Adgrf5, a gene associated with airway inflammation and potentially regulating CCL2-mediated inflammation (Kubo et al., 2019), harbored four types of AS events (excluding RI), and another gene, Lilrb4b, featured four types of AS events (excluded MXE) (Fig. 5D).

In addition to immune response, ferroptosis is one of the key mechanisms in the occurrence of HALI (Yin et al., 2021). We retrieved a list of 33 ferroptosis-related genes from the KEGG database (KEGG ID: map04216). A substantial number of these genes exhibited up-regulation in the context of HALI (Fig. 6A, Table S9), albeit not statistically significant, corroborating previous findings. Intriguingly, several genes demonstrated significant alternative splicing despite the absence of significant differential expressions, such as Gss, Trp53, and Lpcat3 (Fig. 6B). This observation underscores the intricate and precise regulatory mechanisms at the transcriptional level.

Identification of miRNA-target genes and lncRNAs

A total of 248 differentially expressed lncRNAs were identified in HALI, including H19 (Table S1). This finding prompted us to investigate the potential functions of these lncRNAs. One known function of lncRNAs is their role as molecular sponges for miRNAs. By binding to miRNAs, lncRNAs can prevent miRNAs from degrading their target mRNAs. (Salmena et al., 2011). MiRNA target genes were obtained from the multiMiR (Ru et al., 2014) database and annotated as either predicted or experimentally validated (Table S10). Not surprisingly, the majority of lncRNAs harbored one single miRNA targeting a specific mRNA, while Gm10447 harbored 16 miRNAs targeting Cacna2d2, as well as an additional 12 miRNAs targeting another gene, Has2 (Fig. 7A).

A comprehensive lncRNA-miRNA-target gene (ceRNA) network was constructed through the integration of differential analysis and miRNA target investigation (Fig. 7B, Table S10). Examination of the ceRNA network disclosed potential biological functions of specific lncRNAs, such as H19, GM10382, GM10447, and GM10244. Furthermore, it unveiled potential regulatory interactions of HALI-associated genes, including Slc7a11, Sox9, Mfsd2a, Nr1d1, Ptgs2, and Agt. These genes were implicated in the modulation of cellular ketone metabolic processes (GO:0010565), which is pivotal in the pathogenesis of HALI. Collectively, the ceRNA network analysis elucidated the putative roles of several lncRNAs and genes in the development of HALI, providing valuable insights for future investigations.