All experiments involving animals were carried out in accordance with the Guide for the Care and Use of Laboratory Animals and with approval from the Committee on the Ethics of Animal Experiments of Chongqing Medical University (Permit number: SYXK-(YU) 2012-0001). All surgeries were performed under anesthesia and every effort was made to minimize pain.

Female BALB/c mice (6–8 weeks old) were obtained from the Experimental Animal Center of Chongqing Medical University (Chongqing, China) and maintained in a specific pathogen-free environment with a 12 h light-dark cycle at a humidity of 55% ± 5% and temperature of 25 ± 2°C. All mice were given free access to food and water.

The RSV-A2 strain (VR-1540) was obtained from the American Type Culture Collection (ATCC, Manassas, VA, USA). The strain was grown in HEp-2 cells and purified by density gradient (Gias et al., 2008). The virus titer was 1.0 × 10⁸ PFU/ml, as determined by serial dilution plaque assay (McKimm-Breschkin, 2004). The BALB/c mice were infected with 60 μl of RSV-A2 via intranasal inoculation under 2%–5% isoflurane. Based on our previous study of lung tissue from RSV-infected mice, where antiviral interferon-α/β expression peaked at 12 h after infection (Ye et al., 2017), we chose the 12-h time point for sampling. At 12 h post-RSV infection, mice were anesthetized with CO₂ and sacrificed by cervical dislocation. The lungs, brains, spleens, and colons (1–2 cm of colon near the ileocecal region) were immediately collected and frozen with liquid nitrogen for at least half an hour and then stored at -80°C.

Lung, spleen, intestine, and brain samples from the control and RSV-infected groups (three mice each) were ground in liquid nitrogen and randomly sampled for total RNA extraction with TRIzol reagent (Invitrogen, CA, USA) following the manufacturer’s procedures. Total RNA quantity and purity were analyzed using a Bioanalyzer 2100 and RNA 6000 Nano LabChip Kit (Agilent, CA, USA) with an RNA integrity number (RIN) of >7.0. Poly (A) mRNA was isolated from 1 μg of total RNA using poly-T oligo magnetic beads (Invitrogen, CA, USA). Following purification, the mRNA was fragmented into 300-nt long oligonucleotides using divalent cations under elevated temperature. The cleaved mRNA fragments were then reverse transcribed to create a final cDNA library in accordance with the protocols of the mRNA Seq Sample Preparation Kit (Illumina, San Diego, USA) using the dUTP method. Average insert size for the paired-end libraries was 300 bp ( ± 50 bp). Paired-end sequencing (2 × 150 bp) on the Illumina NovaSeq™ 6000 platform was performed at LC-BIO Bio-Tech Ltd. (Hangzhou, China) following the vendor’s recommended protocols.

Raw counts were transformed by vst conversion and differentially expressed genes (DEGs) in the mouse groups were determined using the Bioconductor package “DESeq2” (version 1.26.0) (Anders and Huber, 2010). Statistical significance was set to |log2 fold-change (logFC)| > 1 and raw p < 0.05 for the spleen, intestine, and brain, but to |logFC| > 2 and raw p < 0.05 for the lung. Significance, fold-change, and gene expression patterns of the DEGs were visualized with volcano plots and heat maps using the OmicStudio tool (https://www.omicstudio.cn/tool). To identify common DEGs (co-DEGs) among the four organs, the online Venn tool (https://www.omicstudio.cn/tool/6) was used to draw a Venn diagram, and co-DEGs were retained for further analysis.

The biological functions of DEGs was identified by Gene Ontology (GO) (2017) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al., 2008) pathway enrichment analyses using the R package “clusterProfiler” (Yu et al., 2012) (version 3.14.3). The GO and KEGG enrichment results were visualized using the OmicStudio tool (https://www.omicstudio.cn). Terms and pathways were considered significantly enriched at a cutoff of p < 0.05. For co-DEGs among the four organs, GO was used to investigate variations in biological functions and Metascape (http://metascape.org/gp/index.html) was used to investigate variations in pathways post-RSV infection.

To explore the potential function of key hub genes, samples were divided into high- and low-expression groups according to the median value of the key hub gene. The “clusterProfiler” package was used for GSEA based on gene expression profiles. Previously annotated gene sets “Mm.c2.cp.reactome.v7.1.entrez.gmt” and “Mm.c2.cp.kegg.v7.1.entrez.gmt” were chosen as the reference gene list. A default algorithm with 1 000 permutations was applied to calculate p-values and enrichment scores. The gene set size filters were minimum of 10, maximum of 1 000, p = 0.05, and pAdjustedMethod = “BH”. Cut-off criteria were set to | Normalized Enrichment Score (NES)| > 1 and p < 0.01.

The deconvolution algorithm CIBERSORT contains gene expression reference values from a signature matrix of 547 genes in 22 types of immune cells (Newman et al., 2015). To evaluate abundance of immune infiltrates in various organs, we uploaded the gene expression matrix data to the CIBERSORT (https://cibersort.stanford.edu/) web portal, and the algorithm was run using the LM22 signature for 100 permutations. The percentage of each immune cell type in the samples was calculated and displayed in a bar plot. The “vioplot” package (verison 0.3.7) (Hu, 2020) was then used to draw violin diagrams to visualize differences in immune cell infiltration. The Wilcoxon test was used to analyze differences in immune cell fractions between the RSV-infected mice and healthy controls.

A protein-protein interaction (PPI) network of DEGs was constructed using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) (Snel et al., 2000) database. An interaction score > 0.4 was regarded as statistically significant. The PPI network was then visualized using Cytoscape (version 3.8.2) (Shannon et al., 2003) and the cytoHubba plug-in was employed to calculate the degree of interaction between DEGs, with the top five genes defined as hub genes.

For gene validation, mRNA expression profiles in the lung tissues of RSV-infected mice were obtained from the Gene Expression Omnibus (GEO) database (accession number: GSE111861). GSE111861 contains three RSV-infected intact, three RSV-infected TLR4 mutant, and three wild-type mice. Preprocessed gene expression data and sample annotation were obtained from GEO (Series Matrix Files). Differential analysis was performed for GSE111861 using the “limma” package (version 3.42.2). A boxplot was generated to visualize the gene expression data using the “ggplot2” package (version 3.3.4) (Ginestet, 2011).

Non-infected and 12, 24, and 48 h post-RSV infected mice were sacrificed and whole lungs, intestines, brains, and spleens were removed. Total RNA was extracted from the tissues using TRIzol Reagent (Invitrogen, USA) and an RNA Extraction Kit (TR205, Tianmo Biotech, China). Total RNA concentration and purity were assessed using a NanoDrop 2000, and 1 μg of RNA was used for first-strand cDNA synthesis with a PrimeScript RT Kit (TaKaRa, Japan) according to the manufacturer’s instructions. The cDNA was amplified using a VeriQuest Fast SYBR Green qPCR Kit (Invitrogen, USA) with the following primers: Hpx (forward 5’TAGCTGGCCCATTGCTCATC3’; reverse 5’-GAGGGCTCCCAAGTTCCTTC-3’) and Gapdh (forward 5’-CATCACTGCCACCCAGAAGACTG-3’; reverse 5’-ATGCCAGTGAGCTTCCCGTTCAG-3’). The PCR cycle conditions were: 95°C for 2 min, then 40 cycles at 95°C for 5 s and 60 °C for 30 s. The fold-change was obtained using the 2^-ΔΔCt method with Gapdh calibration. RSV N gene were performed using primers and probes: forward primer, 5’-AGATCAACTTCTGTCATCCAGCAA -3’, reverse primer, 5’-TTCTGCACATCATAATTAGGAGTATCAAT-3’ and probe, 5’-FAM-CACCATCCAACGGAGCACAGGAGAT-BHQ1-3’. The PCR cycle conditions were: 50°C for 2 min, 95°C for 10 min, then 40 cycles at 95°C for 15 s and 60°C for 1 min. For absolute quantification, the RSV N gene copy was calculated using a plasmid DNA standard curve.

Total protein in tissues was extracted using a Total Protein Extraction Kit (KGP2100, KeyGEN, China). Western blotting was performed following previously reported procedures (Ye et al., 2017). The primary antibody was anti Hpx (ab133415, Abcam, USA). Gapdh (200306, ZENBIO, China) was used as a loading control. Image collection was performed using Quantity One v4.6.2 (USA) and densitometry analysis was conducted using ImageJ v2.0.0 (USA).

Pearson correlation analysis of key hub genes and immune cells was performed to analyze the immune mechanism during RSV infection using the “ggstatsplot” R package (Patil, 2021).

All quantitative data are presented as mean ± standard error of the mean (SEM). The proportion of innate immune cells between groups was compared using the unpaired Student’s t-test. Pearson correlation analysis was performed to reveal the relationship between key hub genes and immune cells. Differences in mRNA expression between/among groups were determined using Student’s t-test or one-way analysis of variance (ANOVA) respectively followed by Tukey’s multiple comparisons in GraphPad Prism v7.0.4 (GraphPad Software Inc., USA). A p-value of less than 0.05 was considered statistically significant.

The distribution of DEGs in the different organs is shown in the volcano plot in Figures 1A–D , with the 10 most significant DEGs marked. In total, we identified 2 417, 750, 701, and 1181 DEGs in the lung, intestine, brain, and spleen, respectively ( Figure 1E ). For more accurate functional enrichment, the DEG threshold for the lung was set to |logFC| > 2. Among the tissue types, the lungs harbored the greatest number of DEGs, indicating that the lungs exhibited the most significant mRNA changes in response to RSV. The significantly up- and down-regulated DEGs are displayed in the heat map ( Figure S1 ), which shows different distributions of gene expression patterns between the RSV-infected and control mice. Gene expression in the four organs post-RSV infection, after hierarchical clustering, is shown in the heatmap in Figure 1F , with genes in boxes indicating tissue-specific gene expression signatures. Clustering analysis revealed that the lung and spleen had the most similar gene expression profiles post-RSV infection ( Figure 1F ).

GO and KEGG enrichment analyses were used to evaluate the potential functions of the DEGs in the lung, intestine, brain, and spleen of RSV-infected mice. Figures 2A–H show the top 10 enriched GO and KEGG terms in the four organs. Based on GO analysis, the DEGs in the lung were mainly enriched in virus defense response, cytokine, and innate immune response ( Figure 2A ). For KEGG analysis, the lung DEGs were mainly enriched in the cytokine-cytokine receptor interaction, completement and coagulation cascades and retinol metabolism pathways ( Figure 2E ). In the intestine, DEGs were significantly enriched in heme binding, tetrapyrrole binding, steroid hydroxylase activity, and long-chain fatty acid metabolic process based on GO analysis ( Figure 2B ), and enriched in the steroid hormone biosynthesis, completement and coagulation, and IL-17 pathways based on KEGG analysis ( Figure 2F ). In the brain DEGs, enriched GO terms included neutrophil and leukocyte migration and cell and neutrophil chemotaxis ( Figure 2C ) and enriched KEGG pathways included the completement and coagulation cascade and IL-17 signaling pathways, as found in the lung, as well as the cytokine-cytokine receptor interaction and tumor necrosis factor (TNF) signaling pathways ( Figure 2G ). In the spleen DEGs, enriched GO terms included high density lipoprotein particle and steroid metabolic process ( Figure 2D ) and enriched KEGG pathways included the steroid hormone biosynthesis and completement and coagulation signaling pathways ( Figure 2H ). These results indicate that the different organs showed both divergent and common responses to RSV infection.

Due to the limitations of DEG enrichment analysis, key genes with moderate expression changes may be omitted by DEG criteria and important pathways underlying pathophysiological changes may be overlooked. Thus, we conducted GSEA to further investigate functional changes correlated with RSV infection in different organs. Based on GSEA, we identified 201, 42, 72, and 56 markedly up-regulated pathways ( Figure 2I ), as well as six, 53, 23, and 33 down-regulated pathways in the lung, intestine, brain, and spleen, respectively ( Figure 2I ). Up-regulated processes in the lung mainly involved immune- and virus recognition-related pathways, which were not significantly increased in the other three organs ( Figure 2J ). Moreover, the apoptosis and regulation of autophagy pathways were only up-regulated in the lung. Interestingly, the lung and brain shared several elevated pathways, including interferon- and STING-mediated induction of host immune responses, thus suggesting a relationship between the lung and brain in RSV infection.

Although all the four organs had immune response changes post RSV infection, we wondered if the immune responses observed in these organs were secondary to lung as part of a systemic response or were the intestine, brain and spleen directly infected by RSV. Therefore, we examined RSV N gene copies in the four organs, and found that RSV genome was detected only in lung ( Figure S2 ), indicating RSV can’t directly infect intestine, brain and spleen. Together, these findings suggest that the lungs are a major virus-host battlefield during RSV infection and the intestines, brains and spleens could be influenced secondary to lung.

To further assess immune response changes in the four organs post RSV infection, we conducted immune cell infiltration with CIBERSORT. The percentages of the 22 types of immune cells in the four organs are shown in the violin plots ( Figure 3 ). Compared with the control group, plasma cells, resting natural killer (NK) cells, M0 macrophages, and resting mast cells were decreased in the lungs of RSV-infected mice, whereas T follicular helper (TFH) cells, activated NK cells, M1 macrophages, activated mast cells, and eosinophils were increased ( Figure 3A ). Several other studies have also reported increased TFH cells, activated NK cells, and eosinophils in the lung during RSV infection (Harrison et al., 1999; Li et al., 2012; Pyle et al., 2021). These findings demonstrate the reliability of CIBERSORT in predicting immune cell infiltration. Interestingly, increased M1 macrophages and activated mast cells and decreased resting mast cells have not been reported previously in the lung post-RSV infection. Our results also showed that immune cells were not significantly altered in the intestine or brain ( Figures 3B, C ), but monocytes were decreased and neutrophils were increased in the spleen, respectively ( Figure 3D ). These changes in neutrophils and monocytes in the spleen post-RSV infection have not been reported previously. Our study demonstrated the occurrence of distinct immune cell populations among different organs during infection; however, the specific mechanism and function of some types remain to be elucidated.

To elucidate the common mechanism underlying multi-organ pathology after RSV infection, co-DEGs were obtained for further analysis. As shown in Figure 4A , 16 genes overlapped in the four organs. Subsequently, the PPI networks of these co-DEGs were constructed and 5 hub genes were screened, i.e., Hpx, Apoh, Uox, Cyp2d26 and Serpinf2, with Hpx identified as the most key hub gene by cytoHubba ( Figure 4B ). Expression of the five hub genes in the different organs was shown in Figure 4C and GO analysis of the hub genes was shown in Figure 4D . According to Metascape pathway enrichment analysis, the 16 co-DEGs were enriched in the platelet degranulation and neutrophil degranulation pathways ( Figure 4E ). These results suggest that aberrantly expressed Hpx and activation of neutrophil degranulation and neutrophil degranulation pathways may play crucial role in the RSV induced multi-organ pathology. As Hpx was identified as the most key hub gene in the PPI network, we chose it for further analysis.

We constructed RSV-infection mouse models at 12, 24, and 48 h post-RSV infection to confirm the RNA-Seq results. Western blotting showed that the protein expression level of Hpx was elevated in the four organs post-RSV infection up to 48 h ( Figures 5A–D ), as confirmed by qRT-PCR ( Figure S3 ). To further validate our results, we downloaded the GSE111861 dataset from the GEO database and analyzed the differential expression of Hpx mRNA post-RSV infection in the lung. By comparing the expression levels of Hpx in four groups of the GSE111861 dataset (i.e., TLR4 intact, TLR4 mutant, TLR4 intact post-RSV infection, and TLR4 mutant post-RSV infection), we found that Hpx expression was elevated in TLR4 intact mice at day 1 post-RSV infection, similar to our data ( Figure 5E ). Interestingly, Hpx expression was not elevated in TLR4 mutant mice post-RSV infection ( Figure 5F ). Therefore, Hpx expression in TLR4 intact mice post-RSV infection was significantly increased compared to the TLR4 mutant groups post-RSV infection ( Figure 5G ). Notably, in the research based on GSE111861 dataset, RSV-infected mice with lower Hpx expression (TLR4 mutant groups post-RSV infection) had lower virus titer and lighter lung inflammation, indicating that Hpx expression may associate with severity of RSV infection.

To assess if Hpx expression in the four organs correlate with each other, we visualized the Hpx expression changes post RSV infection in the four organs in one histogram and then showed the correlation of Hpx expression among the four organs with a correlation heat map of fragments per kilobase million (FPKM) values. Results showed the Hpx was significantly increased in all the four organs with the lung had the maximum fold change of Hpx mRNA expression ( Figure 6A ). Hpx mRNA expression levels in lung, intestine, and spleen were significantly positively correlated with each other and correlations of Hpx expression between intestine, spleen and brain were much lower, though it was high between lung and brain ( Figure 6B ).

To evaluate the relationship between Hpx expression and immune cell infiltration status, we further used Pearson correlation analysis to determine the correlations between Hpx expression and immune cell fractions. Correlation analysis between Hpx and immune cell fractions in the lung showed that Hpx was strongly positively correlated with increased TFH cells, activated mast cells, M1 macrophages, activated NK cells, and eosinophils, and negatively correlated with decreased plasma cells, resting NK cells, and resting mast cells ( Figure 6C ). In the spleen, Hpx showed positive association with neutrophils and negative association trends with monocytes ( Figure S4 ).

Finally, to explore the potential function of Hpx post-RSV infection, we conducted GSEA to analyze enriched KEGG and Reactome pathways based on Hpx expression in the four organs. Gene sets associated with the completement and coagulation cascades, JAK-STAT, interferon-α/β, platelet activation, and aggregation signaling pathways were highly up-regulated in the groups with high Hpx expression ( Figures 6D–G ), indicating that Hpx may act together with these pathways. Together, these results suggest that Hpx may be vital in the multi-organ pathogenesis of RSV infection for its potential participation in immune cell infiltration and pathways.