Cytokine systems approach demonstrates differences in innate and pro-inflammatory host responses between genetically distinct MERS-CoV isolates

Topological characterization of the host response shows spatio-temporal transcriptomic differences between MERS-CoV strains

To further assess the magnitude impact of MERS-CoV infection in Calu-3 2B4 cells, we integrated viral genomic RNA levels into a topological assessment of host response differences. MERS-CoV genomic RNA was measured by qRT-PCR using primer (KL200F) located in the leader region (nt 9-28) and primer (KL201) located in the ORF1 region of the CoV genome (Figure 1A and Additional file 1: Table S5). Additionally, we measured viral titers by plaque assay and found MERS-CoV SA 1 had more robust infectious viral particle production than MERS-CoV Eng 1, with significantly different viral titers at 18 and 24 hpi, despite similar viral genomic RNA levels at these late time points (Additional file 2: Table S4). In measuring viral genomic RNA and infectious viral particle production, we assessed two different stages in the viral life cycle and the data support differences in the kinetics of replication for the two MERS-CoV strains. Since our analysis focuses on the host response to infection and its origin in pathogen recognition and activation of innate immune responses, we computationally integrated viral genomic RNA with host mRNA to better assess the differences in host response kinetics. Towards this end, we used a filtered clustering method [17] that reduced the space of 93 samples to a clustering graph of 18 nodes. First, the viral genomic RNA levels were used as a filter to bin the samples into overlapping subsets with similar viral load (i.e., viral gRNA). Second, within each subset, samples that showed a high degree of interconnection were retained. Thus, there was high similarity in gene expression levels within each subset. For a graphical representation, each set of retained samples formed a node and nodes were considered as adjacent whenever they had at least one sample in common (Methods and Figure 1B). Mock-infected samples were all in one isolated cluster. With increasing viral load, the temporal delay between MERS-CoV Eng 1 and MERS-CoV SA 1 started as early as 3 hours post infection (hpi), with 3 and 7 hpi MERS-CoV Eng 1 samples clustering with 7 and 12 hpi MERS-CoV SA 1 samples. Genomic RNA levels at 18 and 24 hpi between the two viruses were not statistically different as determined by the non-parametric Mann-Whitney test. Despite similar viral gRNA at 18 and 24 hpi (Figure 1A), the filtered clustering approach demonstrated MERS-CoV Eng 1 and MERS-CoV SA 1 were distinct at the late time points, with 18 and 24 hpi MERS-CoV Eng 1 samples separated from 24 hpi MERS-CoV SA 1 samples (Figure 1B).

Taken together, these whole transcriptome topological analyses showed robust differences in both kinetics and magnitude of the host response between MERS-CoV SA 1 and MERS-CoV Eng 1, despite similar viral replication. Another example where viral replication may not be the best readout of host response differences is icSARS-CoV. The deletion of the entire ORF6 gene of icSARS-CoV (icSARS-CoV ΔORF6) does not impact viral replication compared to wild-type icSARS-CoV, yet substantial differences in the host response between the two viruses are detected in Calu-3 2B4 cells [18].

MERS-CoV SA 1 and MERS-CoV Eng 1 elicit striking differences in epithelial adherens junction signaling and pro-inflammatory cytokine signaling late in infection

Differential gene expression between MERS-CoV Eng 1 and MERS-CoV SA 1 was directly examined by comparing fold changes with respect to time- and experiment-matched mocks (absolute log2 FC > 1, FDR-corrected p-value < 0.05). We found the most robust transcriptional changes to MERS-CoV SA 1 and MERS-CoV Eng 1 infections occurred at the late times points (18 and 24 hpi). Compared to MERS-CoV SA 1 there was an accelerated host response to MERS-CoV Eng 1, with differential gene expression observed as early as 7 hpi. As shown in Figure 2, MDS representation of MERS-CoV SA 1, MERS-CoV Eng 1 and icSARS-CoV samples for 6432 DE genes changing in at least one virus condition and at one time-point confirmed the early onset of host responses to MERS-CoV Eng 1 infection and the sustained separation between MERS-CoV Eng 1 and MERS-CoV SA 1 beginning at 18 hpi. The MDS representation also showed MERS-CoV Eng 1 samples in an intermediate position between MERS-CoV SA 1 and icSARS-CoV, with the host response to icSARS-CoV being largely delayed until 30 and 36 hpi. There is increasing evidence that interferon (IFN) antagonism and avoidance of interferon-stimulated gene (ISG) effector functions are major contributors to the delayed host response in icSARS-CoV infection [18, 19]. In this context, the intermediate positioning of MERS-CoV Eng 1 between MERS-CoV SA 1 and icSARS-CoV led us to investigate MERS-CoV strain-specific differences in innate immune signaling early in infection and the subsequent impact on the later host response.

The greatest number of gene expression differences between MERS-CoV SA 1 and MERS-CoV Eng 1 were induced at the late time points, with 2160 genes and 2611 genes differentially expressed (DE) at 18 and 24 hpi, respectively (Figure 3A). There was a total of 4861 DE genes between MERS-CoV SA 1 and MERS-CoV Eng 1 in at least one time point (herein known as contrasting genes). In comparison to MERS-CoV SA 1, there was an earlier host response to MERS-CoV Eng 1 at 7 and 12 hpi. The accelerated Calu-3 response to MERS-CoV Eng 1 may be the result of the difference in kinetics of viral gRNA replication, with MERS-CoV Eng 1 more efficiently replicating at 3 and 7 hpi compared to MERS-CoV SA 1 (Figure 1A). Alternatively, the delay in the host response to MERS-CoV SA 1 may be due to the virus more efficiently evading innate immune responses that leads to enhanced viral replication compared to MERS-CoV Eng 1 (Additional file 2: Table S4). Between the two viruses there are amino acid differences in ORF4a and PLpro, viral proteins known to modulate the innate immune response during MERS infection [5–7]. We further examined host gene expression at the early time-points (3, 7 and 12 hpi) and found significant enrichment of genes associated with the STAT3 pathway. STAT3 pathway genes CDC25A, MYC, SOCS3 and SOCS4 were more strongly induced by MERS-CoV Eng 1 compared to MERS-CoV SA 1, particularly at 7 hpi (Figure 3B). Increased SOCS gene expression and decreased expression of PIM1 gene in response to MERS-CoV Eng1 indicated decreased STAT3 activity and possibly differential induction of apoptosis-related pathways. In a direct virus comparison, differential expression of pro-apoptotic BID, BAX, and BIM genes was observed at the early time-points and by 24 hpi there was extensive cytopathic effects caused by both viral infections that was likely the result of caspase-dependent apoptosis, as previously shown for MERS infection [20]. Of the 4861 contrasting genes, 2653 genes were also DE against time- and data set-matched mocks. Hierarchical clustering of the 4861 DE genes resulted in distinct gene clusters, with striking expression pattern contrasts between MERS-CoV SA 1 and MERS-CoV Eng 1 at 18 and 24 hpi. As shown in Table 2, functional analysis of the five most prominent clusters with contrasting gene expression revealed enrichment of integrin linked kinase (ILK) signaling and epithelial adherens junction signaling pathways, glutathione metabolism, and interferon and pro-inflammatory signaling pathways. Genes related to glutathione metabolism included GSTM1 and GSTM3, which were strongly downregulated in response to both viruses. ISGs, IFIT1 and IFIT3, were highly induced in response to both MERS-CoV Eng 1 and MERS-CoV SA 1, with pronounced early up-regulation specifically in response to MERS-CoV Eng 1. Genes associated with ILK signaling and epithelial adherens junction signaling pathways were strongly downregulated in response to MERS-CoV SA 1, whereas MERS-CoV Eng 1 predominantly up-regulated expression of these genes at the late time points (Additional file 3: Figure S1). Within the highly enriched pathways, cellular genes including PVRL1, RHOF and CREBBP with highest expression contrasts between the two infections were chosen for confirmation by qRT-PCR (Additional file 4: Figure S2). Pro-inflammatory CSF3 gene was found more highly induced by MERS-CoV SA 1 (log ratios in MERS-CoV SA1 and MERS-CoV Eng1 respectively: 2.2 and 0.5), whereas CCL5 gene was more highly induced by MERS-CoV Eng 1 (log ratios in MERS-CoV SA1 and MERS-CoV Eng1 respectively: 0.71 and 2.2). Examination of the aforementioned contrasting genes showed that only a small number of those were already differentially expressed at the early time-points. For example, expression of RHOF gene in response to MERS-CoV Eng 1 was increased 4-fold relative to mock at 12 hpi, whereas MERS-CoV SA 1 did not induce RHOF gene expression at this time-point. We therefore focused on differences in the host response mainly at the later time-points that had the highest number of contrasting genes.

MERS-CoV Eng 1 and MERS-CoV SA 1 elicit differences in the timing of cytokine-mediated innate immune and pro-inflammatory responses in Calu-3 cells

While there was no difference in IFN- β gene expression detected between the two MERS-CoV strains, there was differential IFN-α2 and IFN- γ gene expression. Specifically, expression of IFN-γ was downregulated early after infection in MERS-CoV Eng 1-infected Calu-3 cells and up-regulated in MERS-CoV SA 1-infected Calu-3 cells. To some extent, this was also observed for IFN-α2 gene expression, which showed higher up-regulation in response to MERS-CoV SA 1 at 24 hpi, as confirmed by qRT-PCR (Additional file 4: Figure S2). In addition, IFN-λ1 and IFN-λ2 were DE relative to mock with high up-regulation at 18 and 24 hpi for both viruses.

There were no virus-specific differences in expression of pro-inflammatory cytokines, IL-1α or TNFα-IP3, which were highly up-regulated in response to both MERS-CoV Eng 1 and MERS-CoV SA 1. The cell migration-promoting factor, TNFα-IP2 [21], was highly up-regulated for MERS-CoV Eng 1 alone (Additional file 4: Figure S2). Cytokines appear to be important for MERS infection. To explore potential mechanisms regulating cytokine activity in response to MERS-CoV, we designed a microarray experiment to analyze Calu-3 responses to various cytokine treatments and develop signatures that were then examined in the context of MERS-CoV-infected Calu-3 cells. Calu-3 cells were treated with either human recombinant interferon IFN-α, IFN-γ, TNF or IL-1α, and cell lysates collected at different time points post-treatment for microarray. We found 399 DE genes responsive to recombinant IFN-α and 261 DE genes responsive to recombinant IFN-γ, and as expected, the majority of these genes were up-regulated following stimulation [22]. In response to pro-inflammatory cytokine treatment, we found 76 DE genes responsive to recombinant TNF and 383 DE genes responsive to recombinant IL-1α, which were also differentially expressed in response to MERS-CoV SA 1 or MERS-CoV Eng 1. Many of these genes showed cytokine-specific expression for IFN-α (260 DE genes), IFN-γ (107 DE genes), and IL-1α (209 DE genes (Figure 4A). TNF induced the least number of DE genes compared to the other cytokines (8 DE genes). Within the cytokine-stimulated genes we identified a set of 149 genes showing strong contrasts at late time points (18 and 24 hpi) between MERS-CoV SA 1 and MERS-CoV Eng 1 (Figure 4B).

Topological analysis of gene expression restricted to IFN and pro-inflammatory cytokine genes sets revealed IFN induction (score: 0.74) as a major contributor to kinetic differences between MERS-CoV SA 1 and MERS-CoV Eng 1 infection (Table 1). The relative numbers of cytokine genes, which were also DE after MERS-CoV SA 1 or MERS-CoV Eng 1 infection (22% for IFN-α, 31% for IFN-γ, 31% for IL-1α and 45% for TNF), indicated the contribution of these specific gene sets to global gene expression changes in MERS infection. Functional analysis of the cytokine DE genes showed antigen presentation pathway was significantly enriched in MERS-CoV-induced IFN-γ responses, as well as enrichment of genes associated with adaptive immune responses such as OX40 signaling and cytotoxic T lymphocyte-mediated apoptosis pathways. IFN-α response genes were significantly enriched for cellular pathways involved in the antiviral response, such as interferon regulatory factor (IRF) activation by pathogen recognition receptors (PRRs). Pro-inflammatory mediators, IL-1α and TNF, which stimulate genes related to IL-17 signaling, were identified following MERS infection. For example, pro-inflammatory cytokine genes, CXCL1, CCL2, and CCL20, regulated by IL-17A [23] and play a prominent role in airway inflammation and disease, and NF-κB genes, NFKBIA, NFKB2 and NFKBIE, were all found to be DE in response to TNF and MERS-CoV infection.

Predicted role for STAT3 in mediating gene expression changes between MERS-CoV Eng 1 and MERS-CoV SA 1 in Calu-3 infected cells

We also identified STAT3 as a predicted upstream regulator at 24 hpi in response to MERS-CoV Eng 1 (z-score = 2.559), with up-regulation of downstream target genes including BCL6, SOCS3, BCL3, IRF7 and PML (Figure 5, Table 3). Through a separate binding motif prediction analysis using JASPAR and PSCAN, we confirmed the STAT3 prediction based on enrichment of STAT3 binding sequences in DE genes at the late time points (Table 3). Conversely, MERS-CoV SA 1 infection resulted in downregulation of these same target genes, suggesting a predicted inhibition or delay in STAT3 regulator activity, though STAT3 did not reach a statistically significant z-score (z-score = 0.364). Among the STAT3 target genes BCL3 had the highest contrasting expression between the viruses. The BCL3 gene was highly up-regulated in response to MERS-CoV Eng 1 (FC = 1.5) and downregulated in response to MERS-CoV SA 1 (FC = -1.66) at 18 hpi. qRT-PCR analysis of BCL3 mRNA expression, and several additional STAT3 target genes, including BCL6, PML, and IRF7, were in agreement with the microarray findings (Additional file 4: Figure S2). Among these genes, IRF7 plays an important role in the cellular antiviral response and Faure and colleagues reported high levels of IRF7 gene expression in BAL cells from a patient that recovered from MERS infection [12].

STAT3 acts as a critical regulator of cellular repair processes upon acute lung injury and BCL3 has been shown to reduce lung inflammation in mice by regulating granulocytes [27]. In addition to BCL3, several other contrasting STAT3 target genes, including BCL6, IL11, and PML, have been reported to impact pro-inflammatory responses [28–30], epithelial integrity and the severity of lung injury after infection [27, 31]. We postulate the difference in pro-inflammatory cytokine gene expression may be the result of differential STAT3 activity. For example, SARS-CoV directly impacts STAT3 activity by dephosphorylating STAT3 at Tyr-705 after 18 hpi in infected Vero E6 cells [32]. Here, the MERS-CoV strain-specific differences may be a causative factor in leading to differential STAT3 activation and the down-stream effects of pro-inflammatory responses.

Cells and viruses

Calu-3 2B4 cells were cultured in minimal essential media (MEM; Gibco) supplemented with 20% fetal bovine serum (HyClone) and 1% antibiotic antimycotic (Gibco). Human coronavirus MERS-CoV SA 1 (GenBank: JX869059.2) was received from Bart L. Haagmans (Erasmus Medical Center, Rotterdam, Netherlands). MERS-CoV SA 1 was isolated from the sputum of a 60-year-old man in Saudi Arabia who died after developing acute respiratory distress syndrome (ARDS) in June 2012 [1]. Human coronavirus MERS-CoV Eng 1 (GenBank: KC667074.1) was received from the United Kingdom Health Protection Agency (HPA) Imported Fever Service. MERS-CoV SA 1 and MERS-CoV Eng 1 were propagated in VeroE6 cells and supernatants clarified. For questions of written informed consent for participation in the study, we refer to [1] and [2], where isolation of viruses from clinical samples is described. Virus stocks were tittered on VeroE6 cells using standard methods, as previously described [14]. SARS-CoV experiments were derived from the infectious clone of SARS-CoV (icSARS-CoV) as described in [18, 38].

Calu-3 2B4 infections

All work was performed in a biosafety level 3 (BSL3) facility at the University of North Carolina-Chapel Hill. Cells were washed with phosphate-buffered saline (PBS) and inoculated with virus at a multiplicity of infection (MOI) of 5 plaque-forming unit (PFU) per cell or mock diluted in PBS for 40 min at 37°C. Following inoculation, cells were washed 3 times, and fresh medium was added. Triplicate infected Calu-3 2B4 cultures and triplicate time-matched mock-infected controls were harvested at 0, 3, 7, 12, 18 and 24 hpi.

Cytokine treatment of Calu-3 cells

Calu-3 cells were treated with either IFN-α (10 ng/ml) (Sigma I4276), IFN-γ (500 U/ml) (Sigma I32265), IL-1α (0.001 ng/ml) (Peprotech, #200-01A) or TNF (0.05 ng/ml) (Peprotech, #300-01A) in triplicate for each cytokine. Treated cell lysates (n = 3) and time-matched mock-treated cell lysates (n = 3) were harvested at 0, 3, 6, and 18 hr following IFN-α treatment, at 0, 3, 6, and 21 hr following IFN-γ treatment, at 3, 6, and 24 hr following IL-1α treatment, and at 6 hr following TNF treatment. Total RNA was extracted from each sample and gene expression analyzed by microarray. Probe labeling and microarray slide hybridization for pooled replicates for each treatment per time-point was performed using Human Genome CGH Microarray (G4447A; Agilent Technologies).

RNA isolation and microarray processing

RNA isolation from Calu-3 2B4 cells infected with MERS-CoV SA 1 and subsequent hybridization to Agilent 4 × 44K human HG arrays was previously reported in [11]. RNA isolation from Calu-3 2B4 cells infected with MERS-CoV Eng 1 was performed following the same protocol used for MERS-CoV SA 1-infected Calu-3 samples and used for subsequent hybridization to Agilent 8 × 60K human arrays. Calu-3 2B4 microarray experiments with an infectious clone recombinant SARS-CoV (icSARS-CoV) (AY278741) was reported in [18]. Using Agilent QC criteria we removed one 12 hpi replicate infected with MERS-CoV SA 1 and one 7 hpi replicate infected with MERS-CoV Eng 1. Due to the absence of mock samples at 24 hpi in the MERS-CoV Eng 1, we used the mock samples harvested at 18 hpi for the comparison at 24 hpi. In order to merge MERS-CoV Eng 1 and MERS-CoV SA1 data sets, we considered only probes that were mapped to gene names (Refseq IDs) in both arrays. To remove batch effects from the merged data sets we applied the R function ComBat [39]. The resulting data set was then quantile normalized. For each sample, a log2 fold change value was calculated as a difference between log2 normalized data for each sample and the average of log2 normalized data for time- and data set-matched mock-infected samples. The same procedure, but handled as separate data sets, was applied to the cytokine treatment microarray data sets.

Statistical analysis

Statistical analyses for viral genomic RNA and viral titer changes in Calu-3 2B4 cells infected with either MERS-CoV SA1 or MERS-CoV were performed by using two-sided non-parametric Mann-Whitney U test (Additional file 2: Table S4).

Topology analysis

For normalized intensity data, the Euclidean distance matrix between samples was clustered with respect to sample-matched viral genomic RNA measurements using an algorithm described in [17]. We used a freely available MATLAB implementation of this algorithm (http://appliedtopology.org/) with the parameter setting filterSamples = 5, pverlapPct = 65, mFudge = 5. The topology induced by a distance matrix can be described by a filter function, which assigns to every sample a value (phenotypic outcome, viral load, etc.). First, the topological space is decomposed into overlapping subsets using level sets of the filter function. The set of samples falling into a particular subset are clustered, and the cluster tree is partitioned into two parts using statistical criteria of the single linkage length. Only the part of the cluster tree with small linkage lengths (dense part) is retained. If resulting filtered clusters of different subsets have at least one sample in common, then they are defined as adjacent in a global discrete cluster graph. Topological differences between MERS-CoV SA 1 and MERS-CoV Eng 1 gene expression during the course of infection was furthermore assessed using homology persistence barcodes (only 0- and 1-homology were taken into account) as described in [15]. For calculations, we used the freely available R package phom [40] and their maximum weighted bipartite graph matching [41] as a measure of difference (see R script implementation in Additional file 5). Additionally, multi-dimensional scaling (MDS) was performed on the Euclidean distance matrix between samples after differential gene expression statistical analysis [42].

Differential gene expression analysis

Differential expression was determined by comparing MERS-CoV SA 1- and MERS-CoV Eng 1-infected replicates to time- and data set-matched mock-infected controls, based on a linear model fit for each probe using the R package Limma [43]. The same method was applied to determine differential expression between strains using time-matched MERS-CoV SA 1- and MERS-CoV Eng 1-infected samples. Criteria for differential expression were an absolute log2 fold change of 1 and a q value of <0.05 calculated using a moderated t test with Benjamini-Hochberg correction. The cytokine Calu-3 treatment data sets were analyzed separately using the same methods and cut-offs as the virus-infected Calu-3 cells. Differentially expressed genes were determined to be statistically significant in at least on time point and for at least one treatment.

Quantitative reverse transcription PCR (qRT-PCR)

RNAs were reverse transcribed using the QuantiTect Reverse Transcription Kit (Qiagen). The resulting cDNA samples were diluted 50X and run on the 7900HT Fast Real-Time PCR System (Applied Biosystems) using Power SYBR Green PCR Master mix (Life Technologies) and 200 nM primers. Primer sets were designed using Primer3 [44]. The primer sequences for cellular gene targets and the strand-specific primers for quantifying viral genomic RNA are listed in Additional file 1: Table S5. For cellular gene mRNA quantification, relative gene expression in infected samples compared to that in mock-infected samples was calculated using the 2^-ΔΔCT method [45]. The RPL14 gene was selected as an internal control (calibrator) due to non-significant changes in RPL14 gene expression in MERS-CoV SA 1- and MERS-CoV Eng 1-infected Calu-3 2B4 cells, as determined by microarray analysis.

Functional enrichment analysis

Functional analysis of statistically significant gene expression changes was performed using Ingenuity Pathways Knowledge Base (IPA; Ingenuity Systems). For all gene set enrichment analyses, a right-tailed Fisher’s exact test was used to calculate P-values associated with each biological function and canonical pathway. Upstream Regulator Analysis in IPA was used to predict regulators and infer their activation state based on prior knowledge of expected effects between regulators and their known target genes according to the Ingenuity Knowledge Base (IKB). The calculated z-score signifies whether gene expression changes for known targets of each regulator are consistent with what is expected from the literature (z > 2, regulator predicted to be activated, z < -2, regulator predicted to inhibited). In addition to IPA Upstream Regulator Analysis, we also used EnrichR [46], a collaborative gene list enrichment tool, to determine possible upstream regulators via overrepresented transcription factor binding site motifs [47].

Causal Network Analysis in IPA was used to understand gene expression changes and causal relationships between genes and networks of upstream regulators in the experimental dataset. The genes within the causal network represent nodes and the edge that defines the biological relationship between two nodes is represented as an arrow signifying regulation. Dashed arrows represent indirect relationships and solid arrows represent direct relationships. All edges are supported by at least one published reference or from canonical information stored in IKB.

Availability of supporting data

Raw microarray data have been deposited in NCBI’s Gene Expression Omnibus and are accessible through GEO accession GSE33267 (icSARS-CoV), GSE45042 (MERS-CoV SA 1), GSE56677 (MERS-CoV Eng 1), GSE33264 (IFN-α, and IFN-γ, Calu-3 treatments), and GSE56678 (IL-1α, and TNF Calu-3 treatments). All data sets and associated metadata have been submitted to Virus Pathogen Resource (ViPR, http://www.viprbrc.org). Study details can be accessed through the Systems Virology website (http://www.systemsvirology.org).