SARS-CoV-2 accessory proteins involvement in inflammatory and profibrotic processes through IL11 signaling

SARS-CoV-2, the cause of the COVID-19 pandemic, possesses eleven accessory proteins encoded in its genome. Their roles during infection are still not completely understood. In this study, transcriptomics analysis revealed that both WNT5A and IL11 were significantly up-regulated in A549 cells expressing individual accessory proteins ORF6, ORF8, ORF9b or ORF9c from SARS-CoV-2 (Wuhan-Hu-1 isolate). IL11 is a member of the IL6 family of cytokines. IL11 signaling-related genes were also differentially expressed. Bioinformatics analysis disclosed that both WNT5A and IL11 were involved in pulmonary fibrosis idiopathic disease and functional assays confirmed their association with profibrotic cell responses. Subsequently, data comparison with lung cell lines infected with SARS-CoV-2 or lung biopsies from patients with COVID-19, evidenced altered profibrotic gene expression that matched those obtained in this study. Our results show ORF6, ORF8, ORF9b and ORF9c involvement in inflammatory and profibrotic responses. Thus, these accessory proteins could be targeted by new therapies against COVID-19 disease.


Introduction
The coronavirus disease 2019 (COVID-19) is a potentially fatal respiratory disease caused by the new Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2), which rapidly spread worldwide causing more than 670 million reported cases and nearly 7 million deaths globally since the start of the pandemic (https://coronavirus.jhu.edu/map.html). The clinical course of COVID-19 exhibits a broad spectrum of severity and progression patterns. While the infection leads to mild upper respiratory disease or even asymptomatic sub-clinical infection in a significant number of people, others develop symptoms and complications of severe pneumonia that can be fatal. Furthermore, pulmonary fibrosis has been described as one of the most common consequences in COVID-19 patients, even in long COVID-19 (1)(2)(3)(4)(5). Indeed, Fabbri et al. estimated that approximately 20% of patients with COVID-19 had evidence of fibrotic sequels one year after viral infection (6). Since March 2020, many efforts have been done to elucidate COVID-19 pathogenesis, but the complete clinical picture following SARS-CoV-2 infection is not yet fully understood.
SARS-CoV-2 ORF6 is a 61 aa protein that localizes in endoplasmic reticulum and membrane of vesicles such as autophagosomes and lysosomes (9,16). This accessory protein displays multifunctional activities such as blocking nucleopore movement of newly synthetized mRNA encoding immunemodulatory cytokines such as IFN-b and interleukin-6 (IL-6) counteracting those cytokines (17). SARS-CoV-2 ORF8 is a 121 aa protein consisting of an N-terminal signal sequence for endoplasmic reticulum (ER) import. It is a secreted protein, rather than being retained in the ER, and its extracellular form has been detected in the supernatant of cell cultures and sera of COVID-19 patients (9,18). In addition, ORF8´s functions are mediated by its binding to CD16a, decreasing the capacity of monocytes to mediate antibody-dependent cellular cytotoxicity (ADCC) (19). SARS-CoV-2 ORF9b is a 97 aa protein that antagonizes type I and III interferons by negatively regulating antiviral immunity (20). It is localized in the mitochondrial membrane associated with TOM70 (13) inducing proinflammatory mitochondrial DNA release in inner membranederived vesicles (21). SARS-CoV-2 ORF9c is a 73 aa membraneassociated protein that suppresses antiviral responses in cells (22). It also interacts with Sigma receptors that are implicated in lipid remodeling and ER stress response (9,23). SARS-CoV-2 mostly affects the respiratory tract usually leading to pneumonia in most patients, and to acute respiratory distress syndrome (ARDS) in 15% of cases. ARDS is mainly triggered by elevated levels of pro-inflammatory cytokines, such as Interleukin 6 (IL6), referred to as cytokine storm (24). Interleukin 11 (IL11) is a member of the IL6 family of cytokines. IL11 is similar to IL6, and both form a GP130 heterodimer complex to initiate its downstream signaling (25-28), but their respective hexameric signaling complex formation differ (29). While IL6R is expressed most highly on immune cells, IL11RA is expressed in stromal cells, such as fibroblasts and hepatic stellate cells, and also on parenchymal cells, including hepatocytes. Hence, it may be expected that IL6 biology relates mostly to immune functions whereas IL11 activity is more closely linked to the stromal and parenchymal biology (25, 30-32). Since the nineties, high IL11 release during viral infections have been described (33,34), and more recently, several studies have related this interleukin to fibrosis, chronic inflammation and matrix extracellular remodeling (31,(35)(36)(37)(38)(39). It is also known that WNT5A and IL11 have the ability of activating STAT3 signaling (40) and this ability has been postulated as a possible mechanism to link WNT5A gene with immunomodulation. WNT5A is a member of WNT family proteins which plays critical roles in a myriad of processes in both health and disease, such as embryonic morphogenesis, fibrosis, inflammation or cancer (41). Several studies have described a crosstalk between transforming growth factor-beta (TGFb) and WNT signaling pathways during fibrotic processes (42- 45), and more recently with the increase in IL11 production (46). TGFb represents the most prominent profibrotic cytokine by upregulating production of extracellular matrix (ECM) components and multiple signaling molecules (47).
It is known that the underlying cause of severe COVID-19 disease is a cytokine dysregulation and hyperinflammation status (24,48,49), and IL6 was from the beginning involved as it was found to be elevated in serum of COVID-19 patients (50, 51). However, little is known about the involvement of IL11 in lung fibrosis in COVID-19 disease. In this study, A549 lung epithelial cells were individually transduced with accessory proteins ORF6, ORF8, ORF9b or ORF9c from SARS-CoV-2 (Wuhan-Hu-1 isolate), and transcriptomic analysis revealed that both, WNT5A and IL11, were significantly up-regulated. IL11 signaling-related genes, such as STAT3 or TGFb, were also differentially expressed. Subsequently, bioinformatics and functional assays revealed that these four accessory proteins were implicated in both inflammatory and fibrotic responses in different extents, suggesting the involvement of ORF6, ORF8, ORF9b and ORF9c in inflammatory and/or fibrotic responses in SARS-CoV-2 infection.

Immunofluorescence microscopy
Cells were seeded on 24-well plates containing glass coverslips coated with poly-lysine solution (100.000 cells per well). Cells were fixed with 4% PFA in PBS for 15 min, washed twice in PBS, and then permeabilized for 10 min with 0.1% Triton X-100 in PBS. Primary antibodies incubation was carried out for 1h in PBS containing 3% BSA and 0.1% Triton X-100 at 1:100 dilution. Coverslips were washed three times with PBS before secondary anti-mouse antibodies incubation (1:1000 dilution). The antibodies used for immunofluorescence are shown in Supplementary Table 5  (Table S5). Phalloidin was used as a cytoplasmic marker at 1:200, and DAPI (4'6-diamidino-2-phenylindole) (Thermo Fisher Scientific, #62248) was used as a nuclear marker. Coverslips were mounted in Mowiol 4-88 (Sigma-Aldrich, #81381). Images were acquired with a confocal laser microscope Leica TCS SP8 STED 3X.

RNA isolation and sequencing
Cells were seeded (3 x 10 5 ) in 6-well plates and lysed using RLT buffer for RNA isolation (RNeasy mini kit, Qiagen, #74106). Each sample was performed in triplicate. RNA was isolated following manufacturer´s protocol, quantified by nanodrop 1000 (Thermo Scientific) and quality controlled by Bioanalyzer (Agilent). All samples sent for sequencing had a RIN (RNA integrity number) over 9.90. cDNA libraries and sequencing were performed by Novogene Europe, using 400 ng of RNA per sample for library preparation. Samples were sequenced in an Illumina platform using a PE150 strategy.

Gene sets and differential gene expression analysis
Sequencing raw data was quality controlled (error rate, GC content distribution) and filtered, removing bad quality and N-containing sequences and adaptors. Clean data were mapped (HISAT2) to reference genome GRCh38.p13, and gene expression was quantified using FPKM (Fragments Per Kilobase of transcript sequence per Millions of base pairs sequenced). Differential expression analysis was performed using DESeq2 R package (53).
Raw counts were transformed with the vst function in the DESeq2 package (54) of the R software version 3.6.3 (55), and subsequent PCA was performed with the prcomp function. The 500 genes with the highest variance among samples were considered. Finally, the PCA graph was made with GraphPad Prism 5 (GraphPad software, San Diego, CA, USA).

Real time qPCR
RNA samples (500 ng) were reverse transcribed using qScript ™ cDNA synthesis kit (Quanta Biosciences Inc., #95047), following manufacturer's instructions. Primers sequences are available in Supplementary

Western blot
Transduced cells were harvested and lysed in ice-cold Pierce IP Lysis Buffer (Thermo Scientific, #87787) at 4°C. Cell lysates were mixed with 5× SDS-PAGE Sample Loading Buffer (Nzytech, MB11701), and heated at 95°C for 5 min. Protein samples were resolved by SDS polyacrylamide gel electrophoresis and transferred onto a PVDF membrane using Mini Trans-Blot System (Bio-Rad, #1703935), followed by blocking for 1 h with 5% BSA in Trisbuffered saline-Tween20 buffer and probing with corresponding primary and secondary antibodies (Table S5). The proteins were visualized by chemoluminiscence using ChemiDoc Imaging Systems (Bio-Rad). Relative protein expression was calculated by sequentially normalizing against the loading control (GAPDH).

Bazedoxifene treatment and ELISA
Cells were seeded (3 x 10 5 ) in 6-well plates and treated with 5 µM of Bazedoxifene acetate (Sigma, PZ0018) for 24h. To perform ELISA experiments, IL-11 levels in supernatants collected after 24 h treatment from different cell lines were detected with the Human DuoSet ELISA Kits (RD Systems, DY218) according to the manufacturer's instructions.

Cell contraction assay
CytoSelect ™ 24-well Cell Contraction Assay Kit (Cell Biolabs, CBA-5020) was used according to the manufacturer's instructions. Briefly, collagen gel lattice was prepared by mixing 4.5 x 10 6 cells/ mL with a collagen gel solution and added to each well of the 24well cell contraction plate. After collagen polymerization, fresh media was added and wells were monitored for contraction over two days at 37°C and 5% CO2. The change in matrix diameter size (in millimeters) was determined with a ruler each 24h.
Bioinformatics analysis. Database comparison. Pathway enrichment analysis, network and ppi module reconstruction.
Functional pathway analysis of transduced cells was performed with Ingenuity Pathway Analysis (IPA) software. Adjusted p-value less than 0.05 was considered as the cut-off criterion for pathway enrichment analysis. To compare our results in A549 lentivirustransduced expressing individual viral accessory proteins ORF6, ORF8, ORF9b or ORF9c with whole virus-infected cell lines A549-ACE2 or Calu3, transcriptomic data from (56) were used. Subsequently, transcriptomic data from (57) were applied for patient samples comparison. Bioinformatics analysis were carried out making Venn diagrams with Venny 2.1 (58), and heatmaps with Heatmapper program. A further enrichment study was performed with DAVID Functional Annotation Tool where selected genes were clustered according to GO Terms and Reactome Gene Sets. Additionally, another pathway enrichment analysis and the gene network reconstruction were carried out using the online Metascape Tool (http://metascape.org) (59) with the default parameters set. Enrichment analyses were carried out selecting the genomics sources: KEGG Pathway, GO Biological Processes, Reactome Gene Sets, Canonical Pathways, and CORUM. Terms with p < 0.01, minimum count 3, and enrichment factor >1.5 were collected and grouped into clusters based on their membership similarities. P-values were calculated based on accumulative hypergeometric distribution, and q-values were calculated using the Benjamini-Hochberg procedure to account for multiple testing. To further capture the relationship among terms, a subset of enriched terms was selected and rendered as a network plot, where terms with similarity >0.3 are connected by edges. Based on Protein-Protein Interaction (PPI) enrichment analysis, we run a module network reconstruction based on the selected genomics databases. The resulting network was constructed containing the subset of proteins that form physical interactions with at least one other list member. Subsequently, by means of Molecular Complex Detection (MCODE) algorithm, we first identified connected network components, then a pathway and process enrichment analysis were applied to each MCODE component independently and the three best-scoring (by p-value) terms were retained as the functional description of the resulting modules.

Quantification and statistical analysis
Statistical analyses were performed using GraphPad PRISM 5. P-values were determined using two-way ANOVA and Bonferroni test correction was applied. Unless otherwise stated, data are shown as the mean of at least three biological replicates. Significant differences are indicated as: *, p <0.05; **, <0.01; ***, p<0.001, ****, p<0.0001.

Results
Expression of SARS-CoV-2 ORF6, ORF8, ORF9b or ORF9c accessory proteins alter gene expression pattern in A549 cells SARS-CoV-2 uses several strategies to interact and interfere with the host cellular machinery. To explore the function of individual ORFs in such interaction, A549 human lung carcinoma cells were lentivirus transduced expressing individual viral accessory proteins ORF6, ORF8, ORF9b or ORF9c (named ORF-A549 thereafter), with a C-terminally 2xStrep-tag to facilitate detection of their expression ( Figure 1A). GFPlentivirus-transduced or wild-type A549 cells were used as control in each experiment, both giving the same results. ORFs overexpression in A549 transduced cells was verified by immunofluorescence staining using anti-StrepTag antibody which highlighted different patterns of localization in A549 cells, as well as variable levels of ORF expression between cells in each cell line ( Figure 1B). ORF9b and ORF9c seemed to be highly concentrated around the nucleus, while ORF6 and ORF8 were localized mainly in a specific perinuclear region.
Differential gene expression analysis was performed in ORF-A549 cells (Figure 2A). Sample quality control was assessed by principal component analysis (PCA) based on normalized counts from DESeq2. High quality was achieved since samples were clustered ( Figure 2B). Further analysis of transcriptomics data revealed a number of genes commonly expressed in all transduced cells, including WNT5A and IL11 whereas they were not as upregulated in other ORF-A549 (data not shown). These two genes were particularly upregulated, as well as other genes previously related to their signaling pathways (40, 41, 60) ( Figure 2C). qRT-PCR was used to validate transcriptomic data in ORF-A549 cells for CXCL1, IL11, WNT5A, WNT5A-AS and STAT3 ( Figure 2D). Also, IL11 release was significantly increased in cells expressing ORF8, ORF9b and ORF9c ( Figure 2E).

A549 transduced cells differently express genes involved in pulmonary fibrosis idiopathic signaling
Based on the above results, a functional pathway analysis with Ingenuity Pathway Analysis (IPA) software was performed. Both WNT5A and IL11 related canonical pathways were selected and two canonical pathways were found in common between the four transduced cell lines: Cardiac Hypertrophy Signaling and Pulmonary Fibrosis Idiopathic Signaling ( Figure 3A). Genes involved in pulmonary fibrosis of each transduced cell line were obtained and further analysis showed a high gene expression pattern similarity between ORF-A549 cells ( Figure 3B). Subsequent qRT-PCR experiments corroborated differential gene expression for collagen genes such as COL1A1, COL4A1 or COL11A1, or other genes like ADAMTS1, BCL2, IL1B, MMP16, SERPINE1, SNAI1 or TFGB1 ( Figure 3C). To assess whether altered expression of these genes could affect the profibrotic behavior of the cells, a functional assay was performed to test the ability of ORF-A549 cells to contract a collagen matrix ( Figure 3D). After 24h, ORF6 and ORF9b expressing cells were able to significantly shrink the collagen matrix, while ORF8 and ORF9c transduced cells were able to do it only after 48h. Surprisingly, the contractile capacity of ORF9b-A549 cells was significantly higher than the others.
Inhibition of IL11 signaling pathway modulates the effect of ORF6, ORF8, ORF9b and ORF9c expression in A549 cells Involvement of IL11 with fibrosis has been shown in previous studies (31, 35, 37, 61). To decipher IL-11 involvement, ORF-A549 cells were treated with an IL11 receptor inhibitor: Bazedoxifene (BAZ). ORF-A549 cells were treated with 5 mM BAZ during 24h and expression levels of various genes involved in fibrosis were determined ( Figures 4A-G). To facilitate data interpretation, expression changes associated with BAZ treatment in each ORF cell line were represented but they were not related to A549 control cells. 100% of expression was considered for ORF-A549 BAZ untreated cells, and BAZ treated expression was calculated with respect to untreated cells in each case. A549 control cells expression nor expression differences of each ORF cell line related to A549 control cells was depicted ( Figures 4A-G). A decrease in IL11 expression levels was observed in ORF8, ORF9b and ORF9c expressing cells, but not in ORF6-A549 cells ( Figure 4A). These results agreed with those observed for IL11 secretion measured by ELISA ( Figure 4H). As expected, we found changes in WNT5A after BAZ treatment ( Figure 4B), but they were cell dependent. After IL11 signaling inhibition by BAZ, only ORF8-A549 cells showed a decrease in WNT5A expression. Interestingly, ORF8-A549 cells had the smallest increase in WNT5A when validated by RT-qPCR ( Figure 2D). By contrast, ORF9b-A549 cells increased WNT5A expression after BAZ treatment, but no changes were observed in ORF6 or ORF9c expressing cells. Surprisingly, we did not observe any change in TGFb expression in any ORF-A549 cells ( Figure 4C). These results suggest that IL11 involvement in such profibrotic processes might be not mediated by TGFb signaling. A decrease in SERPINE1 expression after BAZ treatment was observed, particularly in ORF8, ORF9b and ORF9c expressing cells ( Figure 4D). On the other hand, a significant increase in IL1B, SNAI1 and ADAMTS1 expression was observed in ORF9b-A549 cells after BAZ treatment (Figures 4E-G). No changes in expression of these genes were shown in ORF6, ORF8 or ORF9c cell lines, suggesting a crosslink between IL11 and IL1B signaling pathways in ORF9b-A549 cells.
IL11 increase after viral infections (33,34) and a relationship between IL11 and WNT5A through STAT3 pathways signaling has been previously described (40). Therefore, STAT3 phosphorylation after IL11 signaling inhibition was analysed by western blot. A significant reduction in STAT3 phosphorylation was observed in cells expressing ORF8 and ORF9c ( Figures 5A, C). It is reported that activation of the TGFb signaling cascade causes phosphorylation and activation of the cytoplasmic effectors such as Smad2 (62). However, we did not observe changes in TGFb expression nor Smad2 phosphorylation in any ORF-A549 cells ( Figures 5A, D, E). These results were consistent with those observed by qRT-PCR ( Figure 4C), where no changes in TGFb expression were observed. Once more, these results suggest that IL11 involvement in this process may not be TGFb dependent. Surprisingly, we did not observe significant changes in WNT5A expression after BAZ treatment (Figures 5B, F). A significant decrease of WNT5A expression was found in ORF9c-A549 cells, but BAZ treatment did not alter such expression ( Figure 5F). Regarding SERPINE1, a reduction in its expression by cells expressing ORF6 and ORF9c after BAZ treatment was observed, but it was only significant in ORF6-A549 (Figures 5B, G). Interestingly, a significant increase of phosphorylated c-jun in Frontiers in Immunology frontiersin.org cells expressing ORF9b and ORF9c was found. In addition, BAZ treatment reduced phosphorylated c-jun in these cell lines ( Figures 5B, H). However, ORF6 and ORF8 cell lines did not show changes in phosphorylated c-jun, and even BAZ treatment significantly augmented phosphorylated c-jun in cells expressing ORF6. Given the fact that expression of these accessory proteins modified their profibrotic capacity, IL11 involvement was analysed by BAZ treatment inhibiting IL11 signaling pathway in the collagen contraction assay (Figures 5I-L). Interestingly, all ORF-A549 cells were able to revert the effect of expressing ORF6, ORF8, ORF9b or ORF9c accessory proteins. After 24h of treatment, we did not find changes in ORF6 and ORF9c cells compared to control cells ( Figures 5I, L), but we did in cells expressing ORF8 and ORF9b ( Figures 5J, K). By contrast, after 48h of BAZ treatment, all ORF-A549 cells recovered similar levels of collagen area when compared with untreated control cells. Therefore, these data indicate that IL11 signaling pathway is directly related to the profibrotic capacity described in ORF-A549 cells.

Profibrotic response of lung epithelial cells to SARS-CoV-2 accessory proteins resemble responses to whole virus infection
In order to investigate the relevance of these profibrotic processes in SARS-CoV-2 virus infection, a bioinformatics comparative study was performed by integrating transcriptomic results from SARS-CoV-2 infected lung cell lines or COVID-19 lung biopsies with those obtained in this study. The aim was to analyse common genes differentially expressed and their possible relationship with host fibrotic response when the whole virus was present. To this end, we grouped the sets of fibrosis-related genes in . Statistical significance is given as follows: *p < 0.05, **p < 0.01, ***p < 0.001 ****p < 0.0001 to A549 control cells.
ORF-A549 cells obtained by IPA analysis, and a single common list of 63 fibrosis-related genes was generated ( Figure 3B). Subsequently, our pre-identified differential expression data list was compared with those obtained from infecting ACE2transfected A549 cells and Calu3 cells with SARS-CoV-2 (NCBI-GEO, GSE147507) (56) ( Figure 6A). Interestingly, we found 4 common genes between the three lung cell lines (IL11, SNAI1, COL4A1 and COL4A2), as well as 4 common genes between our ORF-A549 cells and ACE2-transfected A549 cells (COL11A1, COL21A1, COL5A2 and COL6A1), and 3 common genes between our ORF-A549 cells and Calu3 cells (SERPINE1, THBS1 and MUC1). Gene expression disclosed two genes commonly upregulated (IL11 and SNAI1) among lung cell lines, except in the case of ORF6-A549 cells, where SNAI1 was not differentially expressed ( Figure 6A). When we compared our list of fibrosis-related genes with transcriptomic data from post-mortem COVID-19 lung biopsies (https://github.com/Jiam1ng/COVID-19_Lung_Atlas) (57), 28 common genes among two data lists were found ( Figure 6B). Further analysis of gene expression revealed 4 genes commonly downregulated (COL4A4, COL4A3, WNT9A and COL21A1) among lung biopsies and ORF-A549 cells, except in the case of ORF6-A549, where COL4A3 and COL21A1 were not differentially expressed. At the same, 10 genes were found commonly upregulated, nevertheless, only four of them were upregulated by ORF6-A549 cells (SERPINE1, CDH2, F2 and IL11). Once again, ORF6-A549 cells had the fewest genes in common with the other cell lines and lung biopsies. These data agreed with those previously shown, for example, in terms of heatmap clustering, PCA or IL11 secretion, which show the difference of ORF6-A549 with the rest of ORF-A549 cells (Figures 2A, B, E).
To further investigate differential perturbation of pathways regulated by ORF6, ORF8, ORF9b and ORF9c accessory proteins in SARS-CoV-2 infection, the 28 common genes list was used to perform an enrichment study with DAVID Functional Annotation Tool, where selected genes were clustered according to GO Terms and Reactome pathways ( Figure 6C). We obtained the most statistically significant pathways involved in fibrosis and  Figure 6C). These results were in line with preliminary results obtained by proteomics analysis, in which three enzymes involved in the formation of collagen fibers were found altered (data not shown). All represented terms showed 5 common genes between our ORF-A549 cells and COVID-19 lung biopsies (COL1A1, COL4A3, COL4A4, COL5A1 and COL11A1). Apart from that, other commonly upregulated genes were found in certain terms, such as IL11, IL1B, SERPINE1, SNAI1 and JUN. Interestingly, both cytokines, IL11 and IL1B were localized in Statistical significance is given as follows: *p < 0.05, **p < 0.01 ***p < 0.001 ****p < 0.0001. In all cases data are represented as mean ± SD (n=3). extracellular space. Later, additional enrichment analysis and PPI network were obtained with MCODE network components (Metascape) ( Figure 6D). Top two best p-value terms were retained: MCODE1, related to genes encoding collagen proteins (NABA_COLLAGENS, M3005, -log10 pval=46,7), and MCODE2, related to WNT signaling (M5493, -log10 pval=10,2). Although we did not observe any MCODE component clustering genes such as IL11, IL1B, SERPINE1 or SNAI1, we did notice a relationship between these genes in the PPI network. Interestingly, JUN appeared as a connecting node of this cluster of genes.
A549 lung epithelial cell line has been used in other studies as a model to study SARS-CoV-2 virus infection (56,64,65). In those studies, A549 cells were modified to express ACE2 receptor to facilitate the virus entry into the cell. In this study, A549 lung epithelial cells were individually transduced with lentivirus encoding accessory proteins ORF6, ORF8, ORF9b or ORF9c from SARS-CoV-2 (Wuhan-Hu-1 isolate), and subsequent transcriptomic with bioinformatic analysis disclosed that these accessory proteins can be involved in inflammatory and/or fibrotic responses in SARS-CoV-2 infection. To perform that approach, SARS-CoV-2 viral infection was not required, so ACE2 expression in A549 cells was not performed. Noteworthy, virulent strains such as MERS-CoV, SARS-CoV and SARS-CoV-2 have a significant number of these accessory proteins, while more harmless coronaviruses have less (66, 67). This suggests that accessory proteins play a key role in pathogenesis not observed in less virulent coronavirus infections, although they have been less characterized than other proteins contained in the viral genome. Importantly, mutations in accessory proteins ORF6, ORF8 and ORF9b have been observed in currently circulating SARS-CoV-2 "variants of concern", thus potentially contributing to increasing pathogenesis and transmissibility (https://covariants.org/variants).
At the beginning of the pandemic, high levels of IL6 in COVID-19 patient serum were described to correlate with severe disease (50, 51). Surprisingly, we found a high overexpression of IL11 in epithelial transduced cells ( Figure 2C), while no changes in IL6 expression or release were observed (data not shown). The validation of IL11 by qPCR was partially consistent with the RNAseq results ( Figure 2D), probably due to the different sensitivity of both procedures. Indeed, IL11 release was also increased in three of the four transduced cell lines ( Figure 2E), which corresponds to those with high levels of IL11 expression in RNAseq. These results agree with those reported in several studies, where IL11 was defined as an "epithelial interleukin", while IL6 biology was related mostly to immune functions (25, 29, 31). In addition, several studies have related high levels of IL11 with fibrosis, chronic inflammation and matrix extracellular remodeling (31, 35-39). However, whether this elevation is pathogenic or a natural host response to restore homeostasis remains unanswered for many diseases (36).
We also found several fibrosis related genes differentially expressed. Among them, WNT5A was particularly upregulated in ORF6 and ORF9b expressing cells. Likewise, we found an increase in WNT5A-AS (Figures 2C, D). WNT5A is a member of WNT family proteins which plays critical roles in a myriad of processes in both health and disease (41), and it is known its relationship with IL11 through STAT3 pathways signalling (40). Besides, chemokines CXCL1 and CXCL12 have been found to be upregulated by WNT5A in various studies (41, 60). These results are consistent with those we have observed in transcriptomic analysis, where both chemokines were upregulated together with WNT5A ( Figures 2C, D), although they did not exactly correlate with ORF6 and ORF9b expressing cells. WNT5A-AS has been reported as a long noncoding RNA (lncRNA) located on the antisense strand of chromosome 16 p16, and which overlaps with introns of WNT5A on the sense strand (68). Indeed, Lu et al. and Salmena et al. provided evidence of a positive correlation between the upregulation of lncRNA WNT5A-AS with that of its antisense gene, WNT5A, suggesting that lnRNA WNT5A-AS acts as a competing endogenous RNA to regulate the expression of WNT5A (68,69). In this study, we found high levels of WNT5A gene expression in ORF6 and ORF9b expressing cells, but we did not observe such increase in protein expression (Figures 5E, F). Thus, it is plausible to think that WNT5A-AS was responsible for regulating posttranscriptional expression of WNT5A.
TGFb represents the most prominent profibrotic cytokine by upregulating production of ECM components and multiple signaling molecules (47). There is, indeed, a clear evidence of the relationship between IL11, WNT, TGFb and fibrosis (42-46). In this study, TGFbwas among the genes related to fibrosis, and we observed an increase in TGFb gene expression but our assays did not find any upregulation in its protein expression (Figures 5A, D). When canonical TGFb signaling pathway through Smad2 phosphorylation was analysed, no significant changes were also observed ( Figures 5A, C), indicating that TGFb involvement may be regulated either by non-canonical TGFb signaling pathway or through TGFb-independent manner.
Within the list of genes involved in fibrosis, we also found genes such as ADAMTS1, BCL2, CDH2, FN1, IL1B, JUN, MMP16, SERPINE1, SNAI1 and various collagen genes ( Figures 3B, C). Expression of these genes differed depending on the expressed accessory protein. In organs, such as the lungs, resident cells actively and continuously remodel the extracellular matrix (ECM), forming a dynamic network balanced by cell-ECM bidirectional interactions (70). In order to check how A549 transduced cells responded to and actively remodeled the ECM, we performed a collagen gel contraction assay. A strong ability in ORF9b-A549 cells and moderate ability in ORF6-A549 cells to contract the collagen gel in the first 24 hours was observed, which was also followed by ORF8 and ORF9c transduced cells after 48 hours. These data indicate that expression of these accessory proteins, and particularly ORF9b, trigger a profibrotic process. In addition, preliminary proteomics analysis revealed three altered enzymes involved in collagen fibers formation: PLOD1, PLOD2 and COLGALT1 (data not shown). It is known that PLOD1 and PLOD2 catalyze the lysyl hydroxylation to hydroxylysine, which is critical for the formation of covalent cross-links and collagen glycosylation (71). COLGALT1 acts on collagen glycosylation and facilitates the formation of collagen triple helix (72), and an increase of WNT5A gene expression has been recently correlated with COLGALT1 downregulation (73). These three enzymes tend to be downregulated in ORF-A549 cells, except in ORF6-A459 cells, where PLOD2 was significantly increased. Furthermore, PLOD1 was significantly decreased in ORF8-A549 cells, while PLOD2 and COLGALT1 were significantly decreased in ORF9b-A549 cells, meaning a possible involvement of these enzymes in the increased collagen-contraction ability of these cells (data not shown).
Our hypothesis was that IL11 might be behind the above mentioned profibrotic alterations in ORF-A549 cells, so we used an IL11 receptor inhibitor to block IL11 signaling pathway. Several studies have identified BAZ as a novel small-molecule inhibitor of GP130 (74), and support its therapeutic action targeting IL-11/ GP130 signaling for cancer therapy (75,76). BAZ binds to GP130 heterodimer and inhibits IL6 family members-induced STAT3 phosphorylation (74), blocking interleukins signaling pathways without affecting their release, as we observed ( Figure 4H). In this study, BAZ treatment of ORF-A549 cells reverted their high c o l l a g e n -c o n t r a c t i o n a b i l i t y ( F i g u r e s 5 I -L) . S T A T 3 phosphorylation was also reduced in ORF8 and ORF9c expressing cells (Figures 5A, B). Noteworthy, ORF8 and ORF9c expressing cells showed the lowest levels of WNT5A gene expression, being ORF8-A549 the lowest one ( Figure 2D). Indeed, WNT5A expression only decreased in ORF8-A549 cells after BAZ treatment ( Figure 4B). These data suggest a possible WNT5A signaling pathway regulation by IL11. Inhibiting IL11 signaling, cells with high expression levels of WNT5A may compensate the effect of BAZ, whereas ORF8-A549 may not. We also observed a decrease in IL11 release by all transduced cells ( Figure 4H) and a downregulation in SERPINE1 expression ( Figure 4D), but we did not notice a significant downregulation of others genes altered by accessory proteins expression, such as TGFb, IL1B, SNAI1 or ADAMTS1 (Figures 4C-G). Interestingly, an increase in IL1B, SNAI1 and ADAMTS1 in ORF9b expressing cells after BAZ treatment was found, suggesting a possible crosslink between IL11 and IL1B signaling pathways. Palmqvist et al. provided evidence of enhanced IL11 expression by IL1B by a mechanism involving MAPK in gingival fibroblasts (77). Nevertheless, further investigations must be performed to test this hypothesis and its relationship with high levels of SNAI1 and ADAMTS1 when blocking IL11 signaling pathway.
Finally, data comparison with lung cell lines infected with SARS-CoV-2 and lung biopsies from patients with COVID-19 showed evidence of altered gene expression that matched with results obtained in this study. Firstly, we found common differentially expressed genes with SARS-CoV-2 infected lung cell lines (56). Among these genes, IL11 was commonly upregulated, as well as SNAI1. On the other hand, 28 common genes related with fibrosis were found between our transduced cells lines and COVID-19 lung biopsies (57). We also found IL11 between this cluster of genes. A subsequent enrichment analysis showed that this set of genes was mostly involved in ECM organization and collagen formation ( Figure 6C). Interestingly, both IL11 and IL1B were located in extracellular space. These data were consistent with the fact that we did not find these interleukins in our proteomics study or by western blot assay (data not shown). Five collagen genes were common in all signaling pathways (COL1A1, COL4A3, COL4A4, COL5A1 and COL11A1) and they were also clustered together by MCODE algorithm using Metascape tool ( Figure 6D). Remarkably, JUN was listed in the cellular component organization GO term (Figure 6C), and it was also found in PPI network as a gene connecting node ( Figure 6D). That pointed to a possible c-jun role connecting profibrotic cell responses. Ser73 c-jun phosphorylation was confirmed by western blot in ORF9b and ORF9c expressing cells ( Figures 5E, H). Indeed, this c-jun activation decreased after BAZ treatment. Noteworthy, phosphorylation of cjun increased in ORF6-A549 cells after blocking IL11 signaling pathway, indicating that c-jun activation in ORF9b-A549 and ORF9c-A549 cells was mediated by IL11 expression. C-jun Ser73 is phosphorylated by MAPK8 (78), and JNK-interacting proteins (JIP) are a scaffold proteins group that selectively mediates JNK signaling by aggregating specific components of the MAPK cascade. Among JIP proteins, SPAG9 or JIP4 (also known as MAPK8IP4) is involved in MAPK signaling pathway to regulate cellular activities (79). Del Sarto et al. recently have identified an increase of phosphorylaton at Ser730 in JIP4 after Influenza A virus infection (80). Similarly, other work has recently described that this phosphorylation promotes cell death via c-jun kinase signaling pathway (81). Thus, there is a correlation between c-jun phosphorylation and SPAG9 (JIP4) phosphorylation after Influenza A virus infection. Preliminary phosphoproteomics analysis in this study revealed the presence of phosphorylated form of SPAG9 in the position Ser730 (data not shown). However, further mass spectrometry analysis will be required to define the phosphorylation patterns and abundance changes of phosphorylated SPAG9 in ORF-A549 cells.
On the other hand, ORF6-A549 cells showed the lowest levels of IL11 release, differentially of the rest of ORF-A549 set of cells. Data from ORF6-A549 appeared separated from the other ORF-A549 cells in transcriptomic clusters ( Figure 2B) and showed less genes in common when compared with SARS-CoV-2 infected lung cell lines and COVID-19 lung biopsies (Figures 6A, B). Indeed, only ORF6-A549 cells showed an increase in PLOD2 enzyme, and a decrease in SERPINE protein expression after BAZ treatment. All these data indicate a different mechanism of action by ORF6 that require further investigations.
It is true that fibrotic processes involve more cell populations than just epithelial cells. However, several studies focused on fibrosis have reported the direct involvement of IL11 secretion by lung epithelial cells after an injury, either in lung inflammation or stimulating fibroblast phenotype transformation in adjacent epithelial cells (25, 35, 37), and infer that IL11 signaling may represent a therapeutic target in lung fibrosis (25). This data suggests and important role of this cytokine early in the fibrosis process which might be considered as a target to counteract COVID-19 fibrotic sequels.
Taken together, our findings indicated that SARS-CoV-2 accessory proteins ORF6, ORF8, ORF9b and ORF9c have the ability to trigger profibrotic cell responses in A549 human lung epithelial cells in vitro (Figure 7). This process was particularly evident in cells expressing ORF9b, which showed a greater effect in extracellular matrix remodeling. Interestingly, increased IL11 led to ECM remodeling. Data from SARS-CoV-2 infected lung cell lines and COVID-19 lung biopsies from patients show a similar response to SARS-Cov-2 infection, so these profibrotic responses may underlie COVID-19-fibrotic alterations. Thus, it is plausible to think that these accessory proteins might be used as a target for new therapies for COVID-19 disease against pulmonary fibrosis.

Limitations of the study
This study has several limitations that should be addressed in future studies. Even though the model used in this work allowed us to study individual SARS-CoV-2 proteins and their interactions with cellular components, the possible co-regulation of viral proteins and the effect this may have on the host cell have not been assessed. Nevertheless, the reductionist model carried out enabled us to gain knowledge of specific viral genes functions that may be difficult to unravel in more complex approaches.
A further limitation was the extrapolation of our in vitro results to the lung biopsy samples. Given the fact that access to solid lung biopsies from patients with COVID-19 was not possible, an in silico comparison with patient data that had been obtained in other studies was chosen for validation. Such a comparison would not be appropriate using other types of biopsies, such as bronchoalveolar lavage (BAL) due to the small percentage of epithelial cells present in these biopsies. Considering this limitation, we found it relevant that almost half of the identified profibrotic genes were found differentially regulated in the data from patient samples. Another related limitation is the fact that we used a single gene list in the bioinformatic comparison with those obtained from lung biopsies. The enrichment analysis naturally showed fibrosis-related pathways, as we expected. However, we found it interesting to focus on these pathways and genes involved, highlighting the similarities between our in vitro epithelial model, and the epithelial environment from lung biopsies, as this is a different model of response to infection than the one referred to in other immunological studies. In addition to this, fibrotic processes involve more cell populations than just epithelial cells. However, the particular secretion of the interleukin IL11 by these cells, and its known relationship with these pro-fibrotic processes, suggest an important involvement of these cells early in the process that may then involve adjacent cell populations.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/ bioproject/, PRJNA946640 for A549 transduced cells and PRJNA841835 for A549 control cells.