Evaluation of weight changes in the experimental groups
We observed lethargy, ruffled fur, and loss of appetite, but no death in any IAV-challenged mice, the degree of weight loss of these two groups with 8 mice per group were shown in Fig. 1 (P = 0.006). The viral load and pathological damage of lung tissue 6 day post infection are depicted in Supplementary Figs. 1 and 2, respectively. Anti-PD-1 antibody treatment significantly alleviated the severity of necrosis and inflammation based on the gross and microscopic lesions and decreased the viral titers in the lungs.
Differential expression analysis
In the differential expression analysis of lungs (Fig. 2), 85 differentially expressed mRNAs (dif-mRNAs) were obtained of which 76 were upregulated and 9 were downregulated, functional genes in PD-1 antagonist treatment followed by A/PR8(H1N1) infection group included Nppa, Myl7, Car3, Itprip, CAMP (log2FoldChange of 11.5, 2.5, 5.7, 2.5 and 3.4 respectively) et al.; 36 differentially expressed miRNAs (dif-miRNAs) were identified, of which 19 were upregulated and 17 were downregulated; 90 differentially expressed lncRNAs (dif-lncRNAs) were obtained, including 70 upregulated and 20 downregulated; and 22 differentially expressed circRNAs (dif-circRNAs) were found, of which 13 were upregulated and 9 were downregulated.
In the spleen data, 45 dif-mRNAs were obtained, of which 18 were upregulated and 27 were downregulated, functional genes in PD-1 antagonist treatment followed by A/PR8(H1N1) infection group included CAMP, Ltf, Ly6g, Bola1 (log2FoldChange of 2.1, 2, 2.6, and 2.5 respectively) et al.; 36 dif-miRNAs were identified, of which 19 were upregulated and 17 were downregulated; 57 dif-lncRNAs were obtained, including 22 upregulated and 35 downregulated; and 24 dif-circRNAs were found, of which 18 were upregulated and 6 were downregulated.
Functional enrichment analysis of dif-mRNAs and dif-miRNAs in lungs and spleens
KEGG and GO analyses were used to investigate the functional associations of gene expression changes. Targeted genes of dif-mRNA and dif-miRNAs of lungs and spleens of the two groups: PD-1 antagonist followed with A/PR8(H1N1) infection group vs. Isotype control followed with A/PR8(H1N1) infection were predicted (Figs. 3 and 4). The gene lists used in the dif-mRNAs analysis contained 18,455 and 17,818 genes for lungs and spleens, respectively. 1290 and 1290 genes were analyzed for lungs and spleens for dif-miRNAs. For GO, biological process, cellular component, and molecular function were selected as the annotation categories for clustering. Once the tool identified enriched ontologies for a particular gene list, it clusters those that have a statistically significant overlap in terms of their constituent genes. P-value was set < 0.05, the dif-mRNAs were enriched in 11 pathways in lungs and 6 pathways in spleens. Dif-miRNAs were enriched in 11 pathways in lungs and 26 pathways in spleens. There was little degree of overlap of dif-mRNAs and dif-miRNAs in lungs between the most enriched clusters. The most enriched clusters of dif-mRNAs of lungs were related to muscle and heart biological behavior. More than 85% of the dif-miRNAs enriched clusters in lungs and spleens overlapped with each other, including localization, metabolic process, positive regulation of metabolic process, and regulation of molecular function in the biological process category; intracellular part, cytoplasm, intracellular, and membrane-bounded organelle in the cellular component category; and protein binding, enzyme binding, and molecular function regulator in the molecular function category.
Functional Enrichment analysis of mRNAs and miRNAs in lungs and spleens obtained from IAV infection mice treated with anti-PD-1 antibody clearly highlighted myocardial damage related to viral infection, mitogen-associated protein kinase (MAPK) signaling pathways, RAP1 (Ras-related protein 1) signaling pathway, and Axon guidance.
The top 20 pathways and GO terms (BP (Biological Process), CC (cellular component), and MF (molecular function)) enriched by 45 dif-mRNAs and 36 dif-miRNAs of spleens of the following groups: PD-1 antagonist treatment followed by A/PR8(H1N1) infection group vs. isotype control followed by A/PR8(H1N1) infection group. (A) Top 20 pathways enriched by dif-mRNAs (B) Top 20 pathways enriched by dif-miRNAs (C) Top 20 GO terms enriched by dif-mRNAs (D) Top 20 GO terms enriched by dif-miRNAs.
Enrichment analysis of lncRNA and circRNA-related target genes
KEGG and GO analysis was performed for dif-lncRNA and dif-circRNA-related target genes (Figs. 5 and 6). P-value was set < 0.05, the dif- lncRNA target genes were enriched in 9 pathways in lungs and 12 pathways in spleens. The dif- circRNA target genes were enriched in 7 pathways in lungs and 4 pathways in spleens. There was a little degree of overlap of lncRNAs and circRNAs in lungs and spleens between the most enriched clusters except for Hypertrophic cardiomyopathy, MAPK signaling pathway, and the AMP-activated protein kinase (AMPK) signaling pathway.
Top 20 pathways and GO terms (BP (Biological Process), CC (cellular component), and MF (molecular function)) enriched by 90 dif- lncRNAs and 22 dif-circRNAs of the lungs (A) Top 20 pathways enriched by dif-lncRNAs. (B) Top 20 pathways enriched by dif-circRNAs (C) Top 20 GO terms enriched by dif-lncRNAs (D) Top 20 GO terms enriched by dif-circRNAs.
Top 20 pathways and GO terms (BP (Biological Process), CC (cellular component), and MF (molecular function)) enriched by 57 dif- lncRNAs and 24 dif-circRNAs of spleens (A) Top 20 pathways enriched by dif-lncRNAs. (B) Top 20 pathways enriched by dif-circRNAs. (C) Top 20 GO terms enriched by dif-lncRNAs (D) Top 20 GO terms enriched by dif-circRNAs.
Competing endogenous RNA network construction
According to the dif-lncRNA–dif-miRNA pairs and dif-miRNA–dif-mRNA pairs, differentially expressed lncRNAs and mRNAs regulated by the same miRNA were screened. In total, 77 lncRNA-miRNA-mRNA interactions in lungs were finally obtained (Supplementary Fig. 1), including 35 upregulated lncRNAs and 9 downregulated lncRNAs, 5 upregulated and 5 downregulated mRNAs, and 2 upregulated and 5 downregulated miRNAs. In spleens, 131 lncRNA-miRNA-mRNA interactions were finally obtained (Supplementary Fig. 2), including 29 upregulated lncRNAs and 26 downregulated lncRNAs, 17 upregulated and 8 downregulated mRNAs, and 5 upregulated and 4 downregulated miRNAs.
Two interaction relationships of circRNA-miRNA-mRNA in lungs were obtained (Supplementary Fig. 3), comprising 2 upregulated circRNAs, 2 upregulated mRNAs, and 1 downregulated miRNA. In spleens, 32 interaction relationships of circRNA-miRNA-mRNA were obtained (Supplementary Fig. 4) including 6 upregulated circRNAs and 1 downregulated circRNA, 16 upregulated mRNAs and 2 downregulated mRNAs, 2 upregulated miRNAs and 4 downregulated miRNAs.
Further, differentially expressed circRNAs, lncRNAs, and mRNAs that were regulated by the same miRNA were further screened based on the lncRNA-miRNA-mRNA and circRNA-miRNA-mRNA analysis. Finally, 595 interaction pairs were obtained in lungs (Fig. 7), comprising 135 upregulated and 63 downregulated mRNAs, 5 upregulated and 5 downregulated miRNAs, 5 upregulated and 2 downregulated circRNAs, and 46 upregulated and 38 downregulated lncRNAs. There were 462 interaction pairs in spleens (Fig. 8), comprising 85 upregulated and 64 downregulated mRNAs, 6 upregulated and 6 downregulated miRNAs, 42 upregulated and 36 downregulated circRNAs, and 10 upregulated and 4 downregulated lncRNAs. Downregulated mmu-miR-7043-3p and Vps39–204 were significantly enriched in the ceRNA network.
Validation by qRT-PCR
The expression levels determined by qRT-PCR were in agreement with the changes in transcript abundance determined by RNA-seq analysis, which suggested that our transcriptome profiling data were highly reliable (Fig. 9).