Deconvolution of RNA-Seq Analysis of Hyperbaric Oxygen-Treated Mice Lungs Reveals Mesenchymal Cell Subtype Changes

Hyperbaric oxygen (HBO) is widely applied to treat several hypoxia-related diseases. Previous studies have focused on the immediate effect of HBO-exposure induced oxidative stress on the lungs, but knowledge regarding the chronic effects from repetitive HBO exposure is limited, especially at the gene expression level. We found that repetitive HBO exposure did not alter the morphology of murine lungs. However, by deconvolution of RNA-seq from those mice lungs using CIBERSORTx and the expression profile matrices of 8 mesenchymal cell subtypes obtained from bleomycin-treated mouse lungs, we identify several mesenchymal cell subtype changes. These include increases in Col13a1 matrix fibroblasts, mesenchymal progenitors and mesothelial cell populations and decreases in lipofibroblasts, endothelial and Pdgfrb high cell populations. Our data suggest that repetitive HBO exposure may affect biological processes in the lungs such as response to wounding, extracellular matrix, vasculature development and immune response.


Introduction
Hyperbaric oxygen (HBO) refers to the inhalation of pure oxygen under greater than one atmosphere pressure, and is now widely used as a treatment for hypoxia-related diseases and some ischemic injuries [1]. Despite HBO being frequently administered in various clinical situations, it has several side effects, the main one being oxygen toxicity, with the risk of developing toxicity increasing with pressure and duration of treatment [1]. The lungs are the primary sites of oxygen toxicity, as a result of their direct exposure to the high partial pressure of oxygen during HBO treatment. Several conditions in pre-term and neonatal infants have been reported which are linked to hyperoxia, including hyperoxia-induced impairment of alveolarization, pulmonary inflammation, bronchopulmonary dysplasia complicated with pulmonary hypertension suggesting the developing lung is more sensitive to hyperoxia [2][3][4][5][6]. Conversely, the impact of hyperoxia on the adult lung

Repetitive HBO Exposure Has Little Impact on Mice Lung Morphology and Body Weight
To assess the effect of repetitive HBO exposure on mice lungs, mice were exposed to 2.5 ATA (Atmosphere Absolute) HBO for 90 min, once a day, for 11 consecutive days. This protocol is consistent to those used in clinic, with the difference being a 2-day rest following 5-day exposure, due to weekends in clinical application. Mouse lungs were collected on day 12. Hematoxylin-eosin (H/E) staining of mice lungs showed little difference between the control group and HBO-treated mice, HBO-treated mice displayed well organized alveolar structure similar to control group mice ( Figure 1A,B). Lung injuries can cause fibrosis, and this is characterized by the accumulation of collagen in the lungs, Masson's trichrome staining was performed to visualize the collagen distribution. No obvious accumulation of collagens was observed in HBO-treated mice ( Figure 1C,D). Mice were weighed every three days, as an indicator of the general health of the mouse, and no significant difference was observed between HBO and control group ( Figure 1E).

Repetitive HBO Exposure Shows Mild Effects on Gene Expression in Mice Lungs.
Although we did not observe changes in morphology, we postulated that the gene expression profiles of HBO treated lungs may be altered before changes in morphology and function are observed. Therefore, RNA-sequencing (RNA-seq)-based analysis was conducted to detect the effect of repetitive HBO exposure on gene expressions of lung tissues in an unbiased manner. The samples were collected at the same time point as the morphological observation to rule out the immediate effect of HBO exposure. Four samples of each group were sequenced. A total of 179,623,190 reads of the control group (44,905,797 reads per sample on average, ranging from 38,614,624 to 53,412,322 reads) and 174,456,292 reads of the HBO group (43,614,073 reads per sample on average, range from 39,350,919 to 47,038,994 reads) passed the quality control. 93.66% of overall reads were mapped to the GRCm38 reference. Using the criteria described in the methods, 16,692 genes with more than 10 counts across the eight samples were identified for further analysis.
DEGs were analyzed using DEseq2 [24]. When screened with a commonly used rigid threshold padj < 0.05 (adjusted p-value), only 33 genes were differentially expressed, of which 15 genes were up-regulated and 18 genes were down-regulated (Table S1). Seven out of 15 genes were up-regulated more than 2-fold, including Spink5, Nr1d2, Dbp, Nmrk2, Ddit4, Nr1d1, mt-Atp6. 14 out of 18 genes were down-regulated more than 0.5-fold: Ctrb1, Spon2, Pnlip, Cela2a, Cpa1, Igkv6-14, Prss2, Rag1, Try4, Try5, Ighv5-9, Cel, Ighv3-8, Amy2b. Some of the abundantly expressed genes were validated through real time qPCR, and the expression change of Ddit4, Spink5, mt-ATP6, mt-CO3, Spon2 and Prss2 was consistent with RNA-seq results ( Figure S1). When the threshold was adjusted to p < 0.05 (p-value), 636 differentially expressed genes were selected, of these 364 were up-regulated while 272 were down-regulated (Figure 2A), the differentially expressed genes are listed in Table S2. To compare the differences between control and HBO groups, we calculated the covariance of all samples and built a matrix that represents the overall features of the data for principal components analysis (PCA). The first and second component represented 31.5% and 20.7% of the variance within the transcriptomic data respectively. In the PCA plot, HBO samples were not all well separated from control samples, indicating the difference between groups is less significant than that within groups, the influence of HBO exposure is minimal ( Figure 2B).
To examine processes and pathways which may be altered upon HBO-exposure, GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis based on the p < 0.05 DEGs was conducted. Several biological processes were identified, including ATP metabolism, proton transport and response to oxidative stress. Further signaling pathways related to these processes were shown to be enriched in the up-regulated genes ( Figure 2C,D), which was consistent with previous studies [9,11,25]. No other processes identified appeared to be directly related to the influence of extra oxygen and oxidative stress. In the down-regulated genes, enriched terms for cellular component classification included "collagen containing extracellular matrix" and "extracellular matrix component", suggesting that the extracellular matrix may be affected.
Since the gene expression changes were minimal when compared the HBO exposure group and control group, gene set enrichment analysis (GSEA) [26] may be a more appropriate method to extract biological insight from the RNA-seq data. Through this analysis, biological processes like protein translation, mitochondrial respiratory chain, redox homeostasis, detoxification and immune related processes were up-regulated, while processes related to extracellular matrix were down-regulated (Table 1 and Table S3), which was consistent with GO enrichment analysis. Notes: GO process: the gene set name; NES: normalized enrichment score; FDR q-value: false discovery rate. The significant gene sets were selected by the threshold FDR q-value < 0.05, and only the top 30 up-regulated gene sets were shown, while all the down-regulated gene sets were listed.

Repetitive HBO Exposure Has Little Effect on the Relative Composition of Different Cell Types of Mice Lungs
The cellular composition of the lung tissue varies according to developmental stages, and between normal and pathological status [27]. To further analyze the influence of repetitive HBO exposure on the cell composition of lung tissue, we utilized a recently developed machine learning method called CIBERSORTx analysis [28]. The single-cell RNA sequencing (scRNA-seq) data of adult mouse lung from the Mouse Cell Atlas (MCA) project [29] were used to determine the gene expression signature of different cell types, in which 32 cell types were classified. We successfully reconstructed the cell clusters to be consistent with those reported [29], which was displayed by cell distributions according to their expression profile in two-dimensions by t-distributed stochastic neighborhood embedding (t-SNE) projections ( Figure 3A). We demonstrated that the parameters we used were consistent with the original findings. Furthermore, the mean Pearson r value was 0.967 when comparing our transcriptome data with the scRNA-seq data ( Figure 3B), demonstrating a good correlation and confidence of the analysis. Transcriptome data were then dissected by CIBERSORTx and we found there was no statistically significant difference in the ratio of any type of the cells ( Figure 3C).

Deconvolution Based Upon Fibrotic Lung Cell Profiles Reveals Mesenchymal Cell Subtype Changes
Repetitive micro-injuries to lung as well as abnormal repair may result in extracellular matrix deposition in damaged lung foci and this can result in the development of pulmonary fibrosis [30,31]. Mesenchymal cells are understood to play a major role, as the origins of secretary collagen and thus be key in disease pathogenesis [30,31]. Given the fact that HBO-induced oxidative stress may cause minor injury to the lungs, we utilized another scRNA-seq dataset which focused on mouse lung mesenchymal cells, from uninjured and bleomycin-treated mouse lungs, classified and defined six mesenchymal cell types in the normal lung and seven in the fibrotic lung [32]. To determine the composition of the mesenchymal cells in the transcriptome data of lung tissues in HBO exposure and control groups, data were deconvoluted using CIBERSORTx analysis. We first reconstructed the eight cell clusters identified in the study [32] (Figure 4A), based on the scRNA-seq data of the seven subtypes of fibroblast as well as endothelial cells from fibrotic lung. Then by CIBERSORTx analysis, the transcriptome data of lung tissues in HBO exposure and control groups were dissected to determine the composition of the mesenchymal cells. Mean Pearson r value of our dataset and the scRNA-seq data was 0.897, suggesting a high confidence level of this analysis ( Figure 4B). The inferred ratios of 6 out of 8 the cell types were found to be statistically different (p < 0.05) in the two groups. Among these, the ratios of Col13a1 matrix fibroblast, mesenchymal progenitor and mesothelial were up-regulated, while those of lipofibroblast, mesothelial, and Pdgfrb high cells were down-regulated when comparing HBO to the control group ( Figure 4C).

GO Enrichment Analysis Identifies Different Roles of the Eight Mesenchymal Cell Types
To identify functional consequence of the changes in the mesenchymal cell populations following repetitive HBO exposure, GO enrichment analysis of biological process using the DEGs in each of the mesenchymal cell types was performed.
Several biological processes were identified; we focused on the non-redundant top 10 GO terms enriched (q < 0.01) in the biological process (all the enriched terms were shown if the number of enriched terms was less than 10) (Figure 5A-F and Figure S2). The full enrichment results are shown in Table  S4-S11. Among the up-regulated cell types, the distinct expression genes in Col13a1 matrix fibroblasts were enriched in the terms "response to toxic substance", "extracellular structure organization" and "cell matrix adhesion" ( Figure 5A). However, the term "response to toxic substance" was not enriched in Col14a1 matrix fibroblast ( Figure S2A), which is another type of Col1a1 high expressing matrix fibroblast, and its abundance was not changed after repetitive HBO exposure. GO terms enriched in mesothelial were mainly related to development ( Figure 5B), and those in mes progenitor, the potential mesenchymal progenitors, were related to cell proliferation ( Figure 5C).
On the other hand, among the down-regulated cell types, the enriched GO terms in the lipofibroblasts were focused on the immune response ( Figure 5D). Previously, lipofibroblasts have been reported to be involved in alveolar development, regeneration, surfactant synthesis and vitamin A (retinoic acid) storage [33], and contribute to the alveolar epithelial cell II (AECII) stem cell niche in the adult mouse lung [34]. The enrichment result in endothelial cells was mainly related to vasculature development ( Figure 5E), which is consistent with the knowledge of endothelial cells as vasculature forming cells [35]. Several varied terms were identified among the enriched biological processes of Pdgfrb high (Pdgfrb hi) cells, suggesting that this cell type is a transitional cell type. Pdgfrb hi cells are reported to only be present in the fibrotic lung [32], suggesting it plays an important role during fibrotic pathogenesis.
Finally, the GO enrichment analysis showed that Col14a1 matrix fibroblast might contribute to extracellular matrix organization and response to growth factor, while myofibroblast showed characteristics of response to fibroblast growth factor, and muscle tissue development which is consistent with previous reports [36]. Neither of the two cell types were altered after repetitive HBO exposure.
To compare the potential contribution of different cell types during fibrotic pathogenesis and response to oxygen, we selected four related biological processes and visualized the expression level of the top 15 significant DEGs of the 8 cell types annotated in these four processes ( Figure 5G). We demonstrate that four cell types underlie these biological processes: lipofibroblasts, endothelial cells, and Col13a1 and Col14a1 matrix fibroblasts. Lipofibroblasts captured a signature related to immune response and response to wounding; endothelial cells showed high expression of genes related to vasculature development and response to wounding. Meanwhile, Col13a1 and Col14a1 matrix fibroblasts expressed high level of genes related to extracellular matrix and response to wounding, suggesting their correlation to these biological processes. Among them, lipofibroblast and endothelial cells were down-regulated, and Col13a1 matrix fibroblasts were up-regulated. Together, these results suggested that repeated HBO exposure may have an influence on these four biological processes.

Discussion
According to the central dogma, changes in the gene expression profile occur earlier than those in morphology and functions [37]. Therefore, the present study was carried out to assess the non-instant effects of repetitive HBO-exposure on the gene expression profile of lung tissues through RNA-seq. Although morphologically normal, small changes in lung tissue gene expression were detected after repeated HBO exposure by RNA-seq followed by qPCR validations. Processes such as protein translation, mitochondrial respiratory chain, redox homeostasis, detoxification and immune related processes were up-regulated while that of extracellular matrix were down-regulated through GO enrichment and GSEA. It is critical to understand which cell type(s) is (are) important to mediate these molecular changes. The recently developed scRNA-seq technique can be used to delineate the contributions of different cell subtypes. Several studies based on the scRNA-seq of mouse lung tissue have already reported on the complexity of lung composition [29,32,[38][39][40]. Alternatively, some bioinformatic methods, such as CIBERSORT, can be used to characterize cell composition of complex tissues from their bulk transcriptome data [41] and have previously been successfully used to enumerate immune subsets in the tumor microenvironment [42][43][44][45][46][47]. CIBERSORTx, an updated version of CIBERSORT, is a machine learning method that can infer cell-type-specific gene expression profiles without physical cell isolation, as well as dissect large-scale tissue data using scRNA-seq data [28]. In the present study, by deconvolution of RNA-seq from HBO-exposed mice lungs using CIBERSORTx and the expression profile matrices of eight mesenchymal cell subtypes obtained from bleomycin-treated mouse lungs [32], we identify several mesenchymal cell subtype changes. These include increases in Col13a1 matrix fibroblasts, mesenchymal progenitors and mesothelial cell populations and decreases in lipofibroblasts, endothelial and Pdgfrb high cell populations.
Lung fibroblasts are a heterogeneous group of cells, producing extracellular matrices to build the lung structural framework, while also directing cell growth by cell-cell and cell-matrix contact and/or paracrine mechanisms [30,31,48]. Based on different gene expression profiles and expected functions, subsets of fibroblasts, including myofibroblasts, Col13a1 matrix fibroblasts, Col14a1 matrix fibroblasts, lipofibroblasts, mesenchymal progenitors, and mesothelial cells, and Pdgfrb hi fibroblast subpopulation, have been delineated [32].
Based on our analysis, the distinct expression genes in Col13a1 matrix fibroblast were enriched in response to toxic substances, and this cell type expanded in fibrotic lung [32], suggesting that it was an important contributor to fibrotic pathogenesis. Its ratio increased in the HBO group suggesting that repetitive HBO exposure might induce some toxic stress. Lipofibroblasts have been reported to contribute to AECII stem cell niche [34] and we showed HBO exposure resulted in reduction of lipofibroblasts, which could be detrimental for AECII stemness maintenance. HBO treatment has been shown to alter inflammation [1]; consistent with that, we found that lipofibroblasts acquire a gene signature involved in immune response and the reduction of this cell type might cause inhibition of inflammation. Chronic wounds exist in a perpetually inflammatory state [1], as such, the inhibition of inflammation might be beneficial. A recent study showed Pdgfrb high cells to be responsive to fibrotic injury [32], reduced levels of these cells following HBO exposure may be protective. A decrease in endothelial cells was observed after HBO exposure, these are a major player in vasculature development. Although HBO treatment promotes neovascularization [1], this effect seems to be tissue specific, as HBO exposure may reduce the blood vessel density in some cancer models [49,50]. Furthermore, both high and low oxygen concentrations have been reported to cause pulmonary vascular remodeling in neonatal infants [5,51]. We suggest that the reduction in the endothelial cell population might be an adaptive response to extra oxygen.
Generally, we found the results from CIBERSORTx analysis, traditional DEG analysis and GSEA to be consistent. Reports have demonstrated that HBO can induce oxidative stress [25]; we found increased expression of Ddit4 (REDD1), a gene induced in response to various stresses [52], the enrichment of gene sets like cell redox hemostasis and detoxification in HBO-group through GSEA, and the increase of Col13a1 matrix fibroblasts, which capture a gene expression profile of response to toxic substance. All of which suggested repetitive HBO exposure induced stresses to lung tissue, consistent with the previous report [25]. Another example is the influence of HBO exposure on inflammation. Spon2 (Mindin), which encodes an integrin ligand and pattern recognition molecule that is critical for inflammatory cell recruitment [53,54], was down-regulated by about 50% in the HBO group. Igkv6-14, Ighv5-9, and Ighv3-8, which all encode an immunoglobulin component (through searching in NCBI gene database), were down-regulated significantly in HBO group. Spink5, which encodes the serine protease inhibitor Lympho-epithelial Kazal-type inhibitor (LEKTI) and plays an anti-inflammatory role [55], is increased in the HBO group. All of these changes in gene expression suggest that HBO exposure inhibits inflammation. The gene set of "negative regulation of immune response" was enriched in HBO group suggesting an inflammation inhibition effect too. The immune inhibition information we get from the DEG analysis and GSEA are consistent with the reduction in lipofibroblasts, which characteristically show activation of inflammation, and consistent with the reports that HBO exposure attenuate inflammation in a few of conditions [56][57][58][59][60][61].
However, there are also some inconsistencies among CIBERSORTx, DEG and GSEA analysis. For example, the GO term "extracellular matrix" was enriched in the downregulated genes through GSEA; however, Col13a1 matrix fibroblast was upregulated in the HBO group, while Col14a1 matrix fibroblast remains unchanged, the total ratio of Col13a1 and Col14a1 matrix fibroblast increased. Both of them are cell types expressing signature genes associated with extracellular matrix and cell adhesion [32], suggesting that the extracellular matrix was upregulated. Several studies have reported that HBO promotes fibroblast proliferation and/or extracellular matrix synthesis in vitro and in vivo [1,[62][63][64][65]. In addition, the effect of HBO on angiogenesis has been reported during HBO treatment of chronic wounding and tumor tissue, although its effect is context-dependent [1,49,50]. This was evidenced by CIBERSORTx analysis, but not through DEG and GSEA analysis. More importantly, CIBERSORTx analysis provides information on potential effects of fibroblast differentiation. Compared to other deconvolution methods, one advantage of CIBERSORTx is that it considers and removes the excessive technical variation in signature matrices from single-cell references and bulk RNA expression [66]. Luca et al. [67] applied CIBERSORTx to impute cell type-specific gene expression profiles of cancer cells in The Cancer Genome Atlas (TCGA). They successfully identified stable transcriptional states of 12 cell types, almost all of them can be validated in published single-cell datasets, revealing a comprehensive insight of cellular states in the tumor microenvironment. Shi et al. [68] showed overall consistency between the result of CIBERSORTx and their novel probabilistic deconvolution model which used to predict the contribution of exosomal sources on the tissue microenvironment. The accuracy of deconvolution predictions was further validated at both the pre-and on-treatment time points. However, since the studies on the heterogeneity of fibroblast subtypes in mouse lung are recently emerging, our understanding of fibroblast subpopulations, of cell markers at the protein level, their localization in the lung, signaling programs, and functional significance are limited. Further identification of "pathogenic" fibroblast subpopulations in lung disease will enable us to understand the precise effects of HBO on lung functions.

Animals
Six-to eight-week-old male C57BL/6 mice (20-25 g) were provided by the Animal Center of Nantong University. Mice were raised under a 12 h light/dark cycle, and normal diet and water were provided ad libitum throughout the study. During the study, the mice were weighed every third day. All the animals were provided by the Experimental Animal Center of Nantong University (Institutional license: SYXK(SU)-2012-0030). All procedures involving animals were approved by the Institutional Animal Care and Use Committee (IACUC) of Nantong University (NTU), Nantong, Jiangsu Province, China (approval number: 20140901-001).

HBO Exposure
For hyperbaric oxygen exposure, a hyperbaric chamber designed for small animal research was used. The exposure regimen consisted of a 5-min pressure ramp-up to 2.5 ATA (1.5 atm) in a 100% O 2 environment, sustained for 90 min at this pressure, and followed by a 5-min decompression phase. To maintain a high concentration of oxygen, the chamber was flushed twice for 5 min with pure oxygen before HBO exposure. The concentrations of oxygen and carbon dioxide were monitored in real time by SDA oxygen and carbon dioxide monitors (ANALOX, North Yorkshire, England), respectively. Twenty mice were randomized into control and HBO groups. The mice in the HBO group were subjected to the hyperbaric oxygen exposure once a day for 11 consecutive days, while the mice in the control group were placed in the chamber for the same duration without pure oxygen pressurization.

Morphological Examination
Twenty-four hours after the last HBO exposure, the mice of both groups were sacrificed under anesthesia with composited anesthetic (with 257 mM chloral hydrate, 176 mM magnesium sulfate, 36 mM pentobarbital sodium, 14.25% ethanol, and 33.8% propylene glycol), and the whole lung tissues were collected. The lung tissues were fixed in 4% paraformaldehyde and dehydrated through gradient concentrations of ethanol. The fixed lung tissues were embedded in paraffin, sequentially sectioned into 5 µm sections, and then stained with H/E solution to observe the morphology of lung tissue, and Masson's trichrome staining to detect the presence of collagen. The selected sections were viewed and imaged using a Leica DM4000B microscope (Wetzlar, Germany).

RNA Isolation, Library Construction and Sequencing
RNA was isolated with Trizol (Invitrogen, Carlsbad, California, USA). A total amount of 3 µg RNA per sample was used as input material for library construction. Sequencing libraries were generated using NEBNext ® UltraTM RNA Library Prep Kit for Illumina ® (NEB, Ipswich, Massachusetts, USA) following manufacturer's instruction. Libraries were pooled in equimolar and sequenced using the paired-end strategy (2 × 150) on the Illumina NovaSeq 6000 platform following the standard protocols. RNA-seq data have been deposited in the Gene Expression Omnibus (GEO) database (accession code GSE143348).

Real-Time qPCR Analysis
Real-time qPCR was done to validate the RNA-seq results. Total RNA was reverse transcribed through HiScript III RT SuperMix for qPCR (+ gDNA wiper) (Vazyme, Nanjing, China), and the qPCR was done with ChamQ Universal SYBR qPCR Master Mix (Vazyme, Nanjing, China). Primers used were as follow: Ddit4

Single-Cell Data Analysis
Single-cell datasets were downloaded from GEO (Accessions GSE108097 [29] and GSE104154 [32]). We used Seurat [73] (version 3.0) to filter cells, select the number of principal components, perform t-distributed stochastic neighbor embedding (t-SNE) reduction according to the description in the original papers, ensuring the same classification of cell types. Numerous biomarkers of defining clusters were identified through FindVariableGenes function.

CIBERSORTx Analysis
We used CIBERSORTx [28] to estimate the percentage of various cell populations in each bulk-RNA sample. Signature matrix of each cell type was generated following the default settings. The deconvolution of bulk RNA-seq using the signature matrix with S-batch correction removed variances between different sequence platforms. The permutation value for obtaining robust statistical result was set as 1000. Two-tailed, unpaired Student's t-test was used to compare significant differences in abundances of member cell types between two conditions. p-values were adjusted for multiple testing using Benjamini-Hochberg method.

Pathway Enrichment Analysis
Gene ontology and KEGG pathway enrichment analysis were generated through Metascape [74] website (http://metascape.org). Parameters were set as five minimum overlap genes, p < 0.05 and 1.5 minimum enrichment factor. GSEA (gene set enrichment analysis) [26] was performed using GSEA software (the Broad Institute platform, version 4.0.3) with the default settings and 1000 gene set permutations. It was performed on normalized counts when comparing HBO samples versus control samples. All gene sets with an FDR (false discovery rate) q-value less than 0.1 were considered to be significant.

Statistical Analysis
Each experiment was repeated at least twice. Each transcriptomic condition had four replicates. Unless otherwise noted, data are presented as mean and standard deviation. Two-tailed, unpaired Student's t-test was used to determine statistical significance through Prism 8. p < 0.05 was considered as statistically significant.

Conclusions
Taken together, this study suggests a greater focus may be required on the lung status when using long-term, repetitive HBO treatment especially in relation to wounding, extracellular matrix organization, vasculature development and immune responses, processes that become aberrant in the development of chronic lung diseases. Further experiments are required to confirm its impact.
Supplementary Materials: Supplementary Materials can be found at http://www.mdpi.com/1422-0067/21/4/ 1371/s1. Figure S1: qPCR validation of differential expressed genes. Figure S2: GO enrichment analysis of the differentially expressed genes of Col14a1 matrix cells and myofibroblasts. Table S1: Differentially expressed genes (padj < 0.05). Table S2: Differentially expressed genes (p < 0.05). Table S3: GSEA analysis. Table S4: GO enrichment of differentially expressed genes in Col13a1 matrix cells. Table S5: GO enrichment analysis of differentially expressed genes in Col14a1 matrix cells. Table S6: GO enrichment analysis of differentially expressed genes in endothelial cells. Table S7: GO enrichment analysis of differentially expressed genes in lipofibroblasts. Table S8: GO enrichment analysis of differentially expressed genes in mesothelial cells. Table S9: GO enrichment analysis of differentially expressed genes in mesprogenitors. Table S10: GO enrichment analysis of differentially expressed genes in myofibroblasts. Table S11: GO enrichment analysis of differentially expressed genes in Pdgfrb high cells.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.