Comparative Transcriptome Analysis Reveals the Intensive Early Stage Responses of Host Cells to SARS-CoV-2 Infection

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has caused a widespread outbreak of highly pathogenic coronavirus disease 2019 (COVID-19). It is therefore important and timely to characterize interactions between the virus and host cell at the molecular level to understand its disease pathogenesis. To gain insights, we performed high-throughput sequencing that generated time-series data simultaneously for bioinformatics analysis of virus genomes and host transcriptomes implicated in SARS-CoV-2 infection. Our analysis results showed that the rapid growth of the virus was accompanied by an early intensive response of host genes. We also systematically compared the molecular footprints of the host cells in response to SARS-CoV-2, SARS-CoV, and Middle East respiratory syndrome coronavirus (MERS-CoV). Upon infection, SARS-CoV-2 induced hundreds of up-regulated host genes hallmarked by a significant cytokine production, followed by virus-specific host antiviral responses. While the cytokine and antiviral responses triggered by SARS-CoV and MERS-CoV were only observed during the late stage of infection, the host antiviral responses during the SARS-CoV-2 infection were gradually enhanced lagging behind the production of cytokine. The early rapid host responses were potentially attributed to the high efficiency of SARS-CoV-2 entry into host cells, underscored by evidence of a remarkably up-regulated gene expression of TPRMSS2 soon after infection. Taken together, our findings provide novel molecular insights into the mechanisms underlying the infectivity and pathogenicity of SARS-CoV-2.


Calu-3 Cell Infections and RNA Isolation
All experiments involving infectious virus were performed in approved biosafety level 3 (BSL) laboratories at the National Institute for Viral Disease Control and Prevention, China CDC. Cells were washed with MEM and inoculated with viruses at a multiplicity of infection (MOI) of 5 or mock-diluted in MEM for 2 h at 37 • C. Following inoculation, cells were washed three times with MEM, and fresh medium was added to signify 0 h. Triplicate samples of mock-infected and virus-infected Calu-3 cells were harvested at different times between 0 and 24 h postinfection (hpi). Calu-3 cells were cultured for RNA isolation using TRIzol reagent (Invitrogen, United States) following the manufacturer's protocol.

Virus Titration in Cell Culture Supernatants
The 50% tissue culture infectious dose (TCID50) per ml was determined for SARS-CoV-2 in Vero cells. In short, Vero cells (2 × 10 5 cells/well) were seeded into 96-well plates and infected with cell culture supernatants in a dilution ratio of 1:10. After adsorption for 1 h at 37 • C, the cell-free medium was removed, cells were washed with DMEM, and then fresh medium (DMEM containing 2% FBS) was added to cells. After 72 h, cells were observed to evaluate cytopathic effect. The TCID50 values were calculated using the Reed-Muench equation.

Library Construction and Sequencing
A total amount of 50 ng RNA per sample was used as input material for the total RNA library construction and host rRNA removal according to the instructions of the Trio RNA-Seq kit (Nugen,. Total RNA libraries were sequenced on the Illumina Novaseq using the 2 × 150 bp paired-end read setting.

Data Analysis
Raw reads were filtered to obtain clean data by Trimmomatic (v0.35) (with parameters "ILLUMINACLIP:adapter.fa:2:30:10 HEADCROP:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36") (Bolger et al., 2014). The cleaned data were mapped to the human GRCh38 reference genome using STAR aligner (v2.7.2a) (Dobin et al., 2013). The htseq-count command was used to count reads mapped to each gene (Anders et al., 2015). The R package DESeq2 was applied to further identify differentially expressed genes (DEGs) [false discovery rate (FDR) < 0.05, | log2FC| ≥ 1] (Love et al., 2014). The unmapped reads against the entire human genome were further aligned to the reference genome of SARS-CoV-2 (EPI_ISL_402119). Virus genome annotation was based on our previous work . Gene Ontology (GO) enrichment analysis was performed by Fisher's exact test with 19,932 human protein-coding genes as a background in R. The GO terms of enrichment analysis were generated from the "gene2go" file 1 , in which redundant GO terms were further trimmed by the following criteria. First, if one pair of parent and child GO terms had the same genes, only the parent GO term was kept. Second, if one pair of siblings GO terms had the same genes, one sibling GO term was omitted. Third, only GO terms with high quality evidence codes ("TAS, " "IDA, " "IMP, " "IGI, " "IPI, " "IEP, " "ISS") were used. GO term relationships were extracted from the "gobasic.obo" file 2 . For analysis of microarray data of SARS-CoV and MERS-CoV, normalization and identification of DEGs (FDR < 0.05, | log2FC| ≥ 1) were conducted using the R package limma (Ritchie et al., 2015) by following the steps in our previous study (Sun et al., 2019).

Quantitative Real-Time PCR
The same RNA samples of RNA-Seq were used for quantitative real-time PCR (qRT-PCR). qRT-PCR analysis was conducted on the 7500 Fast (Life Technology) by using HiScript II One Step qRT-PCR SYBR Green Kit (Vazyme) according to the manufacturer's instructions. Among housekeeping genes that were downloaded from this paper (Eisenberg and Levanon, 2013), 10 genes showed no significant differential expression at any time point during SARS-CoV-2 infection, in which ATF4 gene was used as reference gene of qRT-PCR. The gene primers were as follows: ATF4 (F: CTCCGGGACAGATTGGATGTT, R:

Transcriptome Profiling of Virus-Host Interactions Following SARS-CoV-2 Infection
We carried out time-course experiments to identify dynamic changes in transcripts in response to SARS-CoV-2 based on the infected and mock-infected groups across four time points (0, 7, 12, and 24 hpi), in which three biologically independent replicates for each treatment group were used for constructing cDNA libraries. The Calu-3 human airway epithelial cell line, a model of human respiratory disease (Love et al., 2014), was used as the host cell of SARS-CoV-2, subjected to the same MOI and host cell used in the previous analyses of SARS-CoV and MERS-CoV infections. The previous datasets of SARS-CoV and MERS-CoV used a sub-population of Calu-3 (Calu-3 2B4 sorted by ACE2 antibody (Josset et al., 2013)) with better infection ability than Calu-3. In order to allow similar cell entry of SARS-CoV-2, Calu-3 cells were incubated with viruses for 2 h before infection timing. After the total RNA isolation and sequencing, we obtained the host transcriptomes, as well as the genomes and transcripts of viruses. The high-throughput sequencing resulted in an average of 49 million paired-end reads per sample, and the sequencing quality was high with a mean mapping rate of unique reads at approximately 72% among mock samples (Supplementary Table S1). The quality control of all samples was assessed by the principal component analysis (PCA) based on normalized counts from DESeq2, which indicated that high quality was achieved given that the majority of samples were well clustered except for only one sample from the infection group at 24 hpi that was removed before further analysis ( Figure 1A).

Rapid Growth of SARS-CoV-2 Accompanied by Dynamic Changes of Host Genes
To evaluate the growth rate of SARS-CoV-2, we calculated the RNA level of the virus represented by unique reads mapping rates at different time points. Our results showed that in general, the virus reads increased sharply from 1.4 to 61.2%, whereas reads mapped to the host genome dropped rapidly from 67.2 to 11.4% (Figure 1B), suggesting a rapid replication of the virus within 24 h. From the results, at the earliest time point (0 hpi), the virus produced high-levels of viral genome RNA as evidenced by relatively even coverage depth across the whole genome ( Figure 1C). Interestingly, we found that there was a significant active transcription of the 3 end of SARS-CoV-2 at 7 hpi, especially for the M, 6, 7a, 7b, 8b, and N genes ( Figure 1C), which could play important roles in the antagonism with host immune response (Lim et al., 2016;Fung and Liu, 2019). After that, the relatively even depth distribution of reads along viral genome was again observed at panels of 12 and 24 hpi. These time-dependent patterns of virus replication and transcription were most likely to play critical roles in the pathology of SARS-CoV-2.
To elucidate the global changes of host gene expression along with virus growth, we identified the overall up-and down-regulated DEGs during SARS-CoV-2 infection ( Figure 1D and Supplementary Table S2). As shown in Figure 1D, during the early stage of infection before 7 hpi, there were many more up-regulated genes than down-regulated genes (498 vs. 212 at 0 hpi, 71 vs. 11 at 7 hpi), and soon after, the number of downregulated genes significantly exceeds that of up-regulated genes (924 vs. 1,501 at 12 hpi, 2,473 vs. 3,611 at 24 hpi). Although the numbers of DEGs were influenced by cutoffs of fold change, there were still many significantly dysregulated genes at the early stage of SARS-CoV-2 infection by using different cutoffs (Supplementary Figure S1). Most importantly, most of the DEGs at 0 hpi were suppressed at 7 hpi, which simultaneously occurred with active transcription of the 3 end of SARS-CoV-2 genome, demonstrating the critical role of the 3 end in antagonizing host immune response. The suppression of host responses was not likely due to sequencing bias because the three samples from the infected group at 7 hpi were clustered with mock samples ( Figure 1A). Interestingly, there seems to be some correlation between the decrease in the levels of the host transcriptome (compared with the total RNA level of SARS-CoV-2) and the relative number of up-regulated genes (compared with down-regulated genes) (Supplementary Figure S2). This may indicate the complex molecular behavior of the host cell in response to the virus infection.

Comparison of Host Transcriptome Responses to SARS-CoV-2, SARS-CoV, and MERS-CoV
To investigate specific host responses during SARS-CoV-2 infection, we performed a comparative transcriptome analysis by integrating two public host transcriptomes of SARS-CoV (GSE33267) (Sims et al., 2013) and MERS-CoV (GSE45042) (Josset et al., 2013) infected in the same cell line with the same MOI. Overall, a huge divergence was presented in time-specific DEG patterns among SARS-CoV-2, SARS-CoV, and MERS-CoV ( Figure 1D). For SARS-CoV-2, 710 DEGs (498 up-regulated and 212 down-regulated) were immediately induced at the very early stage (0 hpi), and many more DEGs were gradually observed at the late stages (12 and 24 hpi). In contrast, SARS-CoV-and MERS-CoV-infected cells exhibited far fewer DEGs (0, 4, and 3 for SARS-CoV and 0, 6, and 54 for MERS-CoV) at the early stages (0, 7, and 12 hpi). However, more DEGs were clearly detected at 24 hpi during SARS-CoV and especially MERS-CoV infection (268 and 4,302, respectively). These distinct DEG patterns indicated that SARS-CoV-2 actually induced earlier host responses than SARS-CoV and MERS-CoV.
To further delineate differential perturbation of pathways among three viruses, we conducted GO-enrichment analysis based on their respective DEGs. Overall, substantially enriched pathways, such as inflammation, apoptosis, antiviral response, transcription, translation, and mitochondrion-related pathways, were detected at various time points during SARS-CoV-2 infection (Figure 2). At 0 hpi, the up-regulated DEGs were mostly enriched in the inflammation-related pathways including the nuclear factor kappa B (NF-kB) signaling and cytokine-mediated signaling pathways, suggesting that SARS-CoV-2 could induce inflammatory responses at the very early stage of infection. At the same time, SARS-CoV-2 also triggered the cellular apoptosis signaling pathway, implying that early onset cell death happened along with inflammation response. Beginning at 7 hpi, our results showed a significant enrichment in antivirus responserelated pathways until 24 hpi (Figure 2). At the late stages (12 and 24 hpi), down-regulated DEGs were exclusively enriched in fundamental host pathways responsible for RNA processing and transcription, protein translation, and mitochondrial activity (Figure 2). Moreover, based on gene expression levels at 0 hpi and other time points, all of the DEGs during SARS-CoV-2 infection were divided into 10 different gene patterns, in which early and late dysregulated DEGs corresponded to distinct biological functions, respectively (Supplementary Figure S3 and Supplementary Table S3). Different from SARS-CoV-2, at the late stage of SARS-CoV infection (24 hpi), the highly enriched genes were identified to be involved in antivirusrelated pathways, whereas no significantly enriched pathways were found for MERS-CoV infection despite many DEGs existing at 24 hpi (Figure 2). Taken together, the above results indicated that the etiology mechanism of SARS-CoV-2 was different from that of SARS-CoV and MERS-CoV as implicated by the overall differential patterns of the host response against infection.
Quantification of the Capacity for Host Antiviral Immunity and Cytokine Production for SARS-CoV-2, SARS-CoV, and MERS-CoV Infections As mentioned above, SRAR-CoV-2 induced specific patterns of host antiviral and inflammation responses compared with SARS-CoV and MERS-CoV. To quantify host antiviral capacity and inflammation responses during infection of the three viruses, two sets of genes were used as their indicators. First, we used a set of 45 early induced genes in interferon-α treated Calu-3 cell (Sun et al., 2019) as antiviral indicators to quantify the level of host antiviral capacity against SARS-CoV-2, SARS-CoV, and MERS-CoV infections. Our analysis showed that the antiviral capacity of the host against SARS-CoV-2 was gradually increased over the time course of infection ( Figure 3A). In contrast, the host antiviral capacities against both SARS-CoV and MERS-CoV were nearly zero at least during the initial stages of infection (between 0 and 12 hpi), followed by a marginal increase at 24 hpi. The antiviral capacity in SARS-CoV-and especially MERS-CoVinfected cells was much lower than that in SARS-CoV-2-infected cells, which might underpin the disparity in mortality between the three viruses. Despite the observation of the potent early induced host antiviral activity during SARS-CoV-2 infection as compared with SARS-CoV and MERS-CoV infection, our results clearly showed that most of the genes (25/45) were significantly induced among infections of the three viruses ( Figure 3B). In addition, a list of the virus-specific antiviralrelated genes was identified, including PARP10 (Yu et al., 2011) and CMPK2 Figure 3B).
Second, we further used a set of 113 human cytokines to quantify host inflammation responses between three viruses. The 113 cytokines from the CytoReg database were often cited by various publications and play a primary role in the immune system (Carrasco Pro et al., 2018). Our results showed that, for SARS-CoV-2, the level of cytokine production was highly induced at 0 hpi, decreased at 7 hpi, and then slowly recovered thereafter ( Figure 3C). Relatively high levels of cytokine expression only occurred at 24 hpi for SARS-CoV and MERS-CoV. Our analysis also provided evidence that SARS-CoV-2 had more cytokines in common with SARS-CoV than with MERS-CoV ( Figure 3D). Unlike the other two viruses, MERS-CoV specifically induced the expression of dozens of cytokines, such as LTA, IL19, CXCL13, and CCL3, at 24 hpi, which were not observed in the case of the other two viruses. Interestingly, among the 28 up-regulated cytokines at the very early stage (0 hpi) during SARS-CoV-2 infection, eight cytokines including IL-6 (IL6), IL-1b (IL1B), IL-8 (CXCL8), G-CSF (CSF3), GM-CSF (CSF2), IP10 (CXCL10), MCP1 (CCL2), and TNF were reported to exhibit substantially elevated serum levels Qin et al., 2020;Xu et al., 2020), which indicated that early induction of cytokines played critical roles in the pathology of SARS-CoV-2. While most of the eight cytokines were moderately up-regulated at the late stage during SARS-CoV and MERS-CoV infections, up-regulation was not observed at the early stage. Collectively, SARS-CoV-2 induced distinct patterns of host antiviral response and cytokine production.

Regulation of Key Genes From Cell Entry to Type-I Interferon Production
Next, to gain possible explanations for the distinct patterns in host antiviral capacity and cytokine production during SARS-CoV-2 infection, dynamic expressions of four types of key genes were evaluated, including virus receptors for cell entry, pathogen recognition receptors (PRRs) for an innate immune startup, and regulator genes for induction of antiviral-related genes and interferon production (Figure 4). For the three cell entry-related genes [ACE2 as the receptor of SARS-CoV and SARS-CoV-2 (Li et al., 2003;Hoffmann et al., 2020), DPP4 as the receptor of MERS-CoV (Raj et al., 2013), and protease TMPRSS2 for S protein priming of SARS-CoV-2 (Carrasco Pro et al., 2018)], we observed the dramatic changes in TMPRSS2 expression with very early induction during SARS-CoV-2 infection and the slightly down-regulated expression of ACE2 in cells infected with SARS-CoV-2 and SARS-CoV, whereas DPP4 was more upregulated in MERS-CoV (Figure 4). For the two PRRs, DDX58 is a canonical RIG-I-like receptor for RNA virus recognition (Kell and Gale, 2015), and TLR3 is a Toll-like receptor playing important roles in initiating a protective innate immune response to highly pathogenic coronavirus infections (Totura et al., 2015). We observed that all three viruses had a notably upregulated expression of DDX58, whereas only MERS-CoV had a suppressed TLR3 at 24 hpi (Figure 4), which is consistent with the fact that decreased expression of TLR3 contributes to the pathology of highly pathogenic coronavirus infections (Totura et al., 2015). Among the four regulator genes, IRF7 is responsible for the expression of most IFN-α subtypes and the type I IFN amplification loop (Lazear et al., 2013), and IRF9, STAT1, and STAT2 form the ISGF3 complex that binds to interferonstimulated response elements and thereby induces the expression of interferon-stimulated genes (Schreiber and Piehler, 2015). As expected, gradual up-regulation of the four primary regulator genes was observed for all three viruses (Figure 4). At last, we found an apparent difference in the expression of IFNB1 between SARS-CoV-2, SARS-CoV, and MERS-CoV, indicating that IFNB1 likely accounted for the observed variations of the host antiviral capacities among three viruses (Figure 4). Taken together, early induction of TMPRSS2 and gradual increased expression level of IFNB1 were likely responsible for the distinct host immune response patterns of SARS-CoV-2 infection.

qRT-PCR Validation of Gene Expression Levels
To validate the accuracy of RNA-Seq gene expression levels, qRT-PCR experiments were performed for nine genes including TMPRSS2, ACE2, DPP4, DDX58, IFNB1, IFNAR2, IL6, IL1B, and TNF. Here, a housekeeping gene ATF4 without significant differential expression at any time point during SARS-CoV-2 infection was used as reference gene of qRT-PCR. As shown in Figure 5, expression levels of the nine genes exhibited high consistencies with those of RNA-Seq (Figures 3, 4), which supported that gene expression quantification from RNA-Seq was reliable.

CONCLUSION
Using time-series profiling of the virus genome and host transcriptome at the same time during SARS-CoV-2 infection coupled with comparative transcriptome analysis, we found that, compared with SARS-CoV and MERS-CoV, SARS-CoV-2 induces strong host cell responses at the very early stage of infection that not only favor its high infectivity to host cells but also restrict its pathogenesis.

DISCUSSION
Here, we sequenced the transcriptomes of SARS-CoV-2-and virus-infected host cells simultaneously during the early stages of infection, providing a robust reference dataset to speculate the antagonistic pattern between pathogen and host cells. To summarize, our findings showed that SARS-CoV-2 induced the significantly high expression of the cellular serine protease TMPRSS2 at 0 hpi to help the entry of viral particles into cells (Hoffmann et al., 2020 ; Figure 4). At the same time, host cell initiated an immediate response for the invasion of SARS-CoV-2 virus ( Figure 1D). Then, the virus successfully suppressed the acute response of host cells for fast proliferation by increasing the transcripts of its 3 genome end, including M, 6, 7a, 7b, 8, and N genes, which were consistent with their reported regulations to host immune response (Lim et al., 2016;Fung and Liu, 2019). As a response from host cell, a FIGURE 5 | qRT-PCR validation of gene expression. The label "mock" indicates mock samples at 0 hpi. The labels from 0 to 24 h represent infected samples at indicated time points. Statistical significance is tested using t-test. "*" denotes significant difference and "ns" for no significance. The error bars represent mean ± SD. *p < 0.05, **p < 0.01, ***p < 0.001 and ****p < 0.0001. number of antiviral pathways and cytokine productions were upregulated to resist the virus infection (Figures 2, 3). In particular, several metabolism-associated pathways were down-regulated at 12 and 24 hpi (Figure 2). After the antagonistic cycle, a dramatic proliferation of viral genome or transcriptome RNA was detected in the early infection of host cells (Figure 1B), which could possibly be an explanation for the fast spread of SARS-CoV-2 in humans. In addition to suppression of host responses by viral proteins, there was a possibility that the replication speed of SARS-CoV-2 could lead to the distinct host cellular response. However, simultaneous detection of SARS-CoV-2 viruses in cell culture supernatants at indicated infection time points showed that the replication level of SARS-CoV-2 (Supplementary Figure S4) seemed similar to those of SARS-CoV and MERS-CoV (Josset et al., 2013).
As SARS-CoV-2 reported a relatively low risk of mortality (Verity et al., 2020) compared with the other two serious human coronaviruses, SARS-CoV and MERS-CoV, we compared and contrasted the host transcriptomes in response to the viral infections. We found that some cytokines in SARS-CoV-2infected cells were markedly up-regulated at a very early stage, which was not observed for SARS-CoV and MERS-CoV and even less frequently observed for other viruses. The unusual high expression of cytokines at 0 hpi possibly explains why patients with severe clinical symptoms rapidly deteriorated. Although the number of infected cases was very high, the majority of infections displayed mild symptoms that are partly explained by a gradual increase in host antiviral capability from 7 to 24 hpi. In contrast to SARS-CoV-2, both SARS-CoV and MERS-CoV were able to inhibit the antiviral capability of the host significantly, which could explain their observed relatively high mortalities. MERS was associated with a higher mortality than SARS, which could be in part attributed to the higher expression of cytokines suppressing the antiviral responses.
Recently, Blanco-Melo et al. (2020) have published transcriptome data of host responses to SARS-CoV-2 from in vitro cell lines including A549 (MOI of 0.2) and NHBE (MOI of 2) at 24 hpi. This previously published data is complemented by our study designed to investigate the early response phase of cell lines infected with SARS-CoV-2. While the previous work did not observe the elevated levels of IFNB1, IFNL1, and IFNL3, our findings show that not only IFNB1 but also IFNL1 and IFNL3 expressions are up-regulated between 7 and 24 hpi ( Figure 3D). Also, they did not detect the gene expression of ACE2 and TMPRSS2 at 24 hpi, whereas we observed that ACE2 is down-regulated at 24 hpi and TMPRSS2 is only up-regulated at 0 hpi before returning to the normal levels (Figure 4). Our time-series sampling revealed distinct early response features of SARS-CoV-2, which provided a possible explanation for some clinical observations. For example, a recent clinical study (Wolfel et al., 2020) found that SARS-CoV-2 could replicate effectively in upper respiratory tract tissues, and that the viral loads appeared earlier (before day 5) and were substantially more than expected. Findings from the present study have confirmed that, at 7 hpi, the 3 end of SARS-CoV-2 genome starts to express densely, reducing the effectiveness of host immune surveillance, which possibly enables the rapid replication of SARS-CoV-2 in upper respiratory tract tissues.
In spite of the fact that several studies have already demonstrated a consistent correlation between gene expression measured by RNA-Seq and by microarray (Nookaew et al., 2012;Guo et al., 2013;Chen et al., 2017), we still need to exclude the possibility of bias resulting from different methodologies. First, because RNA-Seq can potentially detect more genes than microarrays, we only considered protein-coding genes for the analysis of RNA-Seq results. For SARS-CoV-2, the microarray analysis identified more than 90% of the 6,794 DEGs, including 6,514 DEGs of SARS-CoV and 6,198 DEGs of MERS-CoV. Second, expressions of the 6,800 DEGs were distributed over the four time points from low to high, not only in SARS-CoV-2 but also in SARS-CoV and MERS-CoV (Supplementary Figures  S5, S6), indicating that the silent early host responses to SARS-CoV and MERS-CoV appeared not to be due to technological biases. Lastly, when extending the infection time from 24 to 72 hpi (GSE33267), thousands of DEGs (minimum 1,022 and maximum 2,017 genes), which had been inhibited at the early stages, were actually induced (Supplementary Figure S7).
During SARS-CoV-2 infection, reads mapped to host genomes actually decreased, which possibly biased gene expression quantification. However, enough host reads were still obtained for the infected samples (Supplementary  5,130,244 (infected 24 hpi) vs. 18,567,346 (mock 24 hpi). Notably, in the early stages (0 and 7 hpi), about 20 M reads were sufficient for accurate quantification of gene expression level, which can solidly support our main findings in the early stages of infection. Furthermore, qRT-PCR experiments showed consistent gene expression levels between RNA-Seq and qRT-PCR technologies (Figure 5).
With respect to Calu-3 cells, the two public datasets used a sub-population of Calu-3 cells (Calu-3 2B4) that were sorted by ACE2 antibody in order to help virus entry into host cells, whereas our dataset used the mixed Calu-3 cells with low and high ACE2 expressions. To allow similar cell entry, thus, a little longer incubation time (2 h) was taken for SARS-CoV-2 than those for SARS-CoV and MERS-CoV (40 min) (Josset et al., 2013). Although MERS-CoV used a different receptor, DPP4, comparisons among three viruses in Calu-3 cells were still reasonable for two reasons. Firstly, DPP4 had relatively similar expression levels to ACE2 in Calu-3 cells (Supplementary Figure S8). Secondly, SARS-CoV and MERS-CoV had very similar replication kinetics within 24 h (Josset et al., 2013).

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: National Genomics Data Center (https://bigd.big.ac.cn/) with the accession number PRJCA002617. https://bigd.big.ac.cn/bioproject/ browse/PRJCA002617.

AUTHOR CONTRIBUTIONS
TJ, WT, BH, and ZX devised the experiment and wrote the manuscript. JSu, AW, JSh, and WZ conducted bioinformatics analysis. FY, BH, and RY prepared the samples. MW provided the Calu-3 cell. MP did the RNA sequencing and qRT-PCR. All authors revised the manuscript.  -19 Prevention and Control (2020XGZX033, ZX). This manuscript has been released as a pre-print at bioRxiv (Sun et al., 2020).

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2020. 593857/full#supplementary-material Supplementary Figure 5 | Expression levels of cytokine genes between the mock and infected groups. The left column is for SARS-CoV-2, the middle column for SARS-CoV, and the right column for MERS-CoV. For SARS-CoV-2, the expression level is quantified by log2(TPM + 1). Cytokine genes are highlighted.
Supplementary Figure 6 | Expression levels of antivirus-related genes between the mock and infected groups. The left, middle, and right columns show results for SARS-CoV-2, SARS-CoV, and MERS-CoV, respectively. For SARS-CoV-2, the expression level is quantified by log2(TPM + 1). Antiviral-related genes are highlighted.