Ribosome-Profiling Reveals Restricted Post Transcriptional Expression of Antiviral Cytokines and Transcription Factors during SARS-CoV-2 Infection.

The global COVID-19 pandemic caused by SARS-CoV-2 has resulted in over 2.2 million deaths. Disease outcomes range from asymptomatic to severe with, so far, minimal genotypic change to the virus so understanding the host response is paramount. Transcriptomics has become incredibly important in understanding host-pathogen interactions; however, post-transcriptional regulation plays an important role in infection and immunity through translation and mRNA stability, allowing tight control over potent host responses by both the host and the invading virus. Here, we apply ribosome profiling to assess post-transcriptional regulation of host genes during SARS-CoV-2 infection of a human lung epithelial cell line (Calu-3). We have identified numerous transcription factors (JUN, ZBTB20, ATF3, HIVEP2 and EGR1) as well as select antiviral cytokine genes, namely IFNB1, IFNL1,2 and 3, IL-6 and CCL5, that are restricted at the post-transcriptional level by SARS-CoV-2 infection and discuss the impact this would have on the host response to infection. This early phase restriction of antiviral transcripts in the lungs may allow high viral load and consequent immune dysregulation typically seen in SARS-CoV-2 infection.

1. Introduction

Coronaviruses are enveloped positive sense RNA viruses with aene">n exn class="Chemical">ceptionally large genome encoding the structural proteins envelope (E), spike (S), membrane (M), and nucleocapsid (N), in addition to a plethora of non-structural and accessory proteins. Infections with endemic human coronaviruses (e.g., 229E, NL63, OC43, and HKU1) cause a mild common cold, however three novel coronaviruses have emerged from animal reservoirs in the 21st century, SARS-CoV, MERS and SARS-CoV-2 causing a fatal respiratory syndrome in 34%, 15%, and 3% of cases, respectively with SARS-CoV-2 being the most infectious [1]. The current COVID-19 pandemic caused by SARS-CoV-2 emerged in China in December 2019 [2] and has since spread across the globe, causing more than 100 million confirmed cases and over 2.7 million deaths (https://covid19.who.int/ accessed on 24 March 2021).

SARS-CoV-2 is adept at evading innate immunity [3], the naïve host’s primary defense against a newly emerged n class="Species">coronavirus. It does this using numerous structural and non-structural proteins that inhibit interferon (IFN) production and function. Through blocking IFN, viral replication can proceed unchecked. In addition to antiviral functions, circulating IFN alerts the host to viral infection with symptoms such as fever, pain, and fatigue [4]. This capacity for stealth replication in the absence of symptoms allows asymptomatic transmission [5], making it particularly difficult to control from a public health perspective. Understanding the molecular mechanics of immune evasion strategies employed by SARS-CoV-2 will allow us to develop more effective diagnostics, therapeutics, and host prognostic indicators.

RNA sequencing (RNA-seq) enables disn class="Chemical">covery of molecular pathways involved in SARS-CoV-2 infection. RNA-seq of nasopharyngeal swabs from infected individuals revealed up-regulation of antiviral factors OAS1-3 and IFIT1-3 as well as chemokines like IP-10, the latter of which showed muted expression in older individuals [6]. Added to this, transcriptional profiles in SARS-CoV-2 infected cell, animal and patient serum samples revealed low levels of type I and III interferons but high levels of TNF, IL-6, RANTES and CCL20 compared to influenza A infection, indicating differences in the immune responses to SARS-CoV-2 compared to another respiratory virus [7]. These studies have provided a greater insight into potential mechanisms underpinning COVID-19 immunopathogenesis. However, they do not capture post-transcriptional regulatory mechanisms, many of which control the initiation, magnitude, duration, and resolution of the innate immune response [8] or allow the virus to restrict host gene expression [9].

Here, we have used ribosome profiling to assess post-transcriptional regulation of host gene expression during n class="Disease">SARS-CoV-2 infection. We found restricted translation of type I and III IFNs 24 h post infection as well as numerous transcription factors and antiviral cytokines. These data provide a set of genes likely responsible for the delayed antiviral responses seen in SARS-type coronavirus infections strengthening development of molecular-based strategies to combat this devastating pandemic.

2. Results

2.1. Development of a SARS-CoV-2 Infection Model in Calu-3 Cells

Of the multiple organs that SARS-CoV-2 infen class="Chemical">cts, viral load is highest in the lungs, and infection of this organ is a significant driver of pathogenesis [10]. We therefore chose a lung epithelial cell line known to be productively infected with SARS-CoV-2 to investigate host responses to infection. Calu-3 cells are a human immortalized cell line derived from lung adenocarcinoma which express sufficient levels of ACE2 to allow SARS-Cov-2 entry. We found SARS-CoV-2 infectious titres peaked 24 h post infection in Calu-3 cells (Figure 1A) as per other studies [11]. Both SARS-CoV-2 intracellular genomic RNA (Figure 1B) and N protein (Figure 1C) followed similar kinetics with peak viral load post 24 h infection. Immunofluorescence staining of SARS-CoV-2 N-protein over a time-course from 1 to 48 h showed infection of more than 50% of cells after 24 h (Figure 1D).

Figure 1

SARS-CoV-2 infection of Calu-3 cells (A) TCID50 measurements of virus titres, (B) qRT-PCR measurements of intracellular viral RNA represented by 2-ΔΔCt normalised first to GAPDH and then to inoculum levels of SARS-CoV-2, set to 1, and (C) intracellular viral protein in Calu-3 cells infected with SARS-CoV-2 (MOI 0.1). (D) Immunofluorescence microscopy showing SARS-CoV-2 N protein staining (green) in Calu-3 cells infected with SARS-CoV-2 (MOI 0.1) at various timepoints. Cell nuclei were stained using DAPI (blue).

Having developed an infection model in n class="CellLine">Calu-3 cells, we proceeded to capture transcriptional and post-transcriptional changes in gene expression at 24 h post SARS-CoV-2 infection using a method deemed suitable for high containment infections [12]. Cells were snap frozen in liquid nitrogen to lock ribosomes onto their associated mRNA, avoiding artifacts introduced by translation inhibitors like cycloheximide [13]. Both total RNA and RNA protected by protein from micrococcal nuclease (MNase) digest were purified. The latter, named herein as ribosome footprints were further purified to contain fragments between 25 and 35 nt which spans the possible 80S ribosomal footprint on actively translated mRNAs (Figure S1). More than 30 million reads were obtained from next generation sequencing of the total mRNA libraries, while more than 50 million reads were obtained from the MNase-digested small RNA libraries (Figure S2).

2.2. Transcriptional Response to SARS-CoV-2 Infection of Calu-3 Cells Is Dominated by Antiviral Defence Genes

We first analyzed reads obtained from total RNA representing the transcriptome. 11,504 genes passed the threshold of >1 n class="Chemical">counts per million in at least three samples and were analyzed for differential expression by comparing triplicate mock and SARS-CoV-2-infected samples. There were 2.5 times more up-regulated than down-regulated genes (166 versus 63, Figure 2A) using a log2 fold change cut-off of 1 and p-value cut-off of 0.05. Upon inspection of the genes up-regulated by SARS-CoV-2 infection, many genes were previously annotated antiviral genes. To confirm this, the 229 significantly altered genes were submitted to https://david.ncifcrf.gov/, (accessed on 23 March 2021) for functional annotation clustering using the Uniprot keyword database. The top seven significantly enriched keywords are plotted in Figure 2B. Antiviral defense showed the highest fold enrichment and a p-value < 0.001. Other significantly enriched keywords included innate immunity, RNA-binding, and transcription. Of the 134 Antiviral defense keywords, 27 were up-regulated by SARS-CoV-2 and are listed in Figure 2C. Many of these genes, including IFIT1, IFIT2, IFIT3, RSAD2, OASL, HERC5, DDX58 and MX2 were also upregulated in nasopharyngeal swabs of COVID19 patients [6].

Figure 2

Transcriptional response to SARS-CoV-2 infection of Calu-3 cells is dominated by antiviral defense genes. (A) Volcano plot showing global transcriptional changes of ~11,000 genes in SARS-CoV-2 infected Calu-3 cells. Log2FoldChange infected 24 h versus mock. 229 transcripts were differentially expressed (DE) in SARS-CoV-2 infected cells based on the cut-off of p-value < 0.05, log2fold change > 1. (B) Significantly altered genes were submitted to https://david.ncifcrf.gov/ (accessed on 23 March 2021) for Functional Annotation Clustering. Here, Uniprot Keywords were summarised using fold enrichment of DE genes, coloured by PValue. (C) Common genes between Uniprot “Antiviral Defense” keywords and the 229 SARS-CoV-2 DE genes.

2.3. MNase Digestion of Cell Extracts Yields CDS-Mapped Ribosome Footprints

Next, we proceeded to aene">nalyze the reads generated from sequenn class="Chemical">cing of size-selected MNase-digested RNA to understand the translatome in SARS-CoV-2 infected cells. 50% of the reads (approximately 25 million reads) mapped to protein coding sequence in both mock and infected cells (Figure 3A), indicating close to 50% coverage of the 18 million nucleotides of coding sequence predicted in the human transcriptome [14]. This finding implies that sequencing depth was sufficient to make gene level inferences on differential ribosome density. More than 50% of reads mapped to the first 30 nucleotides of tRNAs mostly carrying the small neutral amino acid Glycine and some carrying the small to medium sized acidic amino acids, Glutamate and Aspartate (Figure S3). These tRNA reads and the less abundant reads mapping to other non-coding RNA were filtered out and remaining reads mapped to different features of human protein coding genes. Most of these filtered reads (60%) mapped to within the CDS (Figure 3B). The 3′UTR and 5′UTR had 20% and 7% of the reads, respectively no significant difference between mock and infected cells was observed. Annotated start and stop codons were present in 2% and 6% of the reads, respectively with no difference between mock and infected. When high-resolution metagene analysis of transcript coverage was performed, a gradual increase in coverage was seen leading up the start codon, whereby coverage increased 8-fold to a peak around 30nt consistent with ribosomes initiating at the start codon (Figure S4).

Figure 3

MNAse digestion of cell extracts yields 30 nt long CDS-mapped ribosome footprints. (A) Percentage of reads mapping to indicated RNA species in mock and SARS-CoV-2 infected cells at 24 h. (B) Percentage of non-coding RNA filtered reads mapping to indicated gene-level features in mock and SARS-CoV-2 infected cells. (C,D) Read-length counts per million read distribution in coding sequence and 5′ untranslated regions (UTR).

Eukaryotic 80S ribosome footprints are ~28–30 nun class="Chemical">cleotides in length [15]. To confirm that CDS-mapped reads were bona fide ribosome footprints, read-lengths were calculated for reads mapping to the CDS and the 5′UTR region as 80S ribosomes scan this region to find the start codon [16]. Read-length distributions were summarized across all genes. Read lengths peaked at 28–30 nucleotides in coding sequence in both mock and infected cells indicative of bona-fide ribosome footprints (Figure 3C). Read-lengths higher than 30 nucleotides showed an exponential shortening towards 40 nucleotides suggesting single nucleotide degradation, while those shorter than 30 nucleotides were staggered suggestive of footprints derived from smaller non-ribosomal RNA binding proteins. In the 5′ UTR read lengths also peaked at 30 nucleotides indicative of translation in the 5′ leader as documented elsewhere [17] but at >6-fold lower density compared to the CDS mapped reads (Figure 3D). There was no significant difference in 5′ UTR ribosome occupancy between mock and infected cells apart from greater variance between replicates in the infected samples.

2.4. Transcription Factors and Cytokines Are Blocked at the Post-Transcriptional Level

Using CDS-mapped read n class="Chemical">counts from MNase-digested RNA libraries, we compared the translatome of mock versus SARS-CoV-2 infected Calu-3 cells. Here, 11,455 genes passed the threshold (>1 counts per million in at least three samples), comparing well to numbers passing threshold expression for the transcriptome (Figure 4A). Three genes (KRT17, CDKN1A and TOM1) showed more than a 2-fold decrease in expression. 18 genes showed more than a 2-fold increase in expression upon infection and are mostly antiviral, i.e., Interferon-induced proteins withtetratricopeptiderepeats (IFIT), oligoadenylate synthase (OAS), DDX58 (RIGI) and Interferon-induced GTP-binding protein Mx1 (MX1). Of note is that many genes seen upregulated at the whole transcript level were not seen to be differentially expressed when ribosome-associated RNA was sequenced. This transcriptome-only change contrasts with other stimulatory contexts where 9% of genes showed transcriptome-only upregulation and 85% of genes showed translation-only upregulation [18].

Figure 4

Absence of transcription factors and cytokines from genes up-regulated by ribosome foot-printing. (A) Volcano plot showing global translational changes of ~11,000 genes in SARS-CoV-2 infected Calu-3 cells. Log2FoldChange infected 24 h versus mock. 21 genes were differentially expressed (DE) in SARS-CoV-2 infected cells based on the cut-off of p-value < 0.05, log2fold change > 1. (B) Functional annotation clustering on the differential expression gene lists for RNA-seq versus Ribo-seq. Uniprot keywords used to classify genes showed an absence of keywords ‘Signal’, ‘Secreted’ and ‘Transcription regulation’ in the Ribo-seq gene list. (C) Heat-map of the log2 fold change of up-regulated RNA-seq transcripts included within ‘Innate immunity’, ‘Transcription’ and ‘Cytokine keywords’ alongside their log2 fold change as determined by Ribo-seq.

We investigated this paucity of up-regulated genes further, using funn class="Chemical">ctional annotation clustering on the differential expression gene lists for RNA-seq versus Ribo-seq. Using Uniprot keywords to classify genes, showed an absence of keywords ‘Signal’, ‘Secreted’ and ‘Transcription regulation’ in the Ribo-seq gene list (Figure 4B). These three annotation groups include genes such as IFNL1, CCL5, IFNB1, IL6, ATF3, JUN, SETD1B and REL all commonly known cytokines and transcription factors. We then focused our attention on three Uniprot keywords; ‘innate immunity’ because many up-regulated genes at the transcript level were in this group, ‘transcription’ because of the discordance between Ribo-seq and RNA-seq, and ‘cytokine’ for the same reason. We selected genes in these groups that were significantly up-regulated by RNA-seq and present their log2 fold change as determined by Ribo-seq (Figure 4C). Within the ‘Innate immunity’ annotation group, most genes upregulated by RNA-seq were also upregulated by Ribo-seq, except for NLRC5, a negative regulator of NFKB and IFN-I signaling [19], which was 2.2-fold and 6.2-fold induced by viral infection for RNA-seq and Ribo-seq, respectively. Within the ‘Transcription’ annotation group, 37 (84%) of the differentially expressed genes by RNA-seq were restricted in expression by Ribo-seq, with one of the most restricted being JUN (3.7-fold increase by RNA-seq and 1.3-fold decreased by Ribo-seq). Those un-restricted at the Ribo-seq level included ZNF107, TXNIP, IRF1 and TRIM22. In the third annotation group ‘Cytokine’, all three Interferon lambdas (IFNL1, IFNL2, IFNL3), Interferon beta, IL6, CCL5 and CXCL11 were restricted at the Ribo-seq level. Cytokines with no restriction at the Ribo-seq level were TNFSF10 (TRAIL), CXCL10 (IP10), CX3CL1 and CSF1.

2.5. Coverage of Individual Genes Confirms Post-Transcriptional Restriction of IL6, CCL5, JUN and IFNB1

To examine ribosome ocn class="Chemical">cupancy pattern across the length of interesting genes we collected sample normalized, per nucleotide coverage statistics for individual transcripts. Coverage values for all transcript variants were summed at each position from the transcription start site. Gene-level coverage for RNA-seq reads was even across the 5′UTR, CDS and 3′UTR with a decrease observed at the 3′-end termini due to terminal exon-skipped transcript variants (Figure 5). Conversely, Ribo-seq coverage was minimal >100 nucleotides up and down from the start and stop codon, respectively, and this was most evident in genes that were well-occupied by ribosomes, i.e., TNFSF10. IRF7 showed noticeable Ribo-seq coverage in the full length of the 5′ and 3′UTRs which were uniquely short. Genes presented here that were upregulated when expression was assayed by RNA-seq and not by Ribo-seq included IL6, JUN, IFNB1 and CCL5 (RANTES). These genes show a noticeable increase in RNA-seq coverage over the CDS from mock to 24h SARS-CoV-2 infection (top right quadrant, Figure 5) consistent with log2 fold change calculated from count data. The transcription factor REL, showed equivalent up-regulation by RNA-seq and Ribo-seq coverage contrary to the count data, which showed a 2.5-fold increase by RNA-seq (p-value = 2.267 × 10−7) versus 1.3-fold increase by Ribo-seq (not significant). Unlike the aforementioned genes, IRF7, CXCL10, TNFSF10 (TRAIL), OAS2 and IRF1 showed up-regulation by RNA-seq and Ribo-seq. This was not due to higher levels of peak coverage as CXCL10 only reached 30 stacked reads per nucleotide (coverage), equivalent to CCL5 which was not concordant between RNA-seq and Ribo-seq in relative changes between mock and 24 h SARS-CoV-2 infected. Neither were Ribo-seq reads in the 3′UTR associated with discordance as both IL6 and IRF1 showed peaks > 250 nucleotides downstream of the reference sequence stop codon. In summary, while transcript levels for genes like IL6, CCL5, IFNB1 and JUN were upregulated, these transcripts were not increased in their association with ribosomes suggesting a restriction in the antiviral response.

Figure 5

Coverage of select genes confirms post-transcriptional restriction of IL6, JUN, IFNB1, CCL5, REL, IRF7, CXCL10, TNFSF10, OAS2 and IRF1. Sample normalized, per nucleotide coverage statistics for indicated transcripts. Coverage values for all transcript variants were summed at each position from the transcription start site. CDS cordinates (blue) were determined using the default transcript per gene as outlined in the Matched Annotation from NCBI and EMBL-EBI (MANE) project. Red indicates 5′ (left) and 3′ (right) untranslated regions.

2.6. mRNA Features Influence the Likelihood of Translation Inhibition

2.6.1. Unstable mRNAs Are More Sensitive to Translation Inhibition in SARS-CoV-2 Infected Calu-3 Cells

Differential expression analysis of Ribo-seq n class="Chemical">count data allowed us to examine the genes which up-regulated in ribosomal occupancy and is a closer approximation to protein output and phenotype than RNA-seq [20]. However, to understand the mechanism through which a restriction in ribosome occupancy occurs, we need to separate the two processes of transcription and translation. Translation efficiency considers the background mRNA levels. We utilized a method, Riborex, whereby RNA-seq count data obtained from matched biological samples is factored into a generalized linear model for the Ribo-seq counts to find genes with differential translational efficiency [21].

Many transcription fan class="Chemical">ctors and secreted proteins, which encompass cytokines were found to have short half-lives in mouse embryonic stem cells [22]. To assist in explaining why some cytokines and transcription factors were sensitive to translation inhibition while others were not, we categorized genes based on their mRNA stability as documented by Sharova et al. [22]. 6878 genes with known mRNA half-lives and which passed the default threshold values for Riborex differential expression analysis are presented in Figure 6A. By separating genes into stable mRNAs (half-life > 5 h) and unstable mRNAs (half-life < 5 h) we found 0.5% of stable mRNAs and 0.2% of unstable mRNAs were increased in translation efficiency. Conversely, more unstable mRNA genes were decreased in translation efficiency, 3.5% compared to 0.74% of stable mRNAs. The genes with the greatest log2 fold decrease, JUN, ZBTB20 and IL6 were all unstable mRNAs. On average, unstable (mRNA half-life < 5 h) genes showed a significant decrease in translation efficiency upon infection with SARS-CoV-2 compared to stable mRNAs (one-way analysis of variance [ANOVA], p < 0.001; Figure 6B).

Figure 6

mRNA features influence sensitivity to translation inhibition. (A) Volcano plot showing global translation efficiency (Riborex engine) changes of 6878 genes in SARS-CoV-2 infected Calu-3 cells with previously documented mRNA stability [22]. Log2FoldChange infected 24 h versus mock. 135 transcripts were differentially expressed (DE) in SARS-CoV-2 infected cells based on the cut-off of p-value < 0.05, log2fold change > 1. Genes were categorized as having half-life more or less than 5 h as indicated. (B) Log2 fold change as presented in (A) of unstable (<5 h mRNA half-life) and stable (>5 h mRNA half-life) genes presented using geom_boxplot function in the ggplot2 R package using the 25th and 75th percentiles to form the box and whiskers no larger than 1.5 times the interquartile range. Data points beyond the whiskers are outliers. *** indicates one-way analysis of variance [ANOVA], p < 0.001. (C) Volcano plot showing global translation efficiency changes of 9751 genes in SARS-CoV-2 infected Calu-3 cells, 2659 of which have one or more high-confidence upstream open reading frames (uORFS) obtained from [26]. (D) Log2 fold change for genes with no uORF or one or more uORF as presented in (C) using geom_boxplot as described in (B). *** indicates one-way analysis of variance [ANOVA], p < 0.001. (E) Venn diagram showing genes significantly decreased in translation efficiency by 24 h of SARS-CoV-2 infection in Calu-3 cells (blue), genes increased in translation during EIF2A repression [25] (yellow) and genes significantly increased in translation efficiency by 24 h of SARS-CoV-2 infection in Calu-3 cells (green).

2.6.2. mRNAs with uORFs Are Not Protected from Translation Inhibition

Infection of n class="CellLine">293/ACE2 cells with SARS-CoV-1 lead to activation of PKR followed by phosphorylation of and hence inactivation of the translation initiation factor EIF2A [23]. This event forms the core of the Integrated Stress Response (ISR) causing global inhibition of translation [24]. A hallmark of this response is the preferential translation of stress-response genes (~200 genes in mouse embryonic fibroblasts) owing to the presence of upstream open reading frames (uORFs) in their 5′UTRs [17,25]. To determine if uORFs conferred protection from translation inhibition seen here, we sorted genes based on the presence of an upstream open reading frame (uORF) in the 5′ UTR as catalogued for 2659 human genes using Ribo-seq datasets from THP-1 and HEK293 cells [26]. JUN, ZBTB20 and IL6 all have no uORFs, while 1.5% of downregulated genes contained no uORF and 1.8% containing at least one uORF, respectively were decreased in translation efficiency (Figure 6C). On average, genes containing an uORF showed a small but significant decrease in translation efficiency upon infection with SARS-CoV-2 compared to genes without uORFs (one-way analysis of variance [ANOVA], p < 0.001; Figure 6D). Since mRNAs predicted to contain at least one uORF are more likely to be decreased instead of increased in translation efficiency, uORFs are not protective against translation inhibition by SARS-CoV-2 infection. To further explore this, we obtained a list of 68 homologous genes known to increase in translation during EIF2A phosphorylation (EIF2A-P) in mouse embryonic fibroblasts [25] and contain a documented uORF in humans [26]. We then compared this list with significantly altered genes from our mock versus SARS-CoV-2 infected comparison of translation efficiency. We found no EIF2A-P enhanced genes in the list of genes upregulated in translation efficiency by SARS-CoV-2 infection (Figure 6E). Eighteen genes were downregulated in translation efficiency by SARS-CoV-2 but known to increase in the presence of EIF2A-P in mouse embryonic fibroblasts. This finding gives further weight to the hypothesis that downregulation of translation efficiency by SARS-CoV-2 is independent of the EIF2A-P stress response.

3. Discussion

Largely driven by signaling molecules n class="Chemical">called interferons (IFNs), the innate immune system first recognizes non-self, pathogen-associated molecular patterns (PAMPs) such as double stranded RNA through host proteins DDX58 (RIG-I) and Toll-like receptor 3/7 (TLR3/7). This recognition initiates phosphorylation and activation of transcription factors IRF3, AP1 complex and NFKB. In infected lung epithelial cells, these transcription factors activate transcription of genes which can be categorized into three groups. (i) Type I and III IFNs (IFNB and IFNL) which increase transcription of interferon stimulated genes (ISGs) including IRF7 which fuels further IFN transcription [27], (ii) Proinflammatory cytokines TNF, IL-6 and IL-1β which prepare bystander immune cells to fight the infection and (iii) chemokines like CCL5 (RANTES) and CXCL10 (IP-10) to attract distal immune cells to the site of infection. IFNs are arguably the most immediate of these three antiviral effectors, binding to their cognate receptors on already infected and neighboring cells, activating kinases JAK and TYK2, which then phosphorylate STAT1, STAT2, MAPK, PI3K allowing their translocation to the nucleus for activation of target genes with diverse antiviral functions [27].

SARS-CoV-2 is remarkably sensitive to type I n class="Gene">IFN pre-treatment [11], suggesting that limiting IFN production, rather than function is the main mechanism through which SARS-CoV-2 achieves high viral loads concomitant with delayed or absent symptoms. Non-structural protein 1 (Nsp1) was found to be one of the most potent inhibitors of IFNB1 promoter induction of 26 viral proteins tested [28]. Nsp1 acts to inhibit host translation by blocking the mRNA entry tunnel on the 40S ribosome [29]. This interaction can also lead to endonucleolytic cleavage in the 5′ UTR preventing the recruitment of more ribosomes to host mRNAs [30,31].

Studies in bronchial epithelial n class="Chemical">cells show posttranscriptional restriction of IFNL1/2, IFNB1, CCL5 and IL6 by SARS-CoV-1 and not Dhori virus, an orthomyxovirus known to also productively infect 2B4 cells [32]. Here, we corroborate these findings in Calu-3 cells using an unbiased genome-wide methodology showing that antiviral cytokines (IFNL1-3, IFNB1, CCL5, CXCL11 and IL6) are restricted in translation at 24 h post infection. This is consistent with the delayed IFN response seen in COVID19 [33]. Importantly, in a mouse model of SARS-CoV, IFN was only detected in the lung after viral load had already peaked rendering this arm of innate immunity entirely ineffective [34].

A host driven delay in IL6 protein produn class="Chemical">ction is unlikely because poly I:C treatment of Calu-3 cells results in a 2-12-fold increase in secreted IL-6 in just 4 h, with the upper range dependent upon a Th1/17 cytokine environment [35]. IL6 is a significant driver of fever [36], which was an initial symptom in 41% of 174 health care workers that tested positive for SARS-CoV-2 infection upon routine screening [37]. Delayed IL6 production could explain the relatively slow onset of symptoms in SARS-CoV-2 infection. Medium time to symptom onset is 5.1 days and up to 14 days [38], compared to influenza A 1.4 days, influenza B 0.6 days, rhinovirus 1.9 days, parainfluenza virus 2.6 days, non-SARS human coronavirus 3.2 days and respiratory syncytial virus 4.4 days, [39]. Moreover, early IL6 and CCL5 signaling enhances recruitment of innate immune cells and resolution lung pathology caused by respiratory syncytial virus [40,41] suggesting that timing is important for the antiviral activity of these pro-inflammatory cytokines.

In addition to the abovementioned antiviral cytokines, we show that selen class="Chemical">ct transcription factors are restricted at the posttranscriptional level, with JUN, ZBTB20, ATF3, HIVEP2 and EGR1 showing the greatest restriction when accounting for the increase in substrate mRNA while REL and CREB5 were significantly upregulated by RNA-seq but not Ribo-seq when analysis was performed on un-normalized Ribo-seq reads. JUN regulates transcription of IFNB1 mRNA as part of the enhanceosome, encompassing ATF2, JUN, IRF3, p50, p65, CBP and p300 [42]. Post-transcriptional restriction of JUN could limit enhanceosome formation and IFN induction. Indeed, weak induction of IFNB1 mRNA was seen by SARS-CoV-2 in Calu-3 cells compared to Sendai virus [43]. ZBTB20 promotes Toll-like receptor innate immune responses by inhibiting transcription of IKBA [44], decreased protein expression of ZBTB20 may therefore delay immune responses.

In contrast, n class="Gene">ATF3 is anti-inflammatory, reducing IL6 and IL12B transcription as part of a negative feedback loop [45]. Restriction in ATF3 protein expression may contribute to un-controlled IL6 production in cells that overcome the translational block, either through competition with ever increasing amounts of IL6 mRNA [46] or through an increase in IL6 mRNA stability enacted by the balance of RNA binding proteins ZC3H12A (Regenase-1) and ARID5A [47]. The latter of which stabilizes IL6 and is negatively regulated by MAPK14 (p38 MAPK) to resolve inflammation [48]. The increase in MAPK14 (p38 MAPK) signaling in aged tissues may contribute to uncontrolled IL6 production [49]. HIVEP2 may also be anti-inflammatory, acting as an NFKB inhibitor, enabling cell survival during memory T cell development [50]. Whether HIVEP2 translational restriction contributes to immune dysregulation or enhances SARS-CoV-2 in infected airway epithelial cells is unclear.

EGR1 does have a n class="Chemical">clear antiviral function, suppressing foot-and-mouth disease virus (FMDV), Sendai virus and Seneca Valley virus through phosphorylation of TBK1, an important kinase that activates IRF3 [51]. Pos-transcriptional restriction of EGR1 may contribute to delayed IFN responses. Transcription factors REL and CREB5 were significantly upregulated by RNA-seq but not Ribo-seq. REL in heterodimers with RELA is a key transcriptional activator of IFNL1 in human airway epithelial cells [52], thus translational blockage could potentially limit IFNL1 transcription. However, Finally, CREB5 downregulation allows viral persistence in nasopharyngeal epithelia of FMDV-infected animals via a decrease in chemokine transcription suggesting it may also be antiviral [51].

RNA binding proteins, PARP14 aene">nd n class="Gene">ZC3HAV1 (PARP13) also showed post-transcriptional restriction. Loss of PARP14 strongly attenuates inducible transcription of IFBN1 [53] and Coronaviruses encode a macrodomain protein (Nsp3) that inhibits PARP14 to counteract the host antiviral IFN response [54]. Translational repression of PARP14 may further limit this anti-coronaviral protein to support viral replication. ZC3HAV1 plays important roles in the posttranscriptional regulation of mRNA during the stress response enacting degradation of host and viral mRNAs, translational repression, and microRNA silencing [55]. ZC3HAV1 is known to be antiviral to a suite of viruses including retrovirus, alphavirus, filovirus and hepadnavirus by degrading viral RNA. It is also pro-apoptotic by inducing degradation of TNFRSF10D (TRAIL Receptor 4) mRNA, a decoy receptor in the TRAIL-induced apoptosis pathway. Thus, post-transcriptional restriction of ZCSHAV1 would likely favor SARS-CoV-2 replication.

Many cytokines, traene">nsn class="Chemical">cription factors and growth factors are time limited in their effects by RNA binding proteins recruited to AU-rich elements (ARE) in their 3′ untranslated region (UTR). Conditions where ARE-containing mRNAs are stabilized for too long can lead to chronic inflammation [56]. For example, decreased expression of the ARE binding protein AUF-1 in chronic obstructive pulmonary disease (COPD) patient samples leads to stabilization of IL-6, CCL2, CCL1 and CXCL8 mRNA [57]. Many of the genes we found to be inhibited in translation were unstable. JUN, ZBTB20, ATF3 and HIVEP2 have half-lives of 1.3, 2.91, 2.17 and 2.03 h in contrast to IRF7 which has a half-life of 15 h, while IL6 has a half-life of only 1 h. Stability is similar in human lung epithelial cells where IL6 half-life was 40 min in the presence of TNF and 80+ minutes in the presence of IL-17 [58]. Similarly, JUN mRNA half-life ranged from 45 to 90 min in a human monocytic cell line [59], while ZBTB20, ATF3 and HIVEP2 mRNA half-life was ≤5 h in HEK293 cells [60].

Assuming a closed loop model of traene">nslation fan class="Chemical">cilitated by the eIF4G (5′ cap) to poly(A)-binding protein (3′ tail) interaction [61], stabilized mRNAs are more likely to have ribosomes that re-initiate on the same mRNA, a process coined closed-loop assisted reinitation (CLAR) [62]. Mammalian ribosomes translate at approximately 6 s per amino acid [63] and must be spaced ~100 nucleotides apart [64]. Thus, IL6 would only have time to reinitiate approximately four ribosomes per mRNA, HIVEP2 only one, while IRF7 would survive long enough to house its maximum of 15. Consistent with this, mRNAs bound by multiple ribosomes have higher mRNA stability than those bound by a single ribosome [65]. mRNAs using CLAR are less dependent upon scanning mechanisms performed by EIF4F/EIF4A and initiate more efficiently [66]. Mechanisms of translation inhibition occurring prior to mRNA engagement, i.e., virally encoded Nsp1 [46] and host encoded 4E-BP1 [67] should be less effective on closed loop mRNAs. While there currently lacks direct evidence that CLAR ribosomes resist translation inhibition by viral Nsp1 or host 4E-BP1, the increased representation of unstable mRNAs in downregulated genes shown here, combined with a potential for reduced access to the mRNA entry tunnel in CLAR ribosomes suggests this might be the case and warrants further investigation.

We found that mRNAs containing at least one high-n class="Chemical">confidence uORF were more likely to be decreased in translation efficiency. Given that genes containing an uORF are preferentially translated in the presence of phosphorylated EIF2A, we propose that the integrated stress response is not directing posttranscriptional restriction in SARS-CoV-2 infected Calu-3 cells at 24 h post infection. Moreover, eighteen genes known to be increased in translation in the presence of phosphorylated EIF2A were significantly decreased in translation efficiency by SARS-CoV-2. Additionally, the best characterized effector of the ISR, ATF4 [68], was not preferentially translated during SARS-CoV-2 infection (Table S3). In fact, none of the ISR genes identified in mouse embryonic stem cells were preferentially translated. This absence of an ISR signature suggests either EIF2A is not being phosphorylated in this context or the effects on translation are being over-ridden by Nsp1. 4E-BP1 remains a possible driver of translation inhibition in this context however Rapamycin, an activator of 4EBP1, leads to inhibition of ribosomal proteins and elongation factor translation [69] which encompass 17% of the upregulated genes (Table S3; RPL32, RPL16, RPL10, RPL23A, RPL6, RPL41, RPL13A). Thus Nsp1 is likely responsible for the restriction and further experiments will aim to confirm this.

Research into the n class="Species">human immune response to SARS-CoV-2 infection is still in its infancy. We have limited patient data on the early molecular response because when patients present to the clinic, they have been infected for days already, yet this response is so important for viral control and transmission reduction. Here, we have used RNA-seq and Ribo-seq in human bronchial epithelial cells to assess the early response to SARS-CoV-2 infection at the level of transcription and translation. We found a robust antiviral host response at the level of transcription with upregulation of IFNB1 and IFNLs, Interferon-induced proteins with tetratricopeptide repeats, oligoadenylate synthase, antiviral cytokines and transcription factors following SARS-CoV-2 infection of Calu-3 cells, complementing the analysis done by Finkel at al on translation of SARS-CoV-2 viral open reading frames [70] and strengthening existing transcriptomic data [71]. The biological effects of genes are mediated by protein not mRNA. We have increased the biological relevance of transcriptome analysis by demonstrating restricted translation of IFNB1, IFNL1-3, CCL5, CXCL11, IL6, JUN, ZBTB20, ATF3, HIVEP2, REL, PARP13, PARP14 and EGR1, which have direct antiviral and/or immune regulatory functions. Moreover, we found that unstable mRNAs are more sensitive to SARS-CoV-2 induced translation inhibition and propose that the inhibition may occur at the level of mRNA engagement with the ribosome. Upstream open reading frames were not protective suggesting a mechanism independent of EIF2A repression, which incidentally occurs post mRNA engagement.

In summary, we identify a selection of aene">ntiviral aene">nd immunologin class="Chemical">cal genes that are restricted early in SARS-CoV-2 infection of human cells and highlight that mRNA stability contributes to gene selectivity. Since cytokine mRNA stability is altered in chronic inflammation, assessment of translation restriction during SARS-CoV-2 infection in these contexts is warranted. This may help us understand the wide range of prognoses during infection. Presently, this study provides detailed molecular insight into the delayed interferon response employed by SARS-CoV-2.

4. Materials and Methods

4.1. Cell Culture

VeroE6 n class="Chemical">cells (ATCC CRL-1586) were maintained in Gibco Dulbecco’s Modified Eagles Medium (DMEM) supplemented with 10% (v/v) foetal calf serum (FCS), 100 U/mL penicillin, and 100 μg/mL streptomycin (Life Technologies, Australia). Calu-3 cells were maintained in Gibco Modified Eagles Medium (MEM) supplemented with 20% (v/v) FCS, 10 mM HEPES, 0.1 mM non-essential amino acids, 2 mM glutamine, 1 mM sodium pyruvate, 100 U/mL penicillin, and 100 μg/mL streptomycin (Life Technologies, Australia). All cells were kept at 37 °C in a humidified incubator (5% CO2).

4.2. Infections

All virology work was condun class="Chemical">cted at the CSIRO Australian Centre for Disease Preparedness at physical containment (PC)-4. The isolate of SARS-CoV-2 (BetaCoV/Australia/VIC01/2020) [72] was received from the Victorian Infectious Disease Reference Laboratory (VIDRL, Melbourne, Australia) and passaged in VeroE6 cells for isolation, followed by passaging in VeroE6 cells for stock generation. All virus stocks were aliquoted and stored at −80 °C for inoculations. The infectious titres of SARS-CoV-2 stocks was determined by TCID50 assays performed as described previously [66]. Samples were titrated in quadruplicate in 96-well plates, co-cultured with VeroE6 cells for four days and monitored for development of cytopathic effects (CPE). Calu-3 cells were fixed for 30 min in 4% paraformaldehyde (PFA) and stained with a polyclonal antibody targeting the SARS-CoV-2 Nucleocapsid (N) protein (Sino Biological, China catalogue number: 40588-T62, used at 1/2000) for 1 h. Cells were subsequently stained with 1/1000 dilution of an anti-rabbit AF488 antibody (Invitrogen catalogue number A11008). Nuclei were counter-stained with diamidino-2-phenylindole (DAPI). Calu-3 cells were imaged using the CellInsight quantitative fluorescence microscope (Thermo Fisher Scientific, Australia) at a magnification of 10×, 49 fields/well, capturing the entire well. The relative viral antigen staining was quantified using the Compartmental analysis bioapplication of the Cellomics Scan software. Calu-3 cells were grown to 80% confluency in T25 culture flasks prior to infection with SARS-CoV2 (MOI 1) or mock infected for RNA sequencing experiments.

4.3. MNase Digest and RNA Extraction

Isolation of whole-cell RNA aene">nd purifin class="Chemical">cation of ribosome-protected fragments (RPFs) was performed based on the protocol described by Reid et al. [12]. Cells were harvested on a dry ice-ethanol slurry in 300 μL ice-cold lysis buffer (1% (v/v) IGEPAL, 200 mM KOAc, 25 mM K-HEPES pH 7.2, 10 mM MgCl, 4 mM CaCl2) and incubated on ice for 10 min. The lysate was then clarified by centrifugation (8600× g, 5 min, 4 °C). 50 μL clarified lysate (supernatant) was set aside for RNA-seq analysis while 200 μL was digested with 300 U/mL micrococcal nuclease S7 (MNase; Sigma Aldrich, Sydeny, Australia) for 30 min at 37 °C to generate ribosome-protected fragments (RPFs) from mRNA. These conditions do not allow rRNA cleavage into 25–35 nt fragments due to its highly structured nature compared to mRNA transcripts. Both digested and undigested RNA was then purified by phenol:chloroform extraction using TRI reagent® (Sigma Aldrich, Australia). The air-dried pellet containing undigested RNA was resuspended in nuclease-free water prior to RNA sequencing.

4.4. Phosphatase Treatment and Size Selection

The air-dried pellet containing the extran class="Chemical">cted MNase-digested RNA was treated with T4 polynucleotide kinase (Promega, Madison, WI, USA) to generate compatible ends for sequencing. RPFs were purified from undigested RNA via electrophoretic separation on a TBE-urea polyacrylamide gel (ThermoFisher, Madison, WI, USA) and excision of the region between 25 and 35 nt using an RNA marker. 400 μL 400 mM NaOAc, pH 5.2, was added to the excised gel and incubated for 10 min at −80 °C. RPFs were then purified from the gel through three successive cycles of incubation at 95 °C for 5 min followed by vortexing for 20 min, before extraction using a 0.45 μm cellulose acetate column (Sigma, Australia) by centrifugation (10 min, 20,000× g). 1 mL ethanol and 30 μg/mL GlycoBlue™ co-precipitant was then added to the flowthrough and incubated for 1 hour at −20°C to precipitate the RPFs, followed by centrifugation (15 min, 20,000× g, 4 °C). The resultant pellet containing RPFs was washed twice in ice-cold 80% (v/v) ethanol, with centrifugation (5 min, 20,000× g, 4 °C) following each wash. The RPF pellet was then air-dried and resuspended in nuclease-free water prior to RNA sequencing. The quality and quantity of RNA was assessed for all samples using a Bioanalyzer (Agilent, Santa Carla, CA, USA).

4.5. RNA Sequencing

RNA-Seq was performed by the Australian Genome Research Fan class="Chemical">cility (AGRF). Illumina TruSeq Stranded mRNA and small RNA libraries were prepared according to manufacturer’s instructions, followed by sequencing on an Illumina Novaseq-6000. Raw data were assessed for overall quality using fastqc v0.11.8. (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/, accessed on 23 March 2021).

4.6. Bioinformatic Analysis of RNA-Seq Reads

Quality and adapter trimming was performed using TrimGalore v0.6.4 (http://www.bioinformatics.babraham.an class="Chemical">c.uk/projects/trim_galore/, accessed on 23 March 2021) with default settings for automatic adapter detection. Adaptor trimmed and quality filtered reads were mapped to the human genome (GENCODE v35 primary assembly of GRCh38.p13) using STAR aligner (version 2.5.3a) [73]. Reads mapping to coding sequence (CDS) were counted to quantify transcripts capable of being translated into protein using featureCounts. The Bioconductor package DESeq2 package in R (version 3.6.3) was used to test for differential expression between different experimental groups [74].

4.7. Bioinformatic Analysis of Ribo-Seq Reads

Adaptor was trimmed from the reads using Cutadapt followed by filtering for quality using FastX-toolkit as appropriate for small read aene">nalysis. Trimmed aene">nd filtered reads were then mapped using Bowtie version 1.2.2 (http://bowtie-bio.sourn class="Chemical">ceforge.net/manual.shtml, accessed on 23 March 2021) to the SARS-CoV-2 genome (accession number MT007544.1), human protein coding transcripts downloaded from BioMart (https://m.ensembl.org/info/data/biomart/index.html, accessed on 23 March 2021), miRNAs from miRbase (http://www.mirbase.org, accessed on 23 March 2021), rRNA from the silva database (https://www.arb-silva.de, accessed on 23 March 2021), snoRNA from the snoRNA-LBME-db database and tRNA from the GtRNAdb database [75]. Percentage of reads with at least one reported alignment were collected from alignment outputs and plotted using ggplot2 in R. Reads mapping to non-coding RNA were filtered out and mapped to the human genome (GENCODE v35 primary assembly of GRCh38.p13) using STAR aligner as per RNA-seq reads with the following parameters–outFilterMismatchNmax 2–quantMode TranscriptomeSAM GeneCounts–outSAMattributes MD NH–outFilterMultimapNmax 1. Reads mapping to CDS were counted and differential expression analysis run as per RNA-seq. For gene feature mapping statistics, Samtools version 1.10.0 [76] was used to generate mapping stats on reads aligned to the indicated features (again obtained from BioMart) after filtering out reads mapping to non-coding RNA. The exception was stop and start codon sequence which was obtained using the STAR alignment files followed by featureCounts as per CDS counts but using stop_codon for the -t flag. For coverage analysis, transcriptome-mapped bam files from the STAR alignment were sorted and indexed using Samtools. Sorted bam files were converted to bed files and read depth at each genome position with 1-based coordinates determined using Bedtools version 2.29.2 [77] using a library normalization factor obtained from DESeq2 analysis of the different experimental groups. For metagene analysis of coverage relative to the start codon, the distance from the transcription start site for each gene was extracted from the GENCODE v35 primary assembly of GRCh38.p13 gtf file then used to summarize coverage for across genes using scripts in R. Read-length distributions were obtained using Samtools stats on position sorted alignment files from reads mapped to human CDS and 3′UTR sequences downloaded from Biomart (https://m.ensembl.org, accessed on 23 March 2021), read-length distributions were summarized across all genes. Differential expression analysis of translation efficiency was performed using Riborex package in R using DESeq2 as the modelling engine [21].