All the experiments related to live SARS-CoV-2 viruses were approved by the Ethics Committee of the Zhejiang Provincial Center for Disease Control and Prevention (ZJCDC) in China, and carried out in the biosafety level 3 (BSL-3) laboratory of ZJCDC.

The SARS-CoV-2 clinical strains F8 and 8X were purified from the pharyngeal swab of a male COVID-19 patient in Hangzhou as mentioned previously (Yao et al., 2020). Calu-3 cells were obtained from the National Collection of Authenticated Cell Cultures and cultured in DMEM (Gibco, United States) with 10% fetal bovine serum (FBS, Every Green, China) at 37°C in a 5% CO₂ incubator. Cells were seeded into 6-well plates at a density of 1*10⁶ cells/well, and were infected by the F8 and 8X strains at a multiplicity of infection (MOI) of 2. Phosphate-buffered saline (PBS) was used as a negative control. Each group had two replicates. One hour post-adsorption at 37°C, the viral inocula were discarded, and the cells were maintained in the virus growth medium (DMEM containing 3% FBS) after washing twice with PBS. Two days post-infection, cells were collected, and the total RNA was extracted by using an RNeasy Plus Mini Kit.

The RNA concentrations and quality of each sample were measured using Nanodrop One. The RNA integrity was detected by Agilent 2,100. Total RNA was used by depleting ribosomal RNA according to the manuscript of the Ribo-Zero rRNA Removal Kit. The rRNA-depleted RNA was fragmented, and the cDNA library was constructed using the TruSeq RNA sample Prep Kit (Illumina, San Diego, CA, USA). The sequencing libraries were sequenced on an Illumina NovaseqTM 6,000 platform according to the manufacturer’s instructions. The sequence data generated from this project has been deposited in NCBI under SRA submission PRJNA909976.

Transcript assembly: First, Cutadapt (Martin, 2011) was used to remove the reads that contained adaptor contamination, low-quality bases and undetermined bases. Then, sequence quality was verified using FastQC (Babraham Bioinformatics, 2022). We used Bowtie2 (Langmead and Salzberg, 2012) and Hisat2 (Kim et al., 2015) to map reads to the genome of Homo sapiens. The mapped reads of each sample were assembled using StringTie (Pertea et al., 2015). Then, all transcriptomes from three samples were merged to reconstruct a comprehensive transcriptome using Perl scripts. After the final transcriptome was generated, StringTie and edgeR (Robinson et al., 2010) were used to estimate the expression levels of all transcripts.

LncRNA identification: First, transcripts that overlapped with known mRNAs and transcripts shorter than 200 bp were discarded. Then, we utilized CPC (Kong et al., 2007) and CNCI (Sun et al., 2013) to predict transcripts with coding potential. All transcripts with CPC scores < −1 and CNCI scores <0 were removed. The remaining transcripts were considered as lncRNAs.

Differential expression analysis of mRNAs and lncRNAs: StringTie was used to determine the expression level for mRNAs and lncRNAs by calculating FPKM (Trapnell et al., 2010). The differentially expressed mRNAs and lncRNAs were selected with log2 (fold change) >1 or log2 (fold change) < −1 and with statistical significance (I value < 0.05) by R package – edgeR (Robinson et al., 2010).

Target gene prediction and functional analysis of lncRNAs: To explore the function of lncRNAs, we predicted the cis-target genes of lncRNAs. LncRNAs may play a cis role by acting on neighboring target genes. In this study, coding genes in 100,000 upstream and downstream were selected by python script. Furthermore, trans-regulation analysis is a genome-wide search for well-associated target genes. Correlation analysis was performed between lncRNAs and the corresponding gene set. Associations with Pearson correlation coefficients greater than 0.4 (p < 0.05) were presumed to have targeted regulatory effects. Then, functional analysis of the target genes for lncRNAs was performed by using the BLAST2GO (Conesa et al., 2005).

Gene Ontology (GO) enrichment analysis of differentially expressed genes or lncRNA target genes was conducted with respect to biological process, molecular function, and cellular component. Kyoto Encyclopedia of Genes and Genomes (KEGG) was used to perform pathway enrichment analysis. The R package clusterprofiler was used to perform the detailed enrichment analysis described above (Wu et al., 2021).

For each lncRNA, the Pearson correlation coefficient of its expression value with that of each protein-coding gene was calculated. Under the conditions of an absolute value of the Pearson correlation coefficient > 0.90 and p < 0.001, the interaction network of the differentially expressed lncRNAs and protein-coding gene co-expression pairs was then constructed using Cytoscape (Shannon et al., 2003).

When constructing the competing endogenous RNA (ceRNA) network, lncRNAs were connected to the differentially expressed (DE) human miRNAs if they were predicted to interact with each other, and the upregulated (downregulated) miRNAs were connected to downregulated (upregulated) mRNAs if the former targeted the latter based on the database miRTarBase (Huang et al., 2020). The lncRNA-miRNA-mRNA network was visualized using Cytoscape (version 3.9.1). The co-expression network was built as follows: A custom database from a combination of public databases, starBase (version 2.0) and miRcode (version 11), was used to predict the interaction between known lncRNAs and miRNAs. The interaction of miRNA with novel lncRNA, circRNA, and mRNA was performed by using TargetScan (version 8.0) and miRanda (version 22.1), using default parameters. Pairs that appeared in both results were considered as prospective interactions. Subsequently, correlation analysis was performed between ceRNAs (circRNA, lncRNA, mRNA) and miRNAs. Based on the correlation between them, a hypergeometric distribution analysis was performed under the threshold of probability less than 0.05. The construction of the ceRNA network followed the following principles: (1) Pearson correlation coefficient within ceRNA pairs or between ceRNAs and miRNAs should be under the threshold of an absolute value of 0.4 (p < 0.05); and (2) At least one miRNA must satisfy the hypergeometric distribution test between two ceRNAs.

All the statistical analyzes in this study were conducted in R (version 4.2.1). The Wilcoxon rank-sum test was used to compare the sample means between different groups, and was conducted with the wilcox.test function. Significance was expressed as a p value < 0.05.

To identify different transcripts in SARS-CoV-2 infected cells, the transcriptomes of the Calu-3 cells infected with F8, 8X, or without SARS-CoV-2 infection as controls were detected using high-throughput RNA sequencing. Robust and reproducible data were obtained from all samples. After quality filtering, clean reads were mapped to the human reference genome (hg38) using HiSat2, and were assembled with StringTie. Then, coverage analysis was performed on these clean reads on different annotated gene types. The distribution of each type of gene was counted according to the expression level (Figure 1). In total, eight categories of RNA were identified according to database annotation of those transcripts, in which the protein-coding genes were highly represented (72.21% in F8 and 74.43% in 8X, respectively).

DEGs were identified from the control vs. SARS-CoV-2 infection groups and used for functional annotation. All DEGs were demonstrated in volcano plots (Figures 2A–C). In the F8 vs. 8X group, there were 149 upregulated DEGs and 41 downregulated DEGs. The comparison between the F8 vs. control or 8X vs. control groups revealed that a total of 1,372 or 1,275 genes were significantly regulated after viral infection. A Venn diagram was plotted to identify the common or unique DEGs in human lung cells with or without SARS-CoV-2 infection (Figure 2D and Support file 2). A total of 14 common DEGs were identified. Given the exposure of cells to F8, the number of DEGs was smaller (345 genes vs. 394 genes) in comparison of the 8X vs. control group. The top 10 DEGs, including the E-selectin (SELE), colony stimulating factor 2 (CSF2) and pentraxin 3 (PTX3), were shown in the heatmap (Figure 2E).

GO enrichment of the DEGs was divided into three types: biological process, cellular component, and molecular function (Figure 3A). Each top ten terms of the three types, such as cytokine activity and membrane raft, were displayed. The top 20 enriched pathways were shown through KEGG functional analysis of the DEGs (Figure 3B). Significant enrichment in the TNF signaling pathway, viral protein interaction with cytokine and cytokine receptor, and cytokine-cytokine receptor interaction were found.

The assembled transcripts were used to identify lncRNAs using the pipeline as described in the Methods. In total, 110 annotated and novel lncRNAs were identified. It has been reported that lncRNAs, in comparison with protein-coding genes, usually share some common genomic features with their sequences. However, they are generally shorter in length, have fewer but longer exons, and have lower expression levels and lower evolutionary sequence conservation (Liu et al., 2019). To further determine the characteristics of the lncRNAs identified in the present study, we compared transcript length, exon number and expression levels of protein-coding genes and lncRNAs. The results revealed that fewer exons and lower expression levels were also found in the lncRNAs, which was consistent with the reported lncRNAs (Figure 4).

We next performed GO and KEGG pathway analyzes for the target genes of lncRNAs between 8X and F8. The top 20 GO terms and KEGG pathways were reported in Figure 5. The GO term analysis divided differentially expressed lncRNAs into the same three types as DEGs. Biological processes such as response to chemokine and protein localization to cytoskeleton, the molecular functions such as cytokine activity and translation regulator activity, and the molecular functions such as membrane raft and transcription regulator complex, were all enriched. KEGG enrichment analysis revealed that pathways related to the immune system, such as viral protein interaction with cytokine and cytokine receptor, chemokine signaling pathway and Toll-like receptor signaling pathways, were preferentially targeted. Therefore, the results of GO and KEGG pathway enrichment analyzes revealed that lncRNAs may act in cis or trans to participate in the regulation of the expression of multiple important genes in different processes, including the immune response and protein localization.

The functional association between regulatory lncRNAs and protein-coding gene transcripts can be determined by performing expression correlation analysis. To further investigate the potential mechanism of the SARS-CoV-2 associated lncRNAs, the DE lncRNAs and their predicted target DE protein-coding genes were investigated by delineating lncRNA-protein-coding gene functional interactions. Here we identified 325 pairs of DE lncRNA-DE protein-coding genes between 8X and F8 (Support file 1), containing 22 lncRNAs and 77 protein-coding genes (p < 0.01). The network of co-expressed lncRNA-protein-coding gene pairs based on a threshold of Pearson’s correlation coefficient of 0.998. Next, the KEGG pathway analysis of co-expressed genes revealed the top 10 pathways, including the TNF signaling pathway, viral protein interaction with cytokine and cytokine receptor, and cytokine-cytokine receptor interaction, which were all significantly enriched (Figure 6A).

The network of lncRNAs with viral protein interaction with cytokine and cytokine receptor (Figure 6B), and TNF signaling pathway-related genes (Figure 6C) were analyzed. Several common cytokines, such as IL-6, CXCL2 and CXCL10, were found in all these pathways. Moreover, the CSF2, one of the top ten DEGs between 8X and F8 groups, was also involved in both the viral protein interaction with cytokine and cytokine receptor pathway and the TNF signaling pathway.

Based on the lncRNA-mRNA co-expression relationship and the regulation relationship of DE-miRNA-DE-mRNA and DE-miRNA-DE-lncRNA, the lncRNAs and mRNAs that were significantly differentially expressed and regulated by the same miRNA were screened. In total, 89 lncRNA-miRNA-mRNA interactions were finally obtained. Furthermore, according to the lncRNA-miRNA-mRNA and circRNA-miRNA-mRNA networks, differentially expressed circRNAs, lncRNAs, and mRNAs that were regulated by the same miRNA were further screened. Finally, 101 interaction pairs were obtained (Figure 7). Among these interactions, there were 2 upregulated lncRNAs and 2 upregulated mRNAs, and 1 upregulated miRNA in Calu-3 cells infected with the F8 variant compared with the 8X strain. There were 12 lncRNAs, 25 upregulated and 11 downregulated mRNAs, and 1 upregulated and 7 downregulated miRNAs in Calu-3 cells infected with F8 or X8 compared with the control.