Proteomics and metabonomics analyses of Covid-19 complications in patients with pulmonary fibrosis

Pulmonary fibrosis is a devastating disease, and the pathogenesis of this disease is not completely clear. Here, the medical records of 85 Covid-19 cases were collected, among which fibrosis and progression of fibrosis were analyzed in detail. Next, data independent acquisition (DIA) quantification proteomics and untargeted metabolomics were used to screen disease-related signaling pathways through clustering and enrichment analysis of the differential expression of proteins and metabolites. The main imaging features were lesions located in the bilateral lower lobes and involvement in five lobes. The closed association pathways were FcγR-mediated phagocytosis, PPAR signaling, TRP-inflammatory pathways, and the urea cycle. Our results provide evidence for the detection of serum biomarkers and targeted therapy in patients with Covid-19.

Chest CT findings. Of the 85 patients, 11 (12.9%) had no abnormality at admission on chest CT, including 10 (24.4%) of 41 nonpulmonary fibrosis patients and 1 (2.3%) of 44 pulmonary fibrosis patients (Table 2 and Fig. 1). Most of the patients (68 [80%]) had the single form of ground-glass opacities in the process of disease, while a few patients (10 [11.8%]) had the mixed form of ground-glass opacities and consolidation. Thirty-eight (44.7%) patients had fibrosis on admission, and 2 had fibrosis before discharge. In the nonpulmonary fibrosis group, the numbers of involved lobes were usually 1, 2, and 5. In contrast, the number of involved lobes in the pulmonary fibrosis group was mainly 5. The total lung severity score of the pulmonary fibrosis group was higher than that of the nonpulmonary fibrosis group at both stages of admission and discharge. Lesions were mainly peripheral and central distribution in the pulmonary fibrosis group (20 [50.0%]) at admission and then turned to peripheral distribution (27 [61.4%]) before discharge, whereas lesions were mainly peripherally distributed in the nonpulmonary fibrosis group. Importantly, almost half of the pulmonary fibrosis patients exhibited linear consolidation, bronchiectasis, and pleural incrassation in both stages of admission and discharge.
Then, we further evaluated the distribution of fibrosis in Covid-19 patients complicated with pulmonary fibrosis, including nonprogressive pulmonary fibrosis patients and progressive pulmonary fibrosis patients (Table 3, Fig. 1c-f.). In the progressive pulmonary fibrosis group, on admission, 6 (20%) patients had no fibrosis, and 14 (46.7%) patients had fibrosis in 5 lobes. However, 23 (76.7%) patients had fibrosis involving 5 lobes before discharge. Pulmonary fibrosis was widely distributed in both groups, especially in the left lower lobe and the right lower lobe.
Proteomic and metabolomic analysis of Covid-19 patients with pulmonary fibrosis. The Q-Exactive HF (Thermo Fisher Scientific, San Jose, CA) instrument collected mass spectrometry data from 27 samples (randomly sampled from each group) in DIA mode (Table S2). A total of 8055 peptide fragments and WBCs (× 10 9 /L; normal range 3.6-10) 5.9 (4.5-6.9) 6.0 (5.1-7.8) 6.1 (4.7-7.6) Increased 9 (10.6) 4 (9.8) 5 (11.4) Decreased 6 (7.0) 2 (2.5) 4 (9.1) Neutrophils (× 10 9 /L; normal range 1.7-7) 3 (Table S3). A heat map of the sample correlation analysis is shown in Fig. S1a. After calculating the fold difference between Covid-19 patients with pulmonary fibrosis (n = 21) and Covid-19 patients without fibrosis (n = 6) by MSstats software, 20 differentially expressed proteins were found, including 13 downregulated and 7 upregulated proteins ( Fig. 2a,b). GO and KEGG 18 enrichment analysis indicated that DEPs were related to the immune system, biological adhesion, and glycosaminoglycan degradation (Fig. S1b-d and Table S4). For metabolomic analysis, the data of positive and negative ions were collected from 28 samples (randomly sampled from each group) to improve the metabolite coverage. The coefficient of variation (CV) plot shows that there was a high degree of repeatability for samples (Fig. 3a,b). PLS-DA (Fig. 3c,d) was used to establish the relationship model between metabolite expression and sample types to predict sample categories and to determine the differential metabolism among groups. A total of 3020 positive compounds with identification information (Table S5) and 978 negative compounds were identified (Table S6), among which 102 positive compounds were upregulated and 21 were downregulated (Fig. 4a), while 38 negative compounds were upregulated and 3 were downregulated (Fig. S2a) between Covid-19 patients with pulmonary fibrosis (n = 22) and Covid-19 patients without fibrosis (n = 6). The cluster analysis of different metabolites is shown in Fig. 4b and Fig. S2b. Based on the KEGG database, pathways including the peroxisome proliferator-activated receptor (PPAR) signaling pathway, D-arginine and D-ornithine metabolism, inflammatory mediator regulation of TRP channels and alpha-linolenic acid metabolism were significantly enriched (Fig. 4c, S2c and Table S7).
Proteomic and metabolomic analyses of Covid-19 patients with progressive pulmonary fibrosis. To further search for factors related to the progression of pulmonary fibrosis, we divided the patients with pulmonary fibrosis into a progressive group and a nonprogressive group and continued to conduct proteomic and metabolomic tests on their serum. It was found that 34 proteins were differentially expressed between the progressive pulmonary fibrosis Covid-19 patients (n = 16) and the nonprogressive pulmonary fibrosis Covid-19 patients (n = 5) (Fig. 5a). After clustering (Fig. 5b) and enrichment analysis of DEPs ( Fig. 5c-e, and Table S4), it was found that DEPs were mainly distributed in the immune system, including primary immunodeficiency,    www.nature.com/scientificreports/ In terms of the metabolomic analysis, we found that the PLS-DA model had a good prediction effect (Figs. 6a and S3a). A total of 511 differential positive metabolites and 160 negative metabolites were screened (Figs. 6b and S3b) between the progressive pulmonary fibrosis Covid-19 patients (n = 16) and the nonprogressive pulmonary fibrosis Covid-19 patients (n = 6) and clustered (Figs. 6c and S3c). Using the KEGG database, metabolic pathway enrichment analysis revealed the following pathways with differential metabolite enrichments: PPAR signaling Table 3. Chest CT evaluations of Covid-19 complicated with pulmonary fibrosis patients. Except where indicated, data are cases of patients, with percentages in parentheses. *The corresponding score of each lobe in the five lobes was 0, 1, 2, 3, 4, respectively, and the corresponding affected area was 0, 1-25%, 26-50%, 51-75%, 76-100%, which was referred to the total lung severity score.   Table S7).

Correlation analysis between proteomics and metabolomics in Covid-19 patients complicated with progressive pulmonary fibrosis.
To further explore the relationship between biological network regulation, we conducted proteome and metabolomic association analyses between the nonprogressive pulmonary fibrosis group and the progressive pulmonary fibrosis group. By performing regularized canonical correlation analysis, correlation clustering of DEPs and differential metabolites can be observed in the heatmap (Fig. 7a). After multivariate dimensionality reduction and relationship visualization of omics (Fig. 7b), it was found that the DEPs and differential metabolites in pulmonary fibrosis group (n = 21) were correlated with each other.

Discussion
To our knowledge, this is the first paper describing an integrated multiomics study in Covid-19 complicated with pulmonary fibrosis patients at the molecular and ionic levels. First, after reviewing the medical records of 85 Covid-19 patients hospitalized in the People's Hospital of Guangxi Zhuang Autonomous Region Yongwu www.nature.com/scientificreports/ branch, we found that as many as 44.7% of the patients had fibrotic stripe shadows in the early onset, and the incidence rate was close to that of SARS outbreaks in 2003 19 . The most common manifestations of chest CT in Covid-19 fibrosis in our 85 cases were summarized as follows: lesions located in bilateral lower lobes, involved in five lobes, and complicated with linear consolidation, bronchiectasis, and pleural incrassation, which partially agreed with a previous report 6 . Second, we described the detailed profile of CT images in Covid-19 complicated with pulmonary fibrosis patients, especially focusing on progressive pulmonary fibrosis patients. Third, using proteomics and metabonomics methods, we explored the mechanism of the information pathway related to the occurrence and development of progressive fibrosis, which is expected to predict disease progression and to provide evidence for targeted therapy.
In recent years, proteomics has become an accurate method to screen disease biomarkers 20 and has also been used in the pathogenesis of SARS-CoV-2 infection 16,17 . In this study, we found significant differential expression levels of immune system and biological adhesion markers and glycosaminoglycan degradation in pulmonary fibrosis patients. Previous experimental and clinical data demonstrated that the accumulation and deposition of glycosaminoglycan in the extracellular matrix contributed to progressive pulmonary fibrosis 21,22 . Moreover, glycosaminoglycan hyaluronan (HA) on AEC2s (stem cells in the adult lung) is critical for restoring renewal capacity and limiting the scope of fibrosis 23 . The general properties of cell adhesion molecules (CAMs) involve cell recruitment; they participate in cell-cell interactions and cell-extracellular matrix interactions 24 and then promote the formation of fibrosis by intervention in the microenvironment 25 . For example, platelet/endothelial cell adhesion molecule-1 (PECAM-1) is expressed at intercellular junctions on platelets and leukocytes, including NK cells, T cells, monocytes and neutrophils 26 , which are involved in chemokine-mediated migration to the inflammatory site 27 . Moreover, PECAM-1 participates in the maintenance and restoration of vascular integrity following alveolar-capillary barrier disruption in ARDS 28 . In addition, CAMs modulate signal transduction by interacting with extracellular matrix proteins 29 . However, the results of this omics analysis are limited. The specific adhesion molecule biomarkers and the occurrence of diseases need further accurate detection. www.nature.com/scientificreports/ It has been well documented that the immune system plays an important role in the occurrence and development of fibrosis 30,31 and Covid-19 32 . The association of FCGR3B (encodes the lgG receptor FcγRIIIb) copy number variations with susceptibility to idiopathic pulmonary fibrosis (IPF) has been previously reported 33 . FcγR binds to serum amyloid P (SAP), inhibits the differentiation of monocytes into fibrocytes 30,31 and affects neutrophil adhesion 30 . Moreover, SAP binds to Ca 2+ -dependent ligands to induce FcγR-mediated phagocytosis in vivo 34 . Recently, it was reported that segregated-nucleus-containing atypical monocytes, which were derived from FcεRI + granulocyte/macrophage progenitors, were involved in fibrosis 35 . NK cells lacking FcεRIγ are associated with lower levels of fibrosis 36 . Furthermore, samples from 18 patients with IPF showed that the proportion of Breg cells (i.e., Breg cells to total B cells) was decreased, even though this proportion was positively correlated with the diffusion capacity of the lungs 37 . Consistent with the results of the above references, our study supported the cross-regulation among pathways of the immune system, calcium signaling and cell adhesion.
As a member of the nuclear receptor and transcription factor family, PPAR is a classical receptor for diabetic therapy 38 . However, recent emerging literature suggests a correlation between PPAR and lung fibrosis [39][40][41] . Activation of the PPARγ pathway provided protection against organ fibrosis in the lung. PPAR agonists are beneficial in treating preclinical models such as bleomycin-induced pulmonary fibrosis 42 and genetically engineered mice 43 . Moreover, inhibiting PPARγ promoted fibroproliferative ARDS, which was observed in clinical bronchoalveolar lavage fluid (BAL) samples 44 . Our results from both pulmonary fibrosis patients and progressive pulmonary fibrosis patients in metabolomic analysis confirmed the differential expression of the PPAR signaling pathway, which suggested that PPAR signaling is a core pathway in the formation and development of lung fibrosis. TRP channels respond to a series of heterogeneous stimuli and are essential for pathophysiological homeostasis 45 . It has been reported that TRP channels participate in lung disease mainly through calcium signaling, recruitment of proinflammatory cells, and neurogenic inflammatory pathways 46 . In our study, pathway enrichment revealed that inflammatory mediator regulation of the TRP channel was an important pathway to regulate the occurrence and development of pulmonary fibrosis. Both arginine and ornithine are essential components in the urea cycle. L-arginine is the precursor of nitric oxide 47 , which is involved in the vascular system 48 . Treatment with L-arginine protects against pulmonary fibrosis progression in a mouse model 49 . D-arginine is an inactive form of L-arginine. The results from our study indicated that lung fibrosis was associated with dysfunction of the urea cycle metabolic pathway. D-arginine and D-ornithine could be used as an early warning and a sign of progression of pulmonary fibrosis disease. Concentric circles of correlation between differential proteins and differential metabolites. If the angle between a differential protein and differential metabolite is an acute angle (less than 90 degrees), the correlation is positive. If the angle between the differential protein and the differential metabolite is an oblique angle (greater than 90 degrees and less than 180 degrees), the correlation is negative. Starting from the center of the circle, the lines are connected to the differential metabolites and differential proteins. The longer the connection length, the stronger the relationship. www.nature.com/scientificreports/ This study has several limitations. First, although the 85 Covid-19 cases were all confirmed patients in our hospital, their numbers are relatively few and cannot represent the whole of China or the international situation or even secondary infection. Second, the CV of proteomics was large because the samples were taken from patients' serum. Third, the immune system, PPAR, and TRP pathways screened by proteomics and metabonomics in this study need to be further verified in vivo and in vitro, while the pathways related to nervous system diseases need to be further explored, although it has been reported that Covid-19 patients have nervous system damage 50 . Fourth, long-term follow-up of the progression of fibrosis is required.

Scientific Reports
In conclusion, in this study, proteomics, metabonomics, and correlation analysis were used to screen the related pathways of fibrosis formation and progression in Covid-19 patients, and the results showed that pathways including the immune system (especially the action of FcγR-mediated phagocytosis), PPAR signaling, TRPinflammatory pathways, and the urea cycle were closely related. Our data provide new ideas for the treatment of Covid-19 and fibrosis patients. Chest CT. Eighty-five patients were imaged with a 64-detector row SOMATOM go. Top (Siemens Healthineers, Germany). The patient was in the supine position, and the images were collected at the end of inspiration. Chest CT was assessed at two time points: on admission and before discharge (or death). All images were independently reviewed by two radiologists with more than 5 years of working experience, and the final results were determined by consensus the two doctors.

Methods
Each image was evaluated for the following features: lesion morphology, fibrosis, number of lobes affected, total lung severity score, distribution of lesions, and other findings, including linear consolidation, discrete pulmonary nodules, pulmonary emphysema, bronchiectasis, pleural effusion and pleural incrassation. The total lung severity score was obtained by adding the involvement of each lobe 52,53 . The corresponding scores of each lobe in the five lobes were 0, 1, 2, 3, and 4, and the corresponding affected areas were 0, 1-25%, 26-50%, 51-75%, and 76-100%.
According to whether there were fibrosis bands on chest CT, the patients were divided into nonpulmonary fibrosis and pulmonary fibrosis groups. The latter group includes the nonprogressive pulmonary fibrosis group and the progressive pulmonary fibrosis group, based on the increased severity of fibrotic pulmonary lobe involvement.
Proteome analysis. Serum samples were placed at 56 °C (water bath) for 30 min to inactivate the SARS-CoV-2 virus and were prepared as previously described with some modifications 16 . One hundred microliters of serum sample was denatured in 900 μL of L3 buffer containing 7 M urea, 2 M thiourea and 20 mM Tris (pH 8-8.5). The serum protein was extracted by the two following reactions: dithiothreitol (10 mM) at 37 °C (water bath) for 30 min and iodoacetamide at room temperature for 30 min (dark condition). The protein was enriched by a C18 column (Waters, USA) at a flow rate of 3 mL/min, washed with 0.1% formic acid (FA) and 75% ACN, and then dried in 75% ACN. Isolated proteins were digested with 1:20 trypsin overnight at 37 °C.
Peptides were fractionated using an LC-20AB system (Shimadzu, Japan) for high-pH reversed-phase fractionation and separated on a Gemini C18 column (5 μm × 4.6 mm × 250 mm; Waters, USA) using a multistep gradient from 5 to 95% buffer B (95% ACN, pH 9.8) at a flow rate of 1 mL/min for 64 min. Eluting peptides were collected every min, which were combined with a chromatographic elution peak diagram. Then, all the fractions were combined into a total of 10 fractions, which were then frozen and dried.
The fractions were redissolved in buffer A (2% ACN and 0.1% FA) and separated on an UltiMate 3000 UHPLC (Thermo Fisher Scientific) equipped with a C18 column (150 μm i.d, 1.8-μm particles, 35 cm long; Waters, USA). Peptides were eluted by a gradient from 5 to 80% buffer B. For data-dependent acquisition (DDA) mode detection, peptides were ionized by a nanoESI source and were then introduced into a Q-Exactive HF tandem mass spectrometer (Thermo Fisher Scientific, San Jose, CA) at a spray voltage of 1.9 kV. The MS spectra (350-1500 m/z) were acquired at a resolution of 120,000, a maximum injection time (MIT) of 100 ms and an AGC target value of 3 × 10 6 . The MS 2 spectra (350-1500 m/z) were acquired at a resolution of 120,000, an MIT www.nature.com/scientificreports/ of 100 ms and an AGC target value of 1 × 10 5 . The data were identified using the Andromeda integrated by Max-Quant, and then the spectrum library was established. For data-independent acquisition (DIA) mode detection, the same nano-LC system and gradient was used as DDA analysis. Briefly, the MS spectra (400-1250 m/z) were acquired at a resolution of 120,000, an MIT of 50 ms and an AGC target value of 1 × 10 6 . The ion fragmentation mode was HCD, the MIT was selected as automatic mode, the fragment ions were detected in Orbitrap, and the resolution was set to 30,000. The fragmentation energy was distributed as follows: 22.5, 25, 27.5. The data were analyzed and controlled by the mProphet algorithm.
Metabolome analysis. The metabolites of serum samples were extracted by internal labeling and with a mixture of methanol and ACN at a ratio of 2:1. Ten microliters of supernatant from each sample was mixed into QC quality control samples to evaluate the repeatability and stability of the LC-MS analysis process. A Waters 2D UPLC (Waters, USA) series Q Exactive HF high-resolution mass spectrometer (Thermo Fisher Scientific, USA) was used for the separation and detection of metabolites.
A 5-μL aliquot of sample solution was injected on a BEH C18 column (100 mm × 2.1 mm, 1.7 μm; Waters, USA). The column oven was maintained at 55 °C. The positive ion mode mobile phases were aqueous solution containing 0.1% formic acid (A solution) and 100% methanol containing 0.1% formic acid (B solution). The negative ion mode mobile phase consisted of aqueous solution containing 10 mM ammonium formate (C solution) and 95% methanol containing 10 mM ammonium formate (D solution). The following gradient was used for elution: 0-1 min, 2% B/D solution; 1-9 min, 2-98% B/D solution; 9-12 min, 98% B/D solution; 12-12.1 min, 98-2% B/D solution; 12.1-15 min, 2% B/D solution. The flow rate was 0.35 mL/min. The range of mass spectrometry scanning with mass nucleus ratio was 70-1050, the first-stage resolution was 120,000, the AGC was 3 × 10 6 , and the MIT was 100 ms. According to the parent ion strength, top3 was selected for fragmentation, and secondary information was collected. The secondary resolution was 30,000, the AGC was 1 × 10 5 , and the MIT was 50 ms. The stepped energy was set to 20, 40, and 60 eV. The ion source was set at the following conditions: sheath gas flow rate, 40; aux gas flow rate, 10; spray voltage positive ion mode, 3.8; spray voltage negative ion mode, 3.2; capillary temperature, 320 °C; aux gas heater temperature, 350 °C.
Statistical analysis. Continuous variables are represented as medians with interquartile ranges, while categorical variables are expressed as counts (%). The statistical analyses were performed using SPSS 13.0. DDA data were identified by MaxQuant 1.5.3.30 (Max Planck Institutes) software and the UniProtKB/Swiss-Prot Homo sapiens database. Identification information that met FDR ≤ 1% was used to build the final library. The parameters were selected as follows: enzyme: trypsin; minimal peptide length: 7; PSM-level FDR: 0.01; protein FDR: 0.01; fixed modifications: carbamidomethyl (C); variable modifications: oxidation (M), acetyl (protein N-term). DIA data were analyzed by Spectronaut (Biognosys) software, which used iRT peptide to correct the retention time. Then, based on the target-decoy model applicable to SWATH-MS, the false-positive control was completed at FDR 1% to obtain significant quantitative results. The R package MSstats (Bioconductor) screened differential proteins according to the fold change > 1.5 and P value < 0.05, while screened differential metabolites according to variable importance for the projection (VIP) ≥ 1, fold change ≥ 1.2 or ≤ 0.83, and P value < 0.05. Meanwhile, differential proteins and differential metabolites were enriched by network pathway analysis database: gene ontology (GO) (http:// www. geneo ntolo gy. org.) and the kyoto encyclopedia of genes and genomes (KEGG). The protein-protein interaction network of the difference proteins was analyzed by STRING database (http:// string-db. org). DEPs and differential metabolites were calculated by regularized canonical correlation analysis and were visualized using the plotVar functions in the mixOmics package.

Data availability
The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium (http:// prote omece ntral. prote omexc hange. org) via the iProX partner repository 54 with the dataset identifier PXD024311. All data are fully available without restriction.