MSCs-derived secretome counteracts the effect of apoptosis inhibition by SARS-CoV-2 S protein on Calu-3 cells mediated lowered expression of PTGS2

Cell culture and apoptosis detection

Human airway epithelial Calu-3 cells (Procell Life Science & Technology Co., Ltd.,China) were cultured at 37℃ with 5% CO2 in α-MEM with 10% fetal bovine serum (FBS), 100 µg/mL streptomycin, 100 U/mL penicillin (All supplements were from Gibco). The cells used for the experiments were 5–8 passages.

Calu-3 cells were treated with recombinant SARS-CoV-2 spike protein for 48 hours at various final concentrations (0, 1, 10, 50, 100, and 500 nM). The recombinant protein used in this study (Sino Biological, Cat PSV001, Sino Biological) corresponds to the receptor-binding domain (RBD) of the S1 subunit. It is expressed in HEK293 cells and purified to >95% purity as verified by SDS-PAGE and high-performance liquid chromatography by the manufacturer. This RBD fragment lacks the S2 subunit and the transmembrane domain, and is therefore provided as a monomeric form that does not recapitulate the native trimeric spike structure. All experiments were performed using a single production lot of this recombinant protein to ensure batch-to-batch consistency. After the 48-hour treatment, some wells were additionally exposed to 0.75 μM actinomycin D (A4262, Sigma) for another 12 hours.

After the indicated treatments, cells were harvested for apoptosis detection with flow cytometry (Cytoflex, Beckman) using double-stained with APC-conjugated Annexin V and 7-aminoactinomycin D (4A Biotech Co. Ltd., Beijing, China) according to the manufacturer's protocol. For all flow cytometry plots, axis limits were uniformly set to include the entire recorded cell population, and the same scaling was applied across all experimental groups to facilitate accurate visual comparison. All flow cytometry experiments were carried out with three independent biological replicates. The apoptosis rate was calculated as the percentage of Annexin V-positive cells (early plus late apoptotic populations) in the total cell population, based on quadrant statistics using the instrument’s analysis software.

EDU assay

Calu-3 cells were seeded in 48-well plates, cultured for 48h. Cell proliferation assay was performed using the BeyoClick™ EdU Cell Proliferation Kit (Beyotime, C0081S). Briefly, various concentration of S protein-treated Calu-3 cells were incubated with 10 μM EdU for 2 hours at 37°C. After fixation, permeabilization and rinse, the cells were exposed to 100 μL of click reaction cocktail for 30 min. The percentage of EdU-positive cells was analyzed by flow cytometry and defined as the proliferation rate

²⁵

.

Western Blot

First, 40 mg total protein lysate was separated on a 10% polyacrylamide gel and then and then transferred to a polyvinylidene difluoride (PVDF) membrane. The membranes were blocked in 5% skim milk for 1 h, then were incubated overnight at 4°C with antibodies against PTGS2 and GAPDH (Cell Signaling, Danvers, MA, USA). Second, the membranes were then washed three times with PBST and incubated with the secondary horseradish-peroxidase (HRP)-conjugated anti-rabbit antibody (Boster, Wuhan, China) at room temperature for 2h. Third, protein bands were exposed to autoradiography film for visualization using enhanced chemiluminescent imaging system (Tanon, Shanghai, China).

RNA extraction, library construction and sequencing

The selected treatment Calu-3 cells were subjected to total RNA extraction using TRIZOL (Ambion). The RNA was further purified with two phenol-chloroform treatments and then treated with RQ1 DNase (Promega, Madison, WI, USA) to remove DNA. The quality and quantity of the purified RNA were determined by measuring the absorbance at 260 nM/280 nM (A260/A280) using Smartspec Plus (BioRad, USA). The integrity of RNA was further verified by 1.5% agarose gel electrophoresis. For each sample, 1 μg total RNAs were used for stranded RNA sequencing library preparation using Ribo-Zero Gold rRNA Removal Kit (human/mouse/rat) (N406-01, Vazyme) and remained RNAs were fragmented and purified, and converted into double strand cDNA. After end repair and A tailing, the DNAs were ligated to VAHTS RNA Adapters (Vazyme). After purification, ligation products corresponding to 200 ~ 500 BPS were digested with heat-labile UDG, and the single strand cDNA was amplified, purified, quantified, and stored at -80℃ before sequencing. For high-throughput sequencing, the libraries were prepared for following the manufacturer's instructions and applied to Illumina HiSeq X Ten system for 150 nt paired-end sequencing.

RT-qPCR experiment

To determine the validity of the RNA‑seq result, RT-qPCR was performed for selected DEGs and DElncRs. Internal control was selected as the housekeeping gene GAPDH. Primers are designed in transcribed regions, which were presented in Table S1. The detailed protocol was fully described in a previously published study

²⁶

.

Process of RNA-seq raw data

The raw reads were trimmed of low-quality bases using a FASTX-Toolkit (v.0.0.13; http://hannonlab.cshl.edu/fastx_toolkit/). Then the clean reads were evaluated using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc). Quality control statistics, including raw reads, clean reads after filtering, and mapping percentages for each sample, are summarized in Supplementary Table S2.

Reads alignment and differentially expressed gene (DEG) analysis

Clean reads were aligned to the human GRch38 genome by tophat2

²⁷

allowing 4 mismatches. Uniquely mapped reads were ultimately used to calculate read number and fragments per kilobase of exon model per pillion mapped fragments (FPKM) for each gene. The expression levels of genes were evaluated using FPKM. The software edgeR

²⁸

,which is specifically used to analyze the differential expression of genes, was applied to screen the RNA-seq data for DEGs. The results were analyzed based on the fold change (FC ≥ 1.5 or ≤ 0.37) and p-value (p-value ≤ 0.01) to determine whether a gene was differentially expressed.

LncRNA prediction and direction identification

To systematically analyze the lncRNA expression pattern, we used a pipeline for lncRNAs identification similar as previously reported

²⁹

, which was constructed based on the cufflinks software

³⁰

.

Co-expression and functional enrichment analysis

To explore the regulatory mode between lncRNAs and their host mRNAs, we calculated the Pearson’s correlation coefficients (PCCs) between them and classified their relation into three class: positive correlated, negative correlated and non-correlated based on the PCCs value (|correlation | > 0.6, p-value < 0.05).

To sort out functional categories of gene sets, Gene Ontology (GO) terms and KEGG pathways were identified using KOBAS 2.0 server

³¹

. Hypergeometric test and Benjamini-Hochberg FDR controlling procedure were used to define the enrichment of each term.

Identification and target prediction of differentially expressed circRNAs

Circular RNAs were predicted from our RNA‑seq data using find_circ software (version 1.2) with default parameters

³²

. Differential expression analysis of circRNAs between paired groups was performed with the TCC package in R

³³

, employing a negative binomial model to estimate significance and false discovery rate (FDR). CircRNAs with an absolute fold change greater than 2 and an FDR below 0.05 were considered differentially expressed. The TCC analysis was configured with the following parameters: normalization method = "tmm", test method = "edger", iterations = 3, and floorPDEG = 0.05. To predict the miRNA targets of these differentially expressed circRNAs, we applied miRanda (v3.3a) with default settings. Candidate circRNA‑miRNA interactions were retained only when the miranda score exceeded 500, ensuring a stringent selection for subsequent ceRNA network construction.

miRNA target prediction

The potential miRNA targets were identified using miRanda (v3.3a) and TargetScan (v6.0) with default settings.All predicted target genes were evaluated by scoring system, and the sequences were considered to be miRNA targets only if miranda score was more than 500 or targetscan context + score was less than 0. To reduce the false identification rate, the gene expression profile of putative targets was required to be negatively correlated with the profile of miRNAs. For each paired samples, we identified the set of differentially expressed target genes and their corresponding miRNAs significantly changing in the opposite way. The mRNA-miRNA pairs whose negatively correlated relationship was detected in at least two paired samples were retained further to plot network by Cytoscape.

Phenotyping and differentiation of hUC-MSCs

Human umbilical cord MSCs (hUC-MSCs, passage 2, cat. no. SNP-H198, SUNNCEL) were purchased and cultured in DMEM-F12 (BI) containing 10% fetal bovine serum (FBS, Gibco) at 37°C in a 5% CO₂ incubator. HLA-DR, CD34, CD45, CD73, CD90, and HLA-ABC (Biolegend) were analyzed by FACS.

Adipogenic and osteogenic differentiation. hUC-MSCs were seeded at a density of 1×105 cells/well in a six-well plate. When the cells reached 100% confluence, the corresponding osteogenic induction solution (hUC-MSC osteogenic differentiation complete medium, HUXUC-90021, Cyagen, China) and adipogenic induction solution (hUC-MSC adipogenic differentiation complete medium, HUXUC-90031, Cyagen, China) were added. The medium was changed according to the operating procedure. After 24 days of induction, Alizarin red staining and Oil Red O staining were carried out for osteogenesis and adipogenesis respectively, then the cells were photographed and observed under a microscope (Olympu，Japan).

Preparation of IFN-γ pre-conditioned MSC-derived culture medium

Passage 3 hUC-MSCs were cultured in six-well plate until subconfluence to obtain culture medium (MSC-CM). Human recombinant IFN-γ (500 U/ml, R & D Systems, Minneapolis, MN, USA) was added in the plate for 12 hours, then culture medium was depleted, and cells were washed twice with HBSS to remove serum and IFN-γ before incubated in 6mL of α-MEM for another 48h. The resultant supernatants were centrifuged at 300 × g for 5 min to remove cells, then at 2,000 × g for 10 min to eliminate cell debris, and finally filtered through a 0.22 µm pore‑size filter to ensure sterility and remove any remaining particulate matter. The clarified medium was collected and used immediately. To minimize variation between preparations, we seeded the same number of cells, used identical culture volumes (6 mL per well), and collected the supernatants after the same incubation period (48 h) 34. The collected medium was not further concentrated or normalized by protein concentration; it was used directly as the conditioned medium. Three independent batches were prepared, and each experiment was repeated with at least two batches to confirm reproducibility.

Preclinical Model Bioinformatics and Transcriptomic Analysis

To investigate the preclinical cellular mechanisms of SARS-CoV-2, secondary bioinformatics analyses were conducted using published datasets. Raw read counts from GSE147507 were utilized to evaluate PTGS2 expression in Calu-3 cells subjected to SARS-CoV-2 infection (MOI=2) versus Mock controls at 24 hours. For clinical tissue-level validation, raw count data from GSE150316 were extracted. Raw count matrices for both datasets were processed and normalized using log2(Counts+1). Statistical significance for differential gene expression between infected/pathological groups and respective baseline controls was determined using the Wilcoxon rank-sum test.

Multi-omics COVID-19 Cohort Analysis and Functional Intersect

Multi-omics data from the Long COVID discovery cohort (GSE150316, GSE147507) were evaluated. Normalized transcriptomic expression values of PTGS2 were extracted and stratified to compare Long COVID subjects against healthy and non-Long-COVID (NLC) post-COVID controls (Panel A). To evaluate the relationship between continuous pathological drivers and clinical decline, Spearman rank correlation tests were applied between PTGS2 expression levels and derived clinical pulmonary function impairment scores (Panel B).Directional network dependencies and multi-omics module intersections were derived strictly from patient-matched data matrices. To map shared biological processes, differentially expressed genes (DEGs) from our in-vitro PTGS2-regulatory network (incorporating IL6, CXCL8, STAT1, PTGS2, ICAM1, and VCAM1) were intersected with clinical DEGs from the Long COVID cohort. Intersecting genes were subjected to Gene Ontology (GO) enrichment analysis focusing on biological processes (Benjamini-Hochberg FDR-adjusted P < 0.05$). All statistical comparisons were two-sided, with multiple testing controlled via Benjamini-Hochberg false discovery rate adjustment where applicable.

Statistical analysis

Principal component analysis (PCA) analysis was performed by R package factoextra (https://​cloud.r-project.org/​package=factoextra) to show the clustering of samples with the first two components. After normalizing the reads by TPM (Tags Per Million) of each gene in samples, in house-script (sogen) was used for visualization of next-generation sequence data and genomic annotations. The pheatmap package (https://cran.r-project.org/web/packages/pheatmap/index.html) in R was used to perform the clustering based on Euclidean distance. Student’s t-test was used for comparisons between two groups.

All quantitative data are presented as the mean ± standard deviation (SD) unless otherwise specified. Prior to parametric comparisons, the normality of data distribution was verified using the Shapiro-Wilk test, and homogeneity of variance was assessed utilizing Levene's test. Data conforming to a normal distribution and equal variance were analyzed using appropriate parametric tests (e.g., unpaired Student's t-test for two groups or one-way ANOVA for three or more groups, followed by Tukey's post-hoc test). For non-parametric data, corresponding non-parametric methods (e.g., the Mann-Whitney U test or Wilcoxon rank-sum test) were applied. Statistical analyses were performed using GraphPad Prism 8. A value of P < 0.05 was considered statistically significant (P < 0.05, *P < 0.01, **P < 0.001; ns, not significant).

SARS-CoV-2 spike protein inhibits apoptosis of Calu-3 cells

Spike glycoprotein (S) is the envelope of coronavirus that mediates host cell infection and a primary determinant of cell tropism and pathogenesis. To systematically explore the signaling functions of the SARS-CoV-2 spike protein, we utilized a commercially purified recombinant RBD fragment of the S1 subunit to treat human lung adenocarcinoma epithelial (Calu-3) cells. Transcriptome sequencing experiment (RNA-seq) was performed to investigate the expressed RNA level changes by spike protein(Fig. 1A). We first examined the effect of S protein treatment on the survival of Calu-3 cells in vitro. As seen in Fig. 1B, with the increase of the concentration of S protein, the apoptotic level of the Calu-3 cells was markedly decreased. Accordingly, treatment at a higher dose (50 to 500 nM) significantly reduced the apoptosis, which showed a steady dose dependent manner (Fig. 1B and C). To some extent, it suggests that SARS-CoV-2 spike protein could have the ability to inhibit lung epithelial cell apoptosis. Extensive studies have shown that the apoptosis signaling pathway plays a critical role in the progression of SARS-CoV virus infection

¹⁴

. For example, a previous study demonstrated that SARS-CoV activates apoptosis of human primary T lymphocytes in vivo

³⁵

. In addition, adenovirus-mediated expression of spike protein of SARS-CoV also induced apoptosis in Vero E6 cells in vitro

³⁶

. However, SARS-CoV-2 is particularly enigmatic, with seemingly contradictory roles in the apoptosis, which was recently found that SARS-CoV-2 virus-like particles (trVLP) can suppressed cell apoptosis

³⁷

. To further confirm our experimental results, we re-performed the experiment using an additional condition. As actinomycin D (Act D) treatment promoted DNA damage and induced apoptosis in lung cancer cells

³⁸

. Therefore, Calu-3 cells were firstly infected with S protein for 48h, then Act-D was added for another 12h. As we expected, Act D clearly increased the Calu-3 cell apoptosis rate (Fig. 1D and E), however, S protein significantly inhibited the pro-apoptotic effect of Act D with a dose dependent manner. In addition, another batched experiment was also performed. Calu-3 cells were pretreated with Act D for 12h, then washed and cocultured with S protein (0, 10, and 100 nM) for 48h. Then, as before, flow cytometry was used for apoptosis detection. Undoubtedly, spike protein still suppressed the Act D-induced apoptosis (Fig. 1F). These results did reveal an interesting phenomenon that SARS-CoV-2 S protein plays an apoptosis-inhibiting role on the human epithelial cells even in the presence of actinomycin D stimulation. Whether the decreased apoptosis was accompanied by an increase in the proliferating rate, EdU proliferation assay was used. As shown in (Fig. S1 A and B), no statistical difference of proliferation rates was observed between control and 1nM, 10nM, 50nM (p > 0.05), whereas there was a statistical proliferation difference between control and 100 or 500nM group, converged an intrinsic apoptosis pathway may be involved.

Spike protein extensively regulates the steady transcriptome profile of Calu-3 cells

The underlying mechanism of inhibition of cellular apoptosis by the S protein is largely unknown. We then performed whole transcriptome sequencing experiment (RNA-seq) for the treated cells by removing ribosomal RNAs (rRNAs) strategy. Following the apoptosis results in Fig. 1B, we selected Calu-3 cells only (control group, SD_0nM) in comparison with 100 nM (SD_100nM) S protein treatment (S-RBD group) both under Act D (0.75 μM) condition. Three biological replicates for each group were prepared for RNA-seq. After filtering low quality reads and aligning reads to human GRCh38 genome sequences, we obtained the genomic region distribution and expression levels for each gene (fragments per kilobase per million, FPKM) by classifying the uniquely aligned reads onto human genome. Genomic region distribution revealed most of the aligned reads were located into the transcribed regions, including 5’UTR, CDS, intron, 3’UTR, and nc_exon (Fig. 2A). There was no difference between SD_0nM and SD_100nM groups (Fig. 2A). We then analyzed whether spike protein could alter global gene expression pattern. Principal component analysis (PCA) showed that SD_0nM samples were clearly separated from SD_100nM samples by the first component, which could explain 32.3% of the total variation (Fig. 2B). Following differentially expressed gene (DEG) analysis revealed 223 upregulated and 140 downregulated genes between SD_100nM vs. SD_0nM with criteria p-value < 0.01 and fold change (FC) > 1.5 (Fig. 2C). Hierarchical clustering heat map of the expression levels of DEGs showed up DEGs and down DEGs were clearly separated (Fig. 2D). We then performed functional enrichment analysis for up and down DEGs based on gene ontology (GO) annotation. For upregulated genes, we found they were enriched in inflammatory response, cell adhesion, GTP catabolic process, axon guidance, small GTPase mediated signal transduction, positive regulation of apoptotic process, carbohydrate metabolic process, signal transduction transmembrane transport, and apoptotic process (Fig. 2E). Several genes involved in apoptotic process were upregulated by S protein treatment, including RHOB, IL24, PTGS2, NGFR, UNC5C, RRAGC, and TMEM102. Meanwhile, down DEGs were enriched in metabolic process related pathways (Fig. 2E and F).

Validation of DEGs associated with apoptosis

To further validate the apoptosis‑related DEGs identified from the RNA‑seq data, we prioritized genes with high fold changes, significant p‑values, and documented involvement in inflammatory or apoptotic signaling, then selected several candidates for RT‑qPCR confirmation, including both up‑regulated and down‑regulated genes. Among the down DEGs, mitochondrially encoded genes were observed among the top significant DEGs sorted by p-value, including MT-ND1, MT-ND2 and MT-ND4. Previous studies have demonstrated MT-ND1, PTGS2 and GTPase play an important role in SARS-CoV-2 affection

^39,

^40,

⁴¹

. We finally selected four genes for validation, including RHOB, PTGS2, RRAGC and MT-ND1, which showed consistent and significant results with the RNA-seq data (Fig. 3A and B), suggesting the authenticity of RNA-seq result. After comparison, we then focused on the PTGS2. Network pharmacology and molecular docking have been revealed that PTGS2 was identified as a potential target in treating COVID-19, which was associated with the suppression of oxidative stress and cell apoptosis

⁴²

. Moreover, other study also confirmed that PTGS2 can inhibit cell apoptosis via the ROS-ERK1/2-NF-κB signaling axis, which enhanced the resistance of non-small cell lung cancer cells against the cisplatin

⁴³

.

MSC-CM reverses the inhibition of Calu-3 cell apoptosis by spike protein, correlating with PTGS2 expression changes

There was research that shows MSCs inhibit the apoptosis of alveolar epithelial cells caused by acute respiratory distress syndrome (ARDS) through CXCL12/CXCR4 axis

⁴⁴

. Considering the importance of immunomodulation in COVID-19, MSCs might modulate the apoptosis of SARS-CoV-2 infected cells. Accumulating research findings confirmed that MSCs pretreated with some inflammatory cytokines, such as IFN-γ, exhibit enhanced immunomodulatory functions

^21,

⁴⁵

. However, whether conditioned medium from IFN-γ-stimulated MSCs (MSC-CM) exerts beneficial effects in terms of SARS-CoV-2 infected cells remains unknown. Based on this knowledge, we initially investigated whether MSC-CM have the potential of reversing the anti-apoptotic potent of S protein on Calu-3 cells.

hUC-MSCs were isolated and identified using procedures described in our previous study. These cells had the ability to differentiate into osteoblasts and adipocytes (Fig. S2 A) after culture in osteogenic and adipogenic media followed by Alizarin red and Oil Red-O staining respectively. Meanwhile, the surface markers of CD73 (99.94%), CD105 (99.86%), and CD90 (99.99%) were positive (>95%), and hematopoietic markers such as CD45 (0.04%), HLA-DR (0.19%) and CD34 (0.12%) were present at low level (Fig. S2 B). MSC-CM was collected from passage 2 hUC-MSCs. Therefore, to determine whether MSC-CM has the ability to attenuate the anti-apoptotic role of S protein, 100 nM concentration of S protein were applied. As we expected, S protein markedly mitigated the apoptosis of Calu-3 cells in normal culture medium. If the cells were cultured in MSC-CM, the apoptotic frequency was significantly increased. More importantly, MSC-CM did diminish the anti-apoptotic effect of S protein (Fig. 4A and B). Additionally, using Act D pretreated Calu-3 cells, we also observed that MSC-CM has a powerful effect of reversing S protein-induced anti-apoptosis ability (Fig. 4C and D). In order to explore the possible mechanism whereby MSC-CM showed strong antagonistic action against S protein, we locked the candidate molecule to PTGS2. A similar result was observed that S protein significantly increased PTGS2 protein level while MSC-CM decreased the PTGS2 production (Fig. 4E and F). In summary, we here confirmed that conditioned medium from IFN-γ-stimulated MSCs strongly induced Calu-3 cell apoptosis, whereby S protein mediated anti-apoptosis effect was also suppressed by MSC-CM. The downregulation of S protein-induced PTGS2 expression by MSC-CM correlates with the restoration of apoptosis, positioning PTGS2 as a candidate molecule for future mechanistic dissection.

Long noncoding RNAs (lncRNAs) were dysregulated by spike (S) protein

Long noncoding RNAs (lncRNAs) are emerging vital regulators in most aspect of cellular functions

⁴⁶

. We then investigated the expression profile of lncRNAs by S protein infection. After predicting novel lncRNAs and identifying expressed lncRNAs, we found the two groups shared a large number of detected lncRNAs, including 1810 co-detected lncRNAs. We also detected 217 and 172 lncRNAs were specifically expressed in SD_100nM and SD_0nM samples, respectively(Fig. 5A). Differentially expressed lncRNAs (DElncRs) analysis revealed 125 upregulated and 58 downregulated lncRNAs in SD_100nM group(Fig. 5B). We observed that both more lncRNAs and mRNAs were detected in S protein treated samples, implying that S protein may stabilize the RNA levels. Hierarchical clustering heat map of the expression levels of DElncRs also showed up lncRNAs and down lncRNAs were clearly separated(Fig. 5C), although the expression pattern among replicates was less consistent than that of mRNAs because of the lower expression levels of lncRNAs. LncRNAs have the ability to regulate gene expression in cis- or trans-manner

⁴⁷

, so we performed lncRNA-mRNA co-expression analysis to predict the potential targets of DElncRs. By calculating the correlated mRNA numbers and statistical significance of co-expressed lncRNAs, we detected several lncRNAs exhibiting more co-expressed mRNAs and significant false discovery rate (FDR), including RP13-14G15.4, AC004988.1, AC003092.1, GAPLINC, RP11-265F19.1, and RP11-459E10.1 (Fig. 5D). Functional analysis of the mRNA that were co-expressed with DElncRs revealed they were highly enriched in inflammatory response, respiratory electron transport chain, cell adhesion, GTP catabolic process, and cellular metabolic process pathways (Fig. 5E). We selected five lncRNAs that were upregulated in SD_100nM samples and highly correlated with mRNAs to show their expression pattern (Fig. 5F). To make a more definite expression changes of these lncRNAs after spike protein treatment, we selected these five lncRNAs for RT-qPCR experiment. The result showed consistent expression changes with RNA-seq result (Fig. 5G).

Circular RNAs acts as competitive endogenous RNAs to participate in Calu-3 cellular functions

Circular RNAs (circRNAs) are emerging as a new class of non-coding RNAs with important molecular functions

⁴⁸

. Recent study demonstrated the important roles of circRNAs during Middle East respiratory syndrome coronavirus (MERS-CoV) infection

⁴⁹

. We used rRNA depletion method to prepare RNA-seq library, enabling us to investigate the circRNA expression profile in our samples. By predicting circRNAs using FindCirc software

⁵⁰

, a total of 5281 circRNAs were obtained and there was no obvious difference of circRNA numbers for the six samples (Fig. 6A). We then used normalized back-splicing reads number to perform differential expression analysis of circRNAs, and finally obtained 39 upregulated and 38 downregulated circRNAs by spike protein treatment (Fig. 6B). Although the differentially expressed circRNAs (DECircRs) only account for 1% of the detected circRNAs, hierarchical clustering heat map showed a high consistency between biological replicates (Fig. 6C). Most of the circRNAs are generated from exon regions of genes by abnormal splicing, which perhaps could affect the expression of their host genes. We then performed functional analysis for the host genes of DECircRs. For upregulated circRNAs, their host genes were enriched in metabolic pathway, lysosome, and several other glycosphingolipid related pathways (Fig. 6D). Host genes of downregulated circRNAs were enriched in ubiquitin mediated proteolysis, hippo signaling pathway, propanoate metabolism and RNA polymerase pathways (Fig. 6E). Two genes involved in lysosome were HEXB and ATP6V1H, which were also upregulated by spike protein. On canonical function of circRNAs is to function as competitive endogenous RNA by sponging miRNAs in cells

⁴⁸

. We then constructed a competitive endogenous RNA network for these dysregulated circRNAs by predicting their binding miRNAs and associated mRNAs. Eight circRNAs (red points), 34 mRNAs, and 14 miRNAs participated in this network (Fig. 6F). One circRNA originated from FAM178A could be served as potential target of eight miRNAs (Fig. 6F). Functional analysis of the involved mRNA genes showed they were enriched in metabolism related pathways (Fig. 6G).

PTGS2 expression across preclinical models and fatal COVID-19 tissues

To validate the clinical relevance of PTGS2-mediated apoptosis inhibition, we interrogated public datasets (GSE147507, GSE150316). Transcriptomic analysis confirmed that PTGS2 is significantly upregulated in Calu-3 cells upon SARS-CoV-2 infection (Figure 7A). This upregulation was consistently observed across various pulmonary cell lines, including NHBE and A549-ACE2, highlighting a robust cellular response (Figure 7B). Furthermore, examining autopsied lung tissues from COVID-19 fatalities revealed markedly higher PTGS2 expression compared to healthy controls (Figure 7C). These findings suggest that the S protein-induced PTGS2 axis is a conserved pathological feature in infected lung tissue, reinforcing its role as a potential driver of viral persistence and tissue-level apoptosis resistance.

PTGS2 as a persistent inflammatory node in Long COVID

Expanding our clinical validation to the post-acute sequelae of COVID-19, we analyzed a multi-omics Long COVID cohort. PTGS2 expression remained persistently elevated in Long COVID patients relative to controls (Figure 8A). Spearman correlation analysis demonstrated a trend where increased PTGS2 levels inversely correlate with pulmonary function markers, suggesting its involvement in impaired tissue repair (Figure 8B). Crucially, functional enrichment of the intersection between our in vitro PTGS2-regulatory network and Long COVID differentially expressed genes highlighted pathways encompassing cytokine signaling, inflammatory response, and prostaglandin biosynthesis (Figure 8C). Collectively, these clinical signatures suggest that the PTGS2-driven inflammatory microenvironment, reversible by MSC-CM in our model, remains a critical therapeutic target for mitigating Long COVID-related chronic immune dysfunction.