**2. Results**

#### *2.1. Aggravation of BLM-Induced Murine Pulmonary Fibrosis Correlated with Increased Intrinsic p38 Activity in the Lungs*

Histopathological assessment revealed that worsening severity of pulmonary fibrosis was associated with increased intrinsic p38 activity in the lungs (Figure 1A,B). At 8 days post-instillation (dpi) of BLM, tissue infiltration of inflammatory cells and thickening of the alveolar interstitium involved in aberrant collagen accumulation were observed. Moreover, these changes proceeded to di ffuse and multifocal distributions at 15 dpi. These histopathological findings of BLM-induced pulmonary fibrosis in the MKK6-constitutive active (MKK6-CA) group were more severe and extensive than those in the wild-type (WT) group, whereas those in the p38-dominant negative (p38-DN) group were less severe and extensive than those in the WT group. Moreover, the distinct severity of pulmonary fibrosis was evident in semi-quantitative evaluation assessed by a modified Ashcroft score and stratified by three mouse groups. In contrast, no apparent inflammatory and fibrotic changes were observed in phosphate-bu ffered saline (PBS)-treated groups.

In addition, total cell counts in bronchoalveolar lavage fluid (BALF) at 8 dpi of BLM tended to increase with increased intrinsic p38 activity in the lungs (Figure 1C). Regardless of mouse genotype, macrophages and lymphocytes accounted for approximately 50% and 30% of the total cells in the BALF of BLM-treated groups, respectively (Supplementary Figure S1). Similarly, comprehensive protein analysis in BALF by western blot array revealed 77 upregulated molecules that were evoked by BLM and were associated with increased p38 activity in the lungs (Supplementary Table S1). These upregulated molecules included many pro-inflammatory and pro-fibrotic mediators such as interleukin (IL)-13, IL-17, stromal cell-derived factor 1(SDF-1)/C-X-C motif chemokine ligand (CXCL)-12, interferon (IFN)-γ, keratinocyte chemoattractant (KC)/CXCL1, monokine induced by gamma interferon (MIG)/CXCL9, macrophage inflammatory protein-1 α (MIP-1 α)/CC chemokine ligand (CCL)-3, and regulated upon activation normal T cell expressed and secreted (RANTES)/CCL5.

Measurements of collagen and static compliance in the lungs also supported the morphological alterations in the three mouse groups treated with BLM (Figure 1D,E). The amount of collagen in the left lungs at 8 dpi in the MKK6-CA group was significantly higher than that in the WT and p38-DN groups, although the di fference between the WT and p38-DN groups was not significant. The decrease in static compliance in the MKK6-CA and WT groups was significantly larger than that in the p38-DN group, although the di fference was not significant between the MKK6-CA and WT groups. Additionally, a reduction in body weight tended to increase with increased p38 signaling in the lungs, exhibiting the systemic e ffect implicated in the severity of BLM-treated mice (Figure 1F).

We verified the di fferences in intrinsic p38 activity in the lungs that underlie the di fferent severity of BLM-induced fibrosis among three mouse genotypes. Immunofluorescence staining revealed the presence of AEC II, macrophages, and other parenchymal cells expressing p38 in the PBS-treated lungs (Supplementary Figure S2A). The di fferences in p38 expression in these lung cells were not observed among three mouse groups. In contrast, the proportion of the lung cells showing the nuclear localization of phospho-p38 (P-p38) was increased by BLM treatment (Supplementary Figure S2B,C). Additionally, the increased proportion corresponded with the theoretical stepwise upregulation of intrinsic p38 activity in the lungs, and this finding was most prominent in AEC II. Although p38 is ubiquitously expressed in the cytoplasm of resting cells, activated p38 is represented by the phosphorylation-dependent nuclear localization of p38 in response to various types of stimulation such as DNA damage [21,22]. Therefore, these results demonstrate the three graded intensities of p38 activation induced by BLM among three di fferent mouse genotypes.

**Figure 1.** Bleomycin (BLM)-induced pulmonary fibrosis in mice bearing three different abilities of p38 in the lungs. The three mouse genotype groups; namely, MKK6 constitutive active group (CA), wild type group (WT), and p38 dominant negative group (DN), were intratracheally administrated BLM and phosphate-buffered saline (PBS). (**A**) Representative histopathological images of the lung sections stained by Masson's trichrome (scale bar = 50 μm). Lungs were collected at 8 days post-instillation (dpi) of BLM and PBS, and 15 dpi of BLM. (**B**) Quantification of the fibrotic severity using modified Ashcroft scoring was evaluated in six different lesions at 8 dpi of BLM and PBS, and 15 dpi of BLM (n = 9). (**C**) The numbers of total cells in bronchoalveolar lavage fluid were measured at 8 dpi of BLM and PBS (n = 7). (**D**) The collagen contents of the left lung lobes were measured at 8 dpi of BLM and PBS and normalized to the weight of each left lung (n = 4). (**E**) The static lung compliances were measured at 8 dpi of BLM and PBS (n = 4). (**F**) Proportions of body weight at 8 dpi of BLM and PBS to that before administration (n = 14). All data are represented as means ± standard error of the mean (SEM). \* *p* < 0.05, n.s., no significant difference (measured by one-way an analysis of variance (ANOVA) followed by Tukey's test or unpaired Student's *t*-test).

#### *2.2. Comparative Transcriptome Analysis of a BLM-Induced Pulmonary Fibrosis Model Exhibiting Di*ff*erent Severity Due to p38 Activity in the Lungs*

RNA sequencing (RNA-seq) was performed using lung samples at 8 dpi when the severity of BLM-induced pulmonary fibrosis was apparently di fferent among the three groups and transcriptomic changes in the BLM-induced fibrosis model are more likely to be correlated with the progression of IPF [23,24]. Principal component analysis (PCA) showed a relationship in the expression of genes among the three mouse groups treated with BLM and PBS, while hierarchical clustering analysis visualized using a heatmap highlighted the trend of di fferentially expressed genes (DEGs) between the BLM- and PBS-treated groups (Figure 2A). In the PCA plot, the BLM-treated groups were all well separated from the PBS-treated groups and the variance in BLM-treated groups was less than that in PBS-treated groups, indicating the assembly of distinct clusters following BLM exposure. Consistent with this observation, hierarchical clustering analysis identified DEGs between the BLMand PBS-treated groups. Gene set enrichment analysis (GSEA) in the p38 MAPK pathway revealed that genes involved in regulating this pathway were significantly upregulated in the BLM-treated WT and MKK6-CA groups compared to those in the PBS-treated group (false discovery rate [FDR] q value < 0.25) but not in the p38-DN group (Figure 2B). Moreover, volcano plots in the three mouse groups showed that the increased number of DEGs between the BLM- and PBS-treated groups was associated with an increase in p38 signaling in the lungs (Figure 2C and Supplementary Tables S2–S4). BLM treatment upregulated approximately two-folds more DEGs and downregulated 2.5-folds more DEGs in the MKK6-CA group than those in the p38-DN group.

Next, we performed transcriptome analysis to detect the enriched functions of DEGs driven by BLM treatment. K-means clustering followed by Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis revealed the enriched pathways of four clusters in DEGs between the BLM- and PBS-treated groups of three mouse genotypes (Figure 3A). These pathways included fibrosis-related pathways such as protein processing in endoplasmic reticulum (ER) (yellow cluster), cytokine–cytokine receptor interaction (purple cluster), and ECM–receptor interaction (purple and green clusters) in addition to pathways related to immune systems such as hematopoietic cell lineage and leukocyte transendothelial migration. In contrast, we identified 493 common DEGs upregulated by BLM among the three mouse groups (Figure 3B and Supplementary Table S5). Regarding 493 common upregulated DEGs, enrichment analysis by gene ontology (GO) revealed three ECM-related annotations among the top five enriched molecular functions and all five terms associated with immune systems among the top five enriched biological processes. Additionally, the cytokine–cytokine receptor interaction was the most significantly enriched pathway among the top five KEGG pathways.

**Figure 2.** *Cont.*

**Figure 2.** Expression profiling of BLM- and PBS-treated lungs in three different mouse genotypes. Three samples from each group were sequenced. (**A**) Principal component analysis of RNA sequencing datasets among the BLM- and PBS-treated three mouse groups (left). Hierarchical clustering shown in a heatmap of gene expression profiles between the BLM- and PBS-treated groups (right). The red and blue strips represent upregulated and downregulated genes in each group, respectively. (**B**) Gene set enrichment analysis of differential expression in the p38 MAPK pathway between the BLMand PBS-treated groups of three mouse genotypes. The normalized enrichment scores (NES), normal *p*-values, and false discovery rate (FDR) q values are indicated. (**C**) Volcano plot of differentially expressed genes altered by BLM treatment among three mouse groups. Upregulated and downregulated genes are discriminated based on log2 fold-change and adjusted FDR *p*-value (<0.05).

**Figure 3.** *Cont.*

**Figure 3.** Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of differentially expressed genes (DEGs) altered by BLM treatment among three mouse groups. (**A**) K-means clustering shown in a heatmap based on the gene expression profiles between the BLM- and PBS-treated groups of three mouse genotypes (left). The colored bars on the right of the diagram indicate clusters. The trees represent enriched KEGG pathways corresponding to each cluster (right). (**B**) Venn diagram showing the overlap of DEGs upregulated by BLM among three mouse groups, with the numbers of DEGs indicated in each area (left). GO and KEGG pathway enrichment analysis of 493 common upregulated DEGs among the three genotypes (right). Top five enriched GO terms associated with molecular function (upper) and biological process (middle), and KEGG pathway analysis (bottom).

#### *2.3. Exploration of Novel Potential Genes Contributing to the Progression of Pulmonary Fibrosis*

To identify the pathogenetically relevant genes in the progression of IPF, we investigated the correlation of upregulated genes between BLM-induced fibrotic lungs showing the three different severity levels and human IPF lungs (Figure 4). In the BLM-treated groups, K-means clustering analysis identified a cluster of 2722 genes that their mean reads per kilobase of exon per million mapped sequence reads (RPKM) values increased along with stepwise elevation of p38 signaling in the lungs (Supplementary Table S6). We verified 137 DEGs that were included in this cluster and upregulated in common with the three BLM-treated mouse groups. Additionally, human RNA-seq data that provided 475 upregulated DEGs in IPF lung tissue compared to healthy lung tissue was obtained from the Gene Expression Omnibus website (accession ID: GSE52463). Finally, comparison of our data with human RNA-seq data identified four overlapping DEGs; namely, EPH receptor A3 (*EPHA3*), POU class 2 homeobox associating factor 1 (*POU2AF1*), SAM domain, SH3 domain and nuclear localization signals 1 (*SAMSN1*), and ectodysplasin A2 receptor (*EDA2R*).

**Figure 4.** Identification of potential target genes by comparison to a publicly available idiopathic pulmonary fibrosis (IPF) dataset. K-means cluster analysis among three mouse groups treated with BLM revealed a cluster of 2722 genes showing correlation between variations of their mean expression values and stepwise changes in intrinsic p38 activity in the lungs (**left upper**). The Venn diagram in the right upper tier shows the overlap of 2722 genes in this cluster and 493 common upregulated DEGs among the three mouse groups. Likewise, the Venn diagram in the right middle tier shows the overlap of 137 genes identified in our study and 475 upregulated DEGs in human IPF lungs from dataset GSE52463. The four overlapping genes identified in these analyses were EPH receptor A3 (*EPHA3*), POU class 2 homeobox associating factor 1 (*POU2AF1*), SAM domain, SH3 domain and nuclear localization signals 1 (*SAMSN1*), and ectodysplasin A2 receptor (*EDA2R*) (**bottom**). Each bar and plot represent mean reads per kilobase of exon per million mapped sequence reads (RPKM) ± SEM and RPKM value of each sample, respectively.
