Comparative transcriptomic analyzes of human lung epithelial cells infected with wild-type SARS-CoV-2 and its variant with a 12-bp missing in the E gene

The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is a novel coronavirus that caused a global outbreak of coronavirus disease 2019 (COVID-19) pandemic. To elucidate the mechanism of SARS-CoV-2 replication and immunogenicity, we performed a comparative transcriptome profile of mRNA and long non-coding RNAs (lncRNAs) in human lung epithelial cells infected with the SARS-CoV-2 wild-type strain (8X) and the variant with a 12-bp deletion in the E gene (F8). In total, 3,966 differentially expressed genes (DEGs) and 110 differentially expressed lncRNA (DE-lncRNA) candidates were identified. Of these, 94 DEGs and 32 DE-lncRNAs were found between samples infected with F8 and 8X. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyzes revealed that pathways such as the TNF signaling pathway and viral protein interaction with cytokine and cytokine receptor were involved. Furthermore, we constructed a lncRNA-protein-coding gene co-expression interaction network. The KEGG analysis of the co-expressed genes showed that these differentially expressed lncRNAs were enriched in pathways related to the immune response, which might explain the different replication and immunogenicity properties of the 8X and F8 strains. These results provide a useful resource for studying the pathogenesis of SARS-CoV-2 variants.

Introduction

The novel coronavirus disease 2019 (COVID-19) is an acute respiratory infectious disease caused by a severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) infection, which was first identified in China (Zhu et al., 2020). All the human beings are generally susceptible to the SARS-CoV-2 virus. By the 4th of December (Kim et al., 2022), it had led to more than 640 million patients and more than 6.6 million deaths worldwide (WHO, 2022). The SARS-CoV-2 virus is a single-stranded RNA virus and can mutate easily. Since the outbreak of COVID-19 in the late of 2019, the SARS-CoV-2 virus has continuously evolved with many variants emerging across the world, such as Alpha, Beta, Delta and Omicron variants (Cao et al., 2022; Iketani et al., 2022). In particular, the Omicron variants BA.4 and BA.5, containing the prominent immune escape characteristics, have spread all over the world and become the dominant strains currently (Cao et al., 2022). In China, a large part of the clustered outbreaks and sporadic cases were caused by the BA.5 infection (Feng et al., 2022; Jiang et al., 2022), which placed great pressure on COVID-19 prevention strategies. Therefore, it is necessary to pay more attention to SARS-CoV-2 variants.

The envelope (E) protein, located at the viral envelope, is the smallest one of four major structural proteins of SARS-CoV-2 virus. It consists of only 76 amino acids, but plays important roles in the viral life cycle, such as viral assembly, budding and so on (Castano-Rodriguez et al., 2018). A SARS-CoV strain lacking the E protein is attenuated in vivo (Regla-Nava et al., 2015). In our previous study, we isolated both the E gene wild-type SARS-CoV-2 strain 8X and the E gene mutant strain F8 from the same specimen (Sun et al., 2020). Although no significant difference in the viral titer and infectivity was found in the E gene mutant SARS-CoV-2 strain F8 and the E gene wild-type strain 8X, F8, which contains a 12-bp deletion in the E gene, could produce a higher S protein content and induce a quicker humoral immune response than 8X. The inactivated SARS-CoV-2 vaccine produced from the F8 strain could trigger higher levels of the IgG titer and neutralizing antibody titer than those from the 8X strain at 1 and 3 weeks post vaccination, respectively. It seemed that the E gene mutation could influence the replication and immunogenicity of SARS-CoV-2. However, the mechanism has not been elucidated yet.

In this study, the whole-transcriptome sequencing was performed based on the F8-and 8X-infected human lung epithelial (Calu-3) cells. Subsequently, differentially expressed genes (DEGs) and lncRNAs of the F8 and 8X groups were analyzed, followed by functional interaction prediction analysis. Although large amounts of data have proven that several lncRNAs are involved in different kinds of viral infections, the underlying mechanisms by which they act are still largely unknown. Based on the whole-transcriptome sequencing analysis, our results may provide novel insights into the molecular basis of SARS-CoV-2 infection.

Materials and methods

Ethics statement

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.

Virus and cells

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.

RNA extraction and next-generation sequencing

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.

RNA-Seq data analysis

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).

Go and KEGG enrichment analysis

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).

Construction of the LncRNA-protein-coding gene co-expression network and competing endogenous RNA (ceRNA) network

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.

Statistical analysis

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.

Results

Transcriptome profiles of SARS-CoV-2 infected cells

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).

FIGURE 1

Figure 1. RNA categories of quantified transcripts in RNA-seq. The percentages of each RNA type are shown in F8/8X infected or mock-infected Calu-3 cells.

Identification of differentially expressed genes

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).

FIGURE 2

Figure 2. Volcano plots showing the DEGs between 8X and control (A), F8 and control (B), and the 8X and F8 (C). (D) Venn diagram of the identified DEGs shared among the above three groups. (E) Heatmap of the top 10 DEGs between 8X and F8 groups.

Go/KEGG analysis of DEGs

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.

FIGURE 3

Figure 3. Dot plots displaying enriched GO terms (A) or KEGG pathways (B) of identified DEGs. Dot size represents the count of the enriched gene within each category, and the x-axis represents the gene ratio (rich factor). The GO terms/KEGG pathways are arranged on the basis of the value of the gene ratio, not on their adjusted p value, for easier visualization purposes.

Identification of lncRNAs and Go/KEGG analysis of lncRNA target genes

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).

FIGURE 4

Figure 4. Comparison of the expression level (A) and exon number (B) of lncRNAs and protein-coding genes.

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.

FIGURE 5

Figure 5. Functional enrichment analysis of identified lncRNAs between 8X and F8. Representative over-represented KEGG pathway (top) and GO terms (bottom) of gene-expression clusters. BP, biological process; MF, molecular function; CC, cellular component.

Expression correlation analysis

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).

FIGURE 6

Figure 6. Co-expression analysis of lncRNAs and protein-coding genes. (A) The top 10 over-represented KEGG pathways of co-expressed genes. The network of lncRNAs with (B) viral protein interaction with cytokine and cytokine receptor, and (C) TNF signaling pathway-related genes (p < 0.01). The solid line represents a positive correlation, while the dotted line represents a negative correlation.

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.

CeRNAs network analysis

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.

FIGURE 7

Figure 7. The competing endogenous RNA (ceRNA) network. mRNAs, circRNAs, host lncRNAs and miRNAs in the networks are represented as squares, circles, triangles and arrows, respectively.

Discussion

Since the outbreak of COVID-19, a large number of studies have focused on the spike protein of SARS-CoV-2. However, few studies have focused on the E gene. Both the E gene and S gene are the structural proteins of SARS-CoV-2. E gene is important for the virus maturation and replication (Regla-Nava et al., 2015). Studies related to the E gene were mainly about its use as a target for the virus detection (Kim et al., 2022; Valadan et al., 2022). It seems that the E gene of SARS-CoV-2 is conserved. However, in our previous study (Sun et al., 2020), we found that a 12-bp deletion in the E gene of the F8 variant was more likely generated as a result of viral adaptation to Vero cells. The E gene of SARS-CoV-2 might have several mutations. Since the E gene mutant strain F8 and E gene wild-type strain 8X were both from the same specimen, these two strains had similar genetic backgrounds. It would be interesting to use them to investigate the function of the E gene.

In this study, by applying transcriptome sequencing, we found 3,966 DEGs and 110 differentially expressed lncRNAs in Calu-3 cells infected with F8 compared with 8X. Direct functionly enrichment analysis of the DEGs showed that these genes were mainly involved in cytokines and chemokines, which are related to the host immune response. What’s more, the co-expression analysis of differentially expressed lncRNAs and DEGs also showed the enrichment of the pathways involved in immune response, which may be critical for viral maturation, such as the TNF signaling pathway, and viral protein interaction with cytokine and cytokine receptor. A previous study showed that lncRNAs play a key role during viral infection (Liu et al., 2019). Additionally, some lncRNAs and miRNAs, such as hsa-miR-335-3p, hsa-miR-92a-1-5p, and hsa-miR-23a-5p, were identified as hub genes in the ceRNA network. Hsa-miR-92a-1-5p was found to enhance the interferon expression by targeting the suppressors of cytokine signaling 5 to inhibit feline panleukopenia virus replication in host cells (Liang et al., 2022). All the results strongly suggest that the genes involved in the immune response, especially the cytokines and chemokines, may play fundamental roles in the pathogenesis of different SARS-CoV-2 variants. Further study might be carried out in the SARS-CoV-2 E gene variant-infected animal models.

Through analysis of DEGs between the F8 and 8X groups, 85 upregulated genes were found. The E-selectin, CSF2 and PTX3 were the top 3 of the most upregulated genes. The E-selectin is important for leukocyte accumulation in inflammatory responses and PTX3 can amplify the immune response (Aslan et al., 2012; Ketter et al., 2016). In our previous research, a different immunogenicity of the inactivated COVID-19 vaccine produced by the F8 strain was found compared with that produced by the 8X strain (Sun et al., 2020). F8 vaccine might induce a higher expression level of E-selectin and PTX3 than 8X vaccine to promote a quicker humoral immune response in mice. The TNF-α signaling pathway is also important for the immune response, and is one of the most upregulated pathways in SARS-CoV-2 infected A549-hACE2 cells (Sun et al., 2021). The release of TNF-α could trigger several cell adhesion molecules, such as selectins and VCAM-1, to induce inflammation in vivo (Kong et al., 2018), which is a protective biological response for eliminating viruses. In our study, the TNF signaling pathway was also enriched significantly through the KEGG functional analysis of DEGs, the co-expression correction analysis of regulatory lncRNAs, and protein-coding gene transcripts. The expression levels of TNF in the F8 group were 24.3-fold higher than those in the 8X group. It seems that the TNF signaling pathway is also related to the pathogenesis of different SARS-CoV-2 variants.

Cytokines and chemokines are important for the immune response (Berantini et al., 2022). Among the viral protein interaction pathway, the immune system pathway and the TNF signaling pathway, cytokines such as IL-6 and CSF2, and the chemokines such as CXCL2 and CXCL10, were all involved. CXCL2 was also found to be one of the common DEGs in SARS-CoV-2 infected Calu-3, hCM and A549-hACE2 (both low and high viral loads) cells (Sun et al., 2021). IL-6 might be a reliable indicator for the COVID-19 severity, and a higher level of IL-6 concentration was found in severe COVID-19 patients than in the moderate or mild groups (Bergantini et al., 2022). The CSF2 protein was inhibited to reduce the production of inflammatory factors and chemotaxis of inflammatory cells (Xiong et al., 2022). CXCL10, also known as interferon gamma-induced protein 10 (IP-10), is a pro-inflammatory chemokine (Coperchini et al., 2020). Higher expression levels of IL-6, CSF2 and CXCL10 were also found in the F8 group than in the 8X group, indicating a more severe inflammatory response in the immunized mice of the F8 group.

In conclusion, our study systematically characterized transcriptome profiles during SARS-CoV-2 infection and provides a comprehensive genome-wide resource for identifying and functionally analyzing the differentially expressed genes and non-coding RNAs. Immune response signaling, such as upregulation of IL-6, CSF2 and CXCL10 cytokines and a higher level of E-selectin and PTX3, as well as the TNF-α signaling pathway, might be the reason to explain the different pathogenesis caused by the E gene mutant strain F8 and E gene wild-type strain 8X. Further functional validations are needed to delineate the exact mechanistic details, and an animal model might be used to study the pathogenesis of both the E gene mutant and wild-type strains. These results will be useful for better understanding the pathogenesis of SARS-CoV-2 variants, and designing better preventive and therapeutic measures against viral infection.

Data availability statement

The data presented in the study are deposited in the NCBI repository, accession number: PRJNA909976.

Ethics statement

The studies involving human participants were reviewed and approved by the Ethics Committee of the Zhejiang Provincial Center for Disease Control and Prevention (ZJCDC) in China. The patients/participants provided their written informed consent to participate in this study.

Author contributions

Y-SS, HS, P-PY, and Y-DL conceived and designed the experiments. Y-SS, H-PZ, FX, H-JL, and AT performed the experiments. G-LL, HS, and J-MJ analyzed the data. Y-SS, J-MJ, and Y-DL drafted the manuscript. All authors read and approved the final manuscript.

Funding

This work was supported by the Basic Public Welfare Research Project of Zhejiang Province (LGF22C080002, LGF22H260024, LDT23H19013H19), the Key R&D Program of Zhejiang Province (2021C03197, 2021C03200), and the Zhejiang Provincial Foundation for Scientific Research in Medicine and Health (WKJ-ZJ-2220).

Acknowledgments

We thank Jian-Hua Li from Zhejiang Provincial Center for Disease Control and Prevention for the technical assistance at the BSL-3 laboratory.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2022.1079764/full#supplementary-material

References

Aslan, M. K., Boybeyi, O., Soyer, T., Senyücel, M. F., Ayva, S., Kısa, U., et al. (2012). Evaluation of omental inflammatory response with P−/E-selectin levels and histopathologic findings in experimental model. J. Pediatr. Surg. 47, 2050–2054. doi: 10.1016/j.jpedsurg.2012.06.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Babraham Bioinformatics. (2022). FastQC. Available at: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/ (Accessed July 20, 2022).

Google Scholar

Bergantini, L., d'Alessandro, M., Cameli, P., Otranto, A., Luzzi, S., Bianchi, F., et al. (2022). Cytokine profiles in the detection of severe lung involvement in hospitalized patients with COVID-19: the IL-8/IL-32 axis. Cytokine 151:155804. doi: 10.1016/j.cyto.2022.155804

PubMed Abstract | CrossRef Full Text | Google Scholar

Cao, Y., Yisimayi, A., Jian, F., Song, W., Xiao, T., Wang, L., et al. (2022). BA.2.12.1, BA.4 and BA.5 escape antibodies elicited by omicron infection. Nature 608, 593–602. doi: 10.1038/s41586-022-04980-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Castano-Rodriguez, C., Honrubia, J. M., Gutierrez-Alvarez, J., DeDiego, M. L., Nieto-Torres, J. L., Jimenez-Guardeno, J. M., et al. (2018). Role of severe acute respiratory syndrome coronavirus Viroporins E, 3a, and 8a in replication and pathogenesis. MBio 9, e02325–e02317. doi: 10.1128/mBio.02325-17

CrossRef Full Text | Google Scholar

Conesa, A., Götz, S., García-Gómez, J. M., Terol, J., Talón, M., and Robles, M. (2005). Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics 21, 3674–3676. doi: 10.1093/bioinformatics/bti610

PubMed Abstract | CrossRef Full Text | Google Scholar

Coperchini, F., Chiovato, L., Croce, L., Magri, F., and Rotondi, M. (2020). The cytokine storm in COVID-19: An overview of the involvement of the chemokine/chemokine-receptor system. Cytokine Growth Factor Rev. 53, 25–32. doi: 10.1016/j.cytogfr.2020.05.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Feng, Z., Shen, Y., Li, S., Li, J., Wang, S., Zhang, Z., et al. (2022). The first outbreak of omicron subvariant BA.5.2 - Beijing municipality, China, July 4, 2022. China CDC Wkly. 4, 667–668. doi: 10.46234/ccdcw2022.136

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, H. Y., Lin, Y. C., Li, J., Huang, K. Y., Shrestha, S., Hong, H. C., et al. (2020). miRTarBase 2020: updates to the experimentally validated microRNA-target interaction database. Nucleic Acids Res. 48, D148–d154. doi: 10.1093/nar/gkz896

PubMed Abstract | CrossRef Full Text | Google Scholar

Iketani, S., Liu, L., Guo, Y., Liu, L., Chan, J. F. W., Huang, Y., et al. (2022). Antibody evasion properties of SARS-CoV-2 omicron sublineages. Nature 604, 553–556. doi: 10.1038/s41586-022-04594-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, H., Wu, C., Xu, W., Chen, H., Fang, F., Chen, M., et al. (2022). First imported case of SARS-CoV-2 omicron subvariant BA.5 - Shanghai municipality, China, may 13, 2022. China CDC Wkly. 4, 665–666. doi: 10.46234/ccdcw2022.104

PubMed Abstract | CrossRef Full Text | Google Scholar

Ketter, P., Yu, J. J., Cap, A. P., Forsthuber, T., and Arulanandam, B. (2016). Pentraxin 3: an immune modulator of infection and useful marker for disease severity assessment in sepsis. Expert. Rev. Clin. Immunol. 12, 501–507. doi: 10.1586/1744666x.2016.1166957

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, D., Langmead, B., and Salzberg, S. L. (2015). HISAT: a fast spliced aligner with low memory requirements. Nat. Methods 12, 357–360. doi: 10.1038/nmeth.3317

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, H. S., Lee, H., Park, J., Abbas, N., Kang, S., Hyun, H., et al. (2022). Collection and detection of SARS-CoV-2 in exhaled breath using face mask. PLoS One 17:e0270765. doi: 10.1371/journal.pone.0270765

PubMed Abstract | CrossRef Full Text | Google Scholar

Kong, D.-H., Kim, Y. K., Kim, M. R., Jang, J. H., and Lee, S. (2018). Emerging roles of vascular cell adhesion Molecule-1 (VCAM-1) in immunological disorders and cancer. Int. J. Mol. Sci. 19:1057. doi: 10.3390/ijms19041057

PubMed Abstract | CrossRef Full Text | Google Scholar

Kong, L., Zhang, Y., Ye, Z. Q., Liu, X. Q., Zhao, S. Q., Wei, L., et al. (2007). CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. 35, W345–W349. doi: 10.1093/nar/gkm391

PubMed Abstract | CrossRef Full Text | Google Scholar

Langmead, B., and Salzberg, S. L. (2012). Fast gapped-read alignment with bowtie 2. Nat. Methods 9, 357–359. doi: 10.1038/nmeth.1923

PubMed Abstract | CrossRef Full Text | Google Scholar

Liang, R., Liang, L., Zhao, J., Liu, W., Cui, S., Zhang, X., et al. (2022). SP1/miR-92a-1-5p/SOCS5: a novel regulatory axis in feline panleukopenia virus replication. Vet. Microbiol. 273:109549. doi: 10.1016/j.vetmic.2022.109549

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, J., Wang, F., Du, L., Li, J., Yu, T., Jin, Y., et al. (2019). Comprehensive genomic characterization analysis of lncRNAs in cells with Porcine Delta coronavirus infection. Front. Microbiol. 10:3036. doi: 10.3389/fmicb.2019.03036

PubMed Abstract | CrossRef Full Text | Google Scholar

Martin, M. (2011). CUTADAPT removes adapter sequences from high-throughput sequencing reads. EMBnet. J. 17:10. doi: 10.14806/ej.17.1.200

CrossRef Full Text | Google Scholar

Pertea, M., Pertea, G. M., Antonescu, C. M., Chang, T. C., Mendell, J. T., and Salzberg, S. L. (2015). StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat. Biotechnol. 33, 290–295. doi: 10.1038/nbt.3122

PubMed Abstract | CrossRef Full Text | Google Scholar

Regla-Nava, J. A., Nieto-Torres, J. L., Jimenez-Guardeño, J. M., Fernandez-Delgado, R., Fett, C., Castaño-Rodríguez, C., et al. (2015). Severe acute respiratory syndrome coronaviruses with mutations in the E protein are attenuated and promising vaccine candidates. J. Virol. 89, 3870–3887. doi: 10.1128/jvi.03566-14

PubMed Abstract | CrossRef Full Text | Google Scholar

Robinson, M. D., McCarthy, D. J., and Smyth, G. K. (2010). edgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139–140. doi: 10.1093/bioinformatics/btp616

PubMed Abstract | CrossRef Full Text | Google Scholar

Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 2498–2504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, G., Cui, Q., Garcia, G. Jr., Wang, C., Zhang, M., Arumugaswami, V., et al. (2021). Comparative transcriptomic analysis of SARS-CoV-2 infected cell model systems reveals differential innate immune responses. Sci. Rep. 11:17146. doi: 10.1038/s41598-021-96462-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, L., Luo, H., Bu, D., Zhao, G., Yu, K., Zhang, C., et al. (2013). Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. 41:e166. doi: 10.1093/nar/gkt646

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, Y.-S., Xu, F., An, Q., Chen, C., Yang, Z.-N., Lu, H.-J., et al. (2020). A SARS-CoV-2 variant with the 12-bp deletion at E gene. Emerg. Microb. Infect. 9, 2361–2367. doi: 10.1080/22221751.2020.1837017

PubMed Abstract | CrossRef Full Text | Google Scholar

Trapnell, C., Williams, B. A., Pertea, G., Mortazavi, A., Kwan, G., van Baren, M. J., et al. (2010). Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat. Biotechnol. 28, 511–515. doi: 10.1038/nbt.1621

PubMed Abstract | CrossRef Full Text | Google Scholar

Valadan, R., Golchin, S., Alizadeh-Navaei, R., Haghshenas, M., Zargari, M., Mousavi, T., et al. (2022). Differential gene expression analysis of common target genes for the detection of SARS-CoV-2 using real time-PCR. AMB Express 12:112. doi: 10.1186/s13568-022-01454-2

PubMed Abstract | CrossRef Full Text | Google Scholar

WHO (2022). WHO Coronavirus Disease (COVID-19) Dashboard. Geneva, Switzerland: World Health Organization. Available at: https://covid19.who.int/. (Accessed December 4, 2022].

Google Scholar

Wu, T., Hu, E., Xu, S., Chen, M., Guo, P., Dai, Z., et al. (2021). clusterProfiler 4.0: a universal enrichment tool for interpreting omics data. Innovation 2:100141. doi: 10.1016/j.xinn.2021.100141

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiong, L., Cao, J., Yang, X., Chen, S., Wu, M., Wang, C., et al. (2022). Exploring the mechanism of action of Xuanfei Baidu granule (XFBD) in the treatment of COVID-19 based on molecular docking and molecular dynamics. Front. Cell. Infect. Microbiol. 12:965273. doi: 10.3389/fcimb.2022.965273

PubMed Abstract | CrossRef Full Text | Google Scholar

Yao, P., Zhang, Y., Sun, Y., Gu, Y., Xu, F., Su, B., et al. (2020). Isolation and growth characteristics of SARS-CoV-2 in Vero cell. Virol. Sin. 35, 348–350. doi: 10.1007/s12250-020-00241-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, N., Zhang, D., Wang, W., Li, X., Yang, B., Song, J., et al. (2020). A novel coronavirus from patients with pneumonia in China, 2019. N. Engl. J. Med. 382, 727–733. doi: 10.1056/NEJMoa2001017

PubMed Abstract | CrossRef Full Text | Google Scholar