- 1Biosafety Level 3 Core Facility, Yong Loo Lin School of Medicine, National University of Singapore, Singapore, Singapore
- 2Department of Microbiology and Immunology, Yong Loo Lin School of Medicine, National University of Singapore, Singapore, Singapore
- 3Infectious Disease Translational Research Programme, Yong Loo Lin School of Medicine, National University of Singapore, Singapore, Singapore
- 4National Public Health Laboratory, National Centre for Infectious Diseases, Singapore, Singapore
- 5Infectious Disease Research and Training Office, National Centre for Infectious Diseases, Singapore, Singapore
- 6Department of Infectious Diseases, Tan Tock Seng Hospital, Singapore, Singapore
- 7Department of Medicine, Yong Loo Lin School of Medicine, National University of Singapore, Singapore, Singapore
- 8Lee Kong Chian School of Medicine, Nanyang Technological University, Singapore, Singapore
- 9Department of Otolaryngology, Yong Loo Lin School of Medicine, National University of Singapore, Singapore, Singapore
- 10Collaborative and Translation Unit for Hand, Foot and Mouth Disease (HFMD), Institute of Molecular and Cell Biology, Agency for Science, Technology and Research, Singapore, Singapore
The on-going COVID-19 pandemic has given rise to SARS-CoV-2 clades and variants with differing levels of symptoms and severity. To this end, we aim to systematically elucidate the changes in the pathogenesis as SARS-CoV-2 evolved from ancestral to the recent Omicron VOC, on their mechanisms (e.g. cytokine storm) resulting in tissue damage, using the established K18-hACE2 murine model. We reported that among the SARS-CoV-2 viruses tested, infection profiles were initially similar between viruses from early clades but started to differ greatly starting from VOC Delta, where the trend continues in Omicron. VOCs Delta and Omicron both accumulated a significant number of mutations, and when compared to VOCs Alpha, Beta, and earlier predecessors, showed reduced neurotropism and less apparent gene expression in cytokine storm associated pathways. They were shown to leverage on other pathways to cause tissue damage (or lack of in the case of Omicron). Our study highlighted the importance of elucidating the response profiles of individual SARS-CoV-2 iterations, as their propensity of severe infection via pathways like cytokine storm changes as more variant evolves. This will then affect the overall threat assessment of each variant as well as the use of immunomodulatory treatments as management of severe infections of each variant.
Introduction
The emergence of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has resulted in the coronavirus disease 2019 (COVID-19) pandemic that caused major healthcare and economic burden worldwide (1). The infection SARS-CoV-2 caused is highly heterogeneous ranging from asymptomatic and mild to severe and even death (2). In addition, the heterogeneity is further enhanced by the rapid mutation of the virus giving rise to different clades and variants that caused varying pathologies over the course of the pandemic (3). Hence, it is imperative to characterize the differential pathology of SARS-CoV-2 clades and variants and their mechanisms to improve management of the infection, especially for development of therapies that are based on the modulation of these mechanisms.
Coronaviruses have been shown to recombine extensively due to their proofreading mechanism that can facilitate homologous recombination events, giving rise to novel variants of varying pathology and severity (4). To date, SARS-CoV-2 is clustered into 9 major clades based on the Global Initiative on Sharing All Influenza Data (GISAID) database (5, 6). Clade L is the first novel viral genome sequence serving as the reference ancestral strain (7), and other representative clades with major genetic mutations include clades S, V, G, GK, GH, GR, GV, and GRY. Among the clades, clade G and its derivatives have rapidly become the dominant globally circulating SARS-CoV-2 that gave rise to the major variants of concern (VOCs) (8), including Alpha (B.1.1.7), Beta (B.1.351), Gamma (P.1), Delta (B.1.617.2) and Omicron (B.1.1.529). The VOCs have caused concerns as they were shown to have increased virulence, greater transmissibility, resistance to antibody neutralization, and evasion to host- and vaccine-mediated immunity (9–11). It is understood that the SARS-CoV-2, especially the VOCs, resulted in different degrees of severity and symptoms in humans and animal models (12, 13). However, due to the rapid evolution of each VOCs and them rarely overlapping in emergence, there have yet to be a systematic comparison of the early triggers of airway and lung damage among the major SARS-CoV-2 variants, especially between all VOCs.
Several studies have shown that the expression of pro-inflammatory cytokines is associated with disease progression in COVID-19 patients (14–17). Hypercytokinemia, also known as cytokine storm, is one of the major triggers for tissue damage, multiple organ dysfunction, and lethal outcomes in SARS-CoV-2 infection (18, 19). Further, it is well documented that the build-up of activated neutrophils and macrophages in the airways and lungs following SARS-CoV-2 infection results in increased pro-inflammatory cytokines that contribute to airway pathology (20). On the other hand, the evolution of SARS-CoV-2 VOCs over the course of the pandemic showed great variation in their severity and pathology, which is especially apparent in Omicron VOC, suggesting differential cytokine activation that contributes to their pathogenesis (21). Henceforth, we employed the established K18-hACE2 mouse model to systematically compare the disease phenotypes across clades and variants and correlate them with their respective early response cytokine profiles, to determine involvement of cytokine storm and selected pathways’ contribution to their pathogenesis. We aim to elucidate the evolutionary progression of SARS-CoV-2 pathogenesis to understand the pathology spectrum of the virus, which will help in current and future management and drug discovery effort against the virus.
Materials and methods
Cell line and SARS-CoV-2 variants
African green monkey kidney cells, (Vero E6; ATCC CRL-1586™) were cultured in Dulbecco’s Modified Eagle’s Medium (DMEM; Cytiva) supplemented with 10% heat-inactivated fetal calf serum (FCS).
SARS-CoV-2 viruses were isolated from nasopharyngeal swabs of qRT-PCR confirmed COVID-19 patients and propagated in Vero E6 for infection experiments, where low passages of not more than 4 passages of the virus were used. Genome sequences from the swab samples uploaded to GISAID by the National Public Health Laboratory, National Centre for Infectious Diseases, Singapore, were used to confirm the variants’ identity. A total of eight viruses were used: 1. Clade L (ancestral), Lineage B (EPI_ISL_574502); 2. Clade V, Lineage B.29 (EPI_ISL_493419); 3. Clade G, Lineage B.1.610 (EPI_ISL_443240); 4. Clade GR, Lineage B.1.1 (EPI_ISL_443199); 5. VOC Alpha variant, Clade GRY (thereafter GRY-α), Lineage B.1.1.7, (EPI_ISL_754083); 6. VOC Beta variant, Clade GH/501Y.V2 (thereafter as GH-β), Lineage B.1.351.3 (EPI_ISL_1173248); 7. VOC Delta variant, Clade GK (thereafter GK-δ), Lineage B.1.617.2 (EPI_ISL_2621925); 8. VOC Omicron variant, Clade GRA (thereafter GRA-o), Lineage B.1.1.529 (EPI_ISL_7195620). All virus experiments were performed in NUS Medicine biosafety level 3 (BSL-3) core facility and all protocols were approved by the BSL-3 biosafety committee (BBC) and the institutional biosafety committee (IBC) of the National University of Singapore (NUS).
Experimental infection
Eight to nine weeks-old female K18-hACE2 transgenic mice (InVivos Ptd Ltd, Singapore) were acclimatized in the ABSL-3 facility for 72 hours prior to infection and were infected with approximately 1 × 103 PFU of SARS-CoV-2 virus suspension in PBS, via the intranasal route. Baseline body weight were measured prior to virus infection. Body weight, physiological conditions, and survival were monitored daily by two personnel for the duration of the experiment (14-days post infection, dpi) or until the humane endpoint was reached. Each infection group was performed at n=6, except variant GRY-α at n=5, and mock infection at n=3. Scoring of the infected mice physiological conditions were based on five criteria: appearance of mouse coat, level of consciousness, activity level, eye condition, and respiratory quality. Conditions were scored on a scale from 1 to 5, using an observation system adapted from Shrum, et al. (22), with 1 denoting normal physiological state and 5 being the most severe. Additional groups of mice at n=7, except GK-δ at n=5, and GRA-o at n=6 were sacrificed on 4 dpi to assess the viral load in the brain, lung, liver, and spleen; and histological analyses of the lung. Tissues were halved and homogenized in DMEM supplemented with antibiotic-antimycotic and the supernatant was used for plaque assay, while the remaining tissue pellets were kept in -80°C for RNA isolation. The other half of the tissue was fixed with 3.7% formaldehyde (10% formalin) solution for at least 19 h before removal from BSL-3 containment for histological analysis. All animal experiments were conducted in NUS Medicine Biosafety Level 3 (ABSL-3) facility in accordance with NUS IACUC protocol no. R20-0504, using NUS IBC and BBC approved SOPs.
Plaque assay
For virus titre determination, viral supernatants from homogenized tissues were serially diluted in 10-fold increments in DMEM supplemented with antibiotic-antimycotic (Gibco). 250 µl of each serially diluted supernatant were added to confluent Vero E6 cells and incubated for 1 h at 37°C with rocking at 15 min intervals. After 1 h of adsorption, virus inoculum was removed, and cells were washed once with phosphate-buffered saline (PBS). Overlay media containing 1.2% microcrystalline cellulose-DMEM supplemented with antibiotic-antimycotic were added to each well and incubated for 3 days at 37°C at 5% CO2 for plaque formation. Cells were fixed in 10% formalin overnight and stained with crystal violet for plaque visualization. Number of plaques were determined, and virus titre of individual samples were expressed in logarithm of plaque forming units (PFU) per organ.
Histopathological analyses
Formaldehyde fixed tissues were routinely processed, embedded in paraffin blocks (Leica Surgipath Paraplast), sectioned at 5 µm thickness, and stained with hematoxylin and eosin (H&E; Thermo Scientific) following standard histological procedures. A multiparametric, semiquantitative scoring system was further used to assess the magnitude of histo-morphological and -pathological changes in lung tissues based on six criteria: inflammatory cell infiltrates, hemorrhage, edema, degeneration of alveolar epithelial cells, parenchymal wall expansion and bronchiole epithelial cell damage (23, 24). For the histopathological parameters, a score of 0 – 3 was ordinally assigned, where 0 indicated normal; 1 indicated less than 10%; 2 indicated 10 – 50%; and 3 indicated more than 50% of lung regions affected. The average of histopathological scores of the mice from each group was taken as the final evaluation index.
Immunohistochemistry (IHC)
The paraffin embedded lung sections were deparaffinized and rehydrated, followed by heat-mediated antigen retrieval process. The sections were then treated with hydrogen peroxide blocking and protein blocking for 10 min, respectively. After blocking the non-specific binding sites, sections were incubated with anti-SARS-CoV-2 nucleocapsid protein monoclonal antibody (1:1000; Abcam), followed by the respective horseradish peroxidase (HRP)-conjugated secondary antibody. The sections were subsequently visualized using DAB solution and counterstained with hematoxylin.
RNA isolation and qRT-PCR
Homogenized lung tissue pellets were subjected to total RNA extraction using QIAGEN RNeasy® mini kit in accordance with manufacturer’s guidelines. RNA purity and concentration were determined using Tecan Infinite® M200 Pro coupled with NanoQuant Plate™. Real-time reverse transcription polymerase chain reaction (qRT-PCR) was performed using QIAGEN RT2 First Strand Kit, QIAGEN RT2 SYBR® Green Mastermix, and QIAGEN RT2 Profiler PCR Array (PAMM-150ZE-4) in accordance with manufacturer’s recommendations.
Gene expression analysis
Murine cytokines mRNA expression levels from cDNA conversion were analyzed using QIAGEN GeneGlobe Design and Analysis hub, and gene expression were expressed as Log2 fold regulations. Expression levels were calculated relative to Gapdh and normalized to the mock infected mice using the ΔΔCT method with the QIAGEN GeneGlobe software. The genes from the QIAGEN panel were further sub-divided into selected pathways that may potentially contribute to tissue damage using DAVID functional annotation tool (25, 26). Pathways were selected from Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases for visualization of enrichment in each pathway across variants.
Statistical analysis
Analyses were calculated using two-tailed Mann-Whitney U test to evaluate the data obtained using Graphpad Prism version 9.0 with * denoting that p < 0.05, ** denoting that p < 0.01, and *** denoting p < 0.001. For gene expression analysis, p-values were obtained from the analysis with QIAGEN GeneGlobe Design and analysis hub. Only genes fulfilling both criteria of >2 fold-regulation compared to uninfected controls and having a p-value of <0.05 were considered significant.
Results
To characterize the pathogenesis of various SARS-CoV-2 viruses from different clades and VOCs, and their respective mechanisms leading to airway and lung damage, we tracked a total of eight major SARS-CoV-2 iterations – L (ancestral), V, G, GR, GRY-α, GH-β, GK-δ, and GRA-o (Figure 1). The repertoire of our coverage includes the major SARS-CoV-2 clades and variants of the pandemic to date. The mutations across the iterations based on their clinical isolates’ sequences obtained from GISAID were listed in Table 1.
Figure 1 Summarized methodology of characterization of eight viral iterations (ancestral to Omicron) of SARS-CoV-2. K18-hACE2 mice were infected the eight different SARS-CoV-2 viruses and divided into survival group and 4 dpi group (sacrificed 4 days post infection). The survival group were used for assessment of survival and physiological assessment. The 4 dpi group were sacrificed and have organs harvested for virologic, histologic and molecular profiling analysis (cytokine and chemokine gene expressions). Figure created with BioRender.com.
Changes in survival, symptoms, and tissue tropism over the course of SARS-CoV-2 evolutionary iterations
We first compared the survival of K18-hACE2 mice following infection with the different SARS-CoV-2 viruses. It was observed that, with the exception of V (34% survival) and GRA-o (83.3% survival), 0% of infected mice survived the infection, where they succumbed between 6 to 8 dpi, with GH-β and GK-δ hitting the end point the earliest at 6 dpi (Figure 2A). When comparing the Kaplan-Meier survival curves, it was shown that only GRA-o, the latest evolutionary iteration, was significantly different in mice survival compared with all other infected groups except for V, another iteration of the virus where there were mice that survived (Figure 2B). The weight loss trends largely followed the survival curves in which rapid weight loss were observed between 4 and 7 dpi for all groups other than V and GRA-o, for which we observed a rebound in weight from surviving mice (Figure 2C). We then compared clinical symptom scores (averaged from appearance, eye condition, consciousness, activity, and respiration) and observed that L and GRY-α showed the most symptom presentation, followed by the others that showed moderate symptoms, while GRA-o was shown to be not presenting any clinical symptom throughout the duration of study (Figure 2D, Figures S1A–E). We then examined viral titre retrieved from the tissues (lung, brain, liver, and spleen) of SARS-CoV-2 infected mice. No infectious viruses were recovered from liver and spleen homogenates using plaque assay (data not shown). For viral titre in the lungs, high titres of SARS-CoV-2 was retrieved from the lungs of L, V, G, GR and GRY-α infected mice, at geometric mean titres of 4.27 × 105, 4.55 × 104, 2.75 × 105, 1.09 × 106, and 3.12 × 105 PFU/organ, respectively. Titres from V are slightly lower, while GR slightly higher when compared with L infected group (p < 0.05). On the other hand, significantly lower titre was observed in GH-β, GK-δ and GRA-o when compared with L, at 6.86 × 102, 7.2 × 101 and 1.1 × 101 PFU/organ, respectively (p < 0.05) (Figure 2E). As for viral titre in the brain, comparatively high titre was found in the brain tissue of L, GRY-α and GH-β infected groups, at 1.93 × 105, 6.02 × 103 and 3.4 × 103 PFU/organ, respectively; while significantly lower titre in the brain was retrieved in brain tissue of G, GR and GK-δ at 1.03 × 102, 4.4 × 101, and 8.7 × 101 PFU/organ, respectively. No virus was retrieved from the brain of GRA-o infected mice (Figure 2F). Notably, unlike in the lung, viral titre in the brain was more variable where about half of the mice did not have viral titre retrieved from the brain in most groups, regardless of the titre observed.
Figure 2 SARS-CoV-2 viruses infected K18-hACE2 transgenic mice. Physiological conditions and weight changes of infected mice group across different SARS-CoV-2 viruses (L, V, G, GR, GRY-α, GH-β, GK-δ, and GRA-o), and uninfected mice over 14 dpi. Mice were infected with approximately 103 PFU of the respective SARS-2-CoV viruses through intranasal delivery. (A) Survival percentage of infected mice over 14 days, (B) Significance comparison of survival curves for SARS-CoV-2 infected mice. Comparison of significance between survival curves of various variant infected mice using log-rank (Mantel-Cox) test. Cut off for p-value using Kaplan-Meier analysis was adjusted to 0.002 for multiple comparison between various survival curves. (C) percentage changes in body weight, (D) changes in physiological condition quantified through scoring of clinical symptoms. Individual mouse was scored across 14 days from 1 (normal) to 5 (worst) for five conditions: coat appearance, level of consciousness, activity, eye condition, and respiration. To determine clinical score of the infected mice, average of the five physiological parameters for each infection group. All infection groups, n=6 with except of variant G, n=5, and mock infected mice, n=3. Error bars denotes range for each time point. (E, F) For viral load comparison after infection, mice were sacrificed at 4 dpi (L, V, G, GR, GRY-α, GH-β: n = 7; GK-δ: n = 5; GRA-o: n = 6) with lungs, brain, liver, and spleen tissues harvested for titration. SARS-CoV-2 virus titre in lung and brain tissue determined via plaque assays. Geometric mean of viral titre is shown. Statistical significance between clade L and variants determined with two-tailed Mann-Whitney U test. *p < 0.05, **p < 0.01 and ***p < 0.001.
Early histopathological damage of SARS-CoV-2 viruses were not clearly distinguishable
From the lung tissue collected at 4 dpi, it was observed that histopathological damage of SARS-CoV-2 variants was largely similar, with histopathological index score ranging between 0.63 to 1.69. Tissues from GR and GK-δ infected mice were on the higher end of the scores while those from L and GRA-o were with lower scores (Figure 3A). The factors that likely contributed to the higher histopathological index were inflammatory cell infiltrates, where the trend was consistent with the overall index compared with other factors (Figure 3B; Figures S2A–F). On the other hand, no histological abnormality was observed in the bronchiole epithelial cells following all variants’ inoculation, and neither bronchitis nor bronchiole epithelial denudation was present in the lung tissues (data not shown). Next, we correlated the extent of infected cells and viral presence in the lungs of the ancestral L infected group, with the VOC infected groups. Using IHC of SARS-CoV-2 nucleocapsid protein, we observed that the viral distribution in the lungs at 4 dpi varied greatly (Figures 3C, D). GK-δ, notably, showed major nucleocapsid protein presence in the lung (93.2% coverage), contrary to the low infectious titre detected from the infected murine lungs while GRA-o had almost no nucleocapsid protein detected in the lung tissues (5.83% coverage). Our observation suggested that there was no correlation between nucleocapsid coverage, which are produced in excess during coronavirus infection (27, 28), in the lungs and the lung histopathological damage index between variants.
Figure 3 Histopathological and immunohistochemical analyses of lung sections from K18-hACE2 mice following infection with SARS-CoV-2 at 4 dpi. (A) Histopathological scoring of lung damage by evaluating the extent and areas of pulmonary involvement based on 6 criteria: inflammatory cell infiltrates, hemorrhage, edema, degeneration of alveolar epithelial cells, parenchymal wall expansion and bronchiole epithelial cell damage, assigned respectively with a score of 0 – 3. The average of the 6 parameters from mice for each group was regarded as the final histopathological scoring. Mice following mock infection are representative of n=3 per group while mice infected with SARS-CoV-2 are representative of n=7 per group (n=5 for GK-δ, and n=6 for GRA-o). (B) Hematoxylin and eosin staining of lung sections from mock- or SARS-CoV-2-infected K18-hACE2 mice. H&E-stained lung sections demonstrated massive inflammatory cell infiltrates, hemorrhage and thickened interalveolar septa in the mice following SARS-CoV-2 inoculation. Hemorrhage was observed in the mock-infected mice could be due to a crushing artifact and damaged capillary walls during tissue handling at necropsy. 50x total magnification, scale bars = 50µm. (C) Percentage of SARS-CoV-2 nucleocapsid protein positivity in the FFPE lung tissues infected with L, GRY-α, GH-β, GK-δ and GRA-o. (D) Detection of SARS-CoV-2 nucleocapsid protein (brown staining) showed the virus was predominantly found in alveolar septa and pneumocytes, but not bronchiole regions. DAB chromogen and hematoxylin counterstain were used. Each image is representative of n=7 per group (L, V, G, GR, GRY-α, GH-β), n=5 (GK-δ), n=6 (GRA-o), n=3 (mock). 5x total magnification, scale bars = 2000µm.
SARS-CoV-2 cytokine profiles between evolutionary iterations revealed differences in lung damage triggers
To further elucidate potential mechanisms in each of the SARS-CoV-2 viruses that led to their respective phenotypes, we performed a qRT-PCR array on 84 genes related to inflammatory responses from lung tissues of SARS-CoV-2 infected K18-hACE2 mice. After taking into account fold-regulation and statistical significance (>2 fold-regulation and p< 0.05), the cytokine profiles between the infected groups differed, some significantly, from each other. We noted the groups with the larger significant changes from the ancestral L virus (18 up; 0 down) were G (23 up; 6 down), GR (21 up; 2 down), GRY-α (26 up; 2 down), GH-β (21 up; 5 down) and GRA-o (26 up; 1 down). On the other hand, V (4 up; 0 down) and GK-δ (6 up; 4 down) had fewer significant changes compared with L (Figure 4). The exact changes in fold-regulation of the genes tested were shown in Table S1. In general, the host responses of the SARS-CoV-2 infection of all groups tested deviated significantly from the ancestral L virus, with some having unique significant responses as observed in G (Lta, Tnf, Pf4), GH-β (Il24, Il1a, Mstn, Il22), GK-δ (Il7, Tnfsf11b, Il18), and GRA-o (Ccl22, Bmp4, Tnfsf10, Mif, Hc, Bmp7, Gpi1, Cx3cl1). In view of the differences present, particularly in GRA-o which has the most of them, further dissection on how the differential cytokine changes led to the pathogenesis of current and future variants is thus warranted.
Figure 4 Comparison of 84 cytokines regulation changes among SARS-CoV-2 viruses at 4 dpi. Differentially regulated cytokines across all SARS-CoV-2 virus tested. Color scale indicating level of averaged fold-regulation changes across the groups, * indicate significant difference in fold-regulation changes (>2-fold, p<0.05) between infected and mock infected mice.
As cytokine storm was perceived to be a major trigger of lung damage in COVID-19 (18, 19), we grouped the panel of genes we tested into pathways associated with cytokine storm using DAVID. We selected pathways involved directly in cytokine storm (19), and associated pathways like MAPK cascade, TNF responses, neutrophil chemotaxis, IL-17 regulation, negative regulation of inflammatory response, and T-cell mediated cytotoxicity to determine how much they played a role in tissue damage for each virus iteration (Figure 5). Firstly, in the cytokine storm pathway, we observed that for iterations that come after L and V, more genes such as Ifng, Il12b and Il10 were found to be significantly differentially expressed in G, GR, GRY-α and GH-β, before being reduced again in GK-δ and GRA-o (Figure 5A). Similar trends were observed in TNF response and neutrophil chemotaxis pathways in genes like Xcl1, Ccl1 and Ccl11 (Figures 5B, C). The MAPK cascade pathway was prominently differentially expressed in G, GR, GRY-α and GH-β, but in this case, L and GK-δ as well (Figure 5D). On the other hand, negative regulation of inflammatory responses early in infection were linked with cytokine storm (29), and we found that G, GR, GRY-α and GH-β further have differentially up-regulated expression of Il10 (Figure 5E). IL-17 pathway, associated with neutrophil infiltration with genes from the pathway like Il15 and Il12b were up-regulated in all groups except V and GK-δ (Figure 5F). Finally, T-cell mediated cytotoxicity pathway was found to be prominently changed for G, GR, GRY-α, GH-β, and GRA-o (Figure 5G). Interestingly, only GRA-o had a further significant up-regulation of Il12a. Overall, genes from cytokine storm and associated pathways (Ccl2, Ccl7 and Ccl12) were found to be differentially regulated in infections of L, G, GR, GRY-α, and GH-β while not as strongly involved in V and GK-δ. GRA-o, while having similar trend of cytokine storm associated pathway gene changes, had more subdued expression, and a slight difference in T-cell mediated cytotoxicity (Il12a).
Figure 5 Comparison of differentially expressed genes in selected pathways relevant to cytokine storm and tissue damage at 4 dpi. Differentially regulated cytokines across all groups in pathways related to (A) cytokine storm; (B) cellular response to TNF (C) neutrophil chemotaxis (D) MAPK cascade (E) negative regulation of inflammatory responses (F) positive regulation of IL-17 production (G) positive regulation of T-cell mediated cytotoxicity. Color scale indicate level of averaged fold-regulation changes across the groups, * indicate significant difference in fold-regulation changes (>2-fold, p<0.05) between infected and mock infected mice.
Discussion
The use of K18-hACE2 murine model to understand SARS-CoV-2 infection has been instrumental to the understanding of SARS-CoV-2 pathogenesis mechanisms since the start of the pandemic, including comparisons of early clades’ differential severity (13, 30–34). Moreover, the K18-hACE2 mice were shown to be simulating severe SARS-CoV-2 infection useful in the identification of immune responses leading to severe outcome (35). Hence, here we reported results using the model to correlate early cytokine profiles (4 dpi) of early SARS-CoV-2 evolutionary iterations (L, V, G, GR), and the VOCs Alpha (GRY-α), Beta (GH-β), Delta (GK-δ) and Omicron (GRA-o), with their infection outcome. Interestingly, the majority of the pathogenesis mechanism and infection outcome remained largely similar up until the Beta VOC, before major changes occur in VOCs Delta and Omicron that emerged later. This may be partially due to the larger number of mutations accumulated (Table 1) for Delta and Omicron VOC, which potentially altered their biology significantly.
Our study demonstrated that early cytokine profiles, viral titre and histopathological damage worked in concert to influence and determine the infection outcome. Overall, the infection outcome remained largely similar from L to GH-β Beta VOC, and then deviated in VOCs that emerged late in the pandemic i.e. GK-δ Delta and GRA-o Omicron VOCs. Among the early clades and VOCs up to Beta, even the milder clade V virus with ~30% survival rate had similar rates of neurotropism based on viral titre from brain tissue, sufficiently high viral titre in the lungs, and cytokine profiles consistent with previous report (36). While clade V virus had fewer significantly differentially expressed genes, the major genes involved in cytokine storm; TNF pathway and neutrophil chemotaxis, Ccl2 and Ccl7 remained, suggesting similar, but milder presentation of pathology for its infection outcome. The subsequent clade and VOCs in the earlier part of the pandemic, G, GR, GRY-α and GH-β had further inflammatory suppression gene Il10 differentially up-regulated that potentially augmented their pathology slightly, as seen in the histopathological index. Furthermore, the viruses from earlier in the pandemic, up to VOC GH-β, showed higher clinical presentation scores, which may be linked to neurotropism, prominent in these groups of infections.
On the other hand, VOCs that emerged later in the pandemic, GK-δ Delta and GRA-o Omicron, followed a different pathogenesis mechanism that resulted in tissue damage. Firstly, neurotropism was observed to be less apparent in Delta, and completely absent in the Omicron variant in the K18-hACE2 mice. It was observed that mechanisms that may contribute to cytokine storm were also less apparent in GK-δ and GRA-o infection. In addition, immune suppression, a mechanism common in SARS-CoV-2 (37), was observed strongly in the GK-δ Delta variant, albeit different from Il10 mediated suppression observed in clades and VOCs that came before. The Delta variant in our study instead showed an almost universal suppression of host responses early post infection at 4dpi (with very little significant differential regulation). This suppression may potentially contribute to the increased nucleocapsid expression and inflammatory infiltrates in the lungs, which may suggest diffusely infected cells in the lungs. On the other hand, the GRA-o variant showed the most varied response among the variants where we observed the least tissue damage, low viral titre and lack of neurotropism. This observation is in line with the finding that Omicron variant has more preference to the upper airway than the lungs (38). In addition, the cytokine profile of our Omicron GRA-o variant infection showed a response largely different than that of its predecessors. It had generally weaker changes in expression of cytokine storm and related pathway genes, coupled with moderately stronger T-cell mediated cytotoxicity gene expression in Il12a and Il12b. The Omicron cytokine profile suggests a different pathogenesis that may contribute to the overall milder tissue damage congruent with other studies to date (21, 39). The effects of the major changes in Omicron cytokine profiles require further studies to validate their correlation with disease outcome, as well as potential interaction with risk factors and chronic conditions (e.g. asthma) due to their major differences and lower lethality.
The finding from our study is crucial in the understanding of SARS-CoV-2 pathogenesis mechanisms. It showed the evolutionary potential of SARS-CoV-2 where its ability to accumulate high number of mutations may result in highly variable pathogenesis and infection severity in the variants tested and may further change in future variants (4, 40). This is particularly evident in our elucidation of Delta and Omicron cytokine profiles suggesting a major change in the biology of the VOCs late in the pandemic (41, 42). Indeed, the emergence of Delta and Omicron variants both caused major outbreaks that dwarfed those that of their predecessors (43, 44). Studies have shown that Delta and Omicron VOCs has evolved into their distinct evolutionary group and therefore may explain the significant differences of their infection profiles and pathogenesis mechanisms (41). In addition, other reasons such as gestation and adaptation in immunocompromised patients (45), spill-over and spill-back to and back from animal reservoirs (45–47), and potential vaccine induced selection pressure (48), may all contribute to the accumulation of these mutations and the differential pathogenesis mechanisms. Overall, our study inferred that the adaptation of the virus may likely lead to emergence of future variants having mutations that drives their pathogenesis differently than existing SARS-CoV-2 viruses. Therefore, this implied that management strategies based on immunomodulation and mediation of host responses have to be revisited in different variants, mirroring the management of influenza, where propensity of cytokine storm differs (49). Finally, our finding may also partially explain the continued efficacy of broad spectrum immunomodulators like corticosteroids; but not more targeted therapies targeting specific pathways due to the potential differences in different SARS-CoV-2 iterations’ pathogenesis mechanisms.
Our study, however, is not without its limitation. We only elucidated the histopathology and cytokine profile at 4 dpi, that only gave us insights of host response early during infection. This also created potential disconnect such as the case between the high Delta nucleocapsid detection in IHC compared to its low lung infectious viral titre. Nevertheless, high viral protein, especially in the case of nucleocapsid, which were produced in excess during infection (28), may not necessarily translate to infectious viral titre. In addition, we only use female mice for our comparative study and thus may not capture any potential sex specific differences between the cytokine profiles of different virus iterations (50). However, a previous study showed that cytokine and chemokine responses between male and female were largely similar, where the major differences lie at the magnitude and timing of the responses (51). Therefore, our study remained representative of the differences in cytokine profiles between variants despite only testing it in female mice. Nevertheless, future studies with male mice can be performed to investigate the extent of differential cytokine expression levels that the different SARS-CoV-2 virus iteration can induce. Finally, we also did not manage to assess the protein levels of the cytokines in our study. Despite the limitations however, our data showed that, by systematically comparing the profiles across the major evolutionary iterations of SARS-CoV-2, a clear heterogeneity of SARS-CoV-2 evolution and adaptation was observed, especially in VOCs Delta and Omicron that emerged later
In conclusion, our study provided insights on the evolution of pathogenesis and its potential mechanisms of SARS-CoV-2 throughout the course of the pandemic. We have observed that adaptation throughout the pandemic can give rise to variants of very different pathogenesis mechanisms that contribute to tissue damage. Our study thus emphasized the importance in differentiating host responses to different SARS-CoV-2 variants. This is especially true as we enter the endemic phase, where the surveillance of differential host responses between variants will become increasingly crucial when assessing future variants, for their threat and impact, as well as when devising immunomodulatory treatments for severe infections.
Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.
Ethics statement
The animal study was reviewed and approved by NUS Institutional Animal Care and Use Committee (IACUC) under protocol no. R20-0504.
Author contributions
CKM and JC conceived and designed the experiments. CKM, ZA, TM, DL, and RL contributed to the virus isolation and sequencing. CKM, ZA, HC, and YW performed the experiment. CKM, YW, and KT performed the histological analyses. ZA, CKM, and KT performed the cytokine analyses. CKM, ZA, YW, KT, and JC performed the data analyses. CKM, ZA, YW, KT, and JC wrote the manuscript. All authors have read and approved the final version of the manuscript prior to submission.
Funding
This research was supported by the following grants: NUHSRO/2020/066/NUSMedCovid/01/BSL3 Covid Research Work, NUHSRO/2020/050/RO5+5/NUHS-COVID/4, Ministry of Education, Singapore MOE2017-T2-2-014, Singapore NMRC Centre Grant Program – Diabetes, Tuberculosis and Neuroscience CGAug16M009, Ministry of Health MOH-COVID19RF2-0001.
Acknowledgments
We are grateful to the National University of Singapore, Yong Loo Lin School of Medicine BSL-3 Core Facility for their support with this work.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2022.950666/full#supplementary-material
Supplementary Table 1 | Murine lung mRNA fold regulation changes of 84 selected cytokines and chemokines induced by different SARS-CoV-2 infections. Expression changes in Log2 fold-regulation and p-values of 84 mRNA isolated from infected mice lung tissues at 4 dpi separated into L, V, G, GR, GRY-α, GH-β, GK-δ, and GRA-o respectively, when compared to uninfected control.
Supplementary Figure 1 | Physiological and symptomatic comparison for SARS-CoV-2 infected mice. Individual physiological condition changes of all eight SARS-CoV-2 infected mice groups over 4 to 6 dpi. Breakdown of individual parameter changes over time from 4 to 6 dpi for (A) appearance, (B) eye condition, (C) level of consciousness, (D) activity, and (E) respiration respectively.
Supplementary Figure 2 | Histopathological comparison of the extent and areas of lung damage inflicted with different SARS-CoV-2 variants. Blinded scoring and quantification of lung injury in respect of (A) inflammatory cell infiltrates (B) hemorrhage (C) pulmonary edema (D) degeneration of alveolar epithelial cells, and (E) parenchymal wall expansion. Data are presented as average of the histopathological scoring from mice for each group. (F) Hematoxylin and eosin-stained images correspond to the mice following mock infection or infection of SARS-CoV-2 (L, GRY-α, GH-β, GK-δ and GRA-o). Images show low-power (top; 5x total magnification; scale bars = 2000µm) and medium-power magnification (bottom; 50x total magnification; scale bars = 50µm). SARS-CoV-2-infected mice are representative of n=7 per group (variants L, V, G, GR, GRY-α, GH-β), n=5 (GK-δ), and n=6 (GRA-o) while mock-infected mice are representative of n=3.
References
1. Hu B, Guo H, Zhou P, Shi ZL. Characteristics of SARS-CoV-2 and COVID-19. Nat Rev Microbiol (2021) 19:141–54. doi: 10.1038/s41579-020-00459-7
2. Mizrahi B, Shilo S, Rossman H, Kalkstein N, Marcus K, Barer Y, et al. Longitudinal symptom dynamics of COVID-19 infection. Nat Commun (2020) 11:6208. doi: 10.1038/s41467-020-20053-y
3. Tao K, Tzou PL, Nouhin J, Gupta RK, de Oliveira T, Kosakovsky Pond SL, et al. The biological and clinical significance of emerging SARS-CoV-2 variants. Nat Rev Genet (2021) 22:757–73. doi: 10.1038/s41576-021-00408-x
4. Gribble J, Stevens LJ, Agostini ML, Anderson-Daniels J, Chappell JD, Lu X, et al. The coronavirus proofreading exoribonuclease mediates extensive viral recombination. PloS Pathog (2021) 17:e1009226. doi: 10.1371/journal.ppat.1009226
5. Maxmen A. One million coronavirus sequences: Popular genome site hits mega milestone. Nature (2021) 589:21. doi: 10.1038/d41586-021-01069-w
6. Elbe S, Buckland-Merrett G. Data, disease and diplomacy: GISAID’s innovative contribution to global health. Glob Chall (2017) 1:33–46. doi: 10.1002/gch2.1018
7. Wu F, Zhao S, Yu B, Chen YM, Wang W, Song ZG, et al. A new coronavirus associated with human respiratory disease in China. Nature (2020) 579:265–9. doi: 10.1038/s41586-020-2008-3
8. World Health, O. COVID-19 weekly epidemiological update. edition 68. Geneva: World Health Organization (2021).
9. Yang TJ, Yu PY, Chang YC, Liang KH, Tso HC, Ho MR, et al. Effect of SARS-CoV-2 B.1.1.7 mutations on spike protein structure and function. Nat Struct Mol Biol (2021) 28:731–9. doi: 10.1038/s41594-021-00652-z
10. Tegally H, Wilkinson E, Giovanetti M, Iranzadeh A, Fonseca V, Giandhari J, et al. Detection of a SARS-CoV-2 variant of concern in south Africa. Nature (2021) 592:438–43. doi: 10.1038/s41586-021-03402-9
11. Wang P, Nair MS, Liu L, Iketani S, Luo Y, Guo Y, et al. Antibody resistance of SARS-CoV-2 variants B.1.351 and B.1.1.7. Nature (2021) 593:130–5. doi: 10.1038/s41586-021-03398-2
12. Boehm E, Kronig I, Neher RA, Eckerle I, Vetter P, Kaiser L, et al. Novel SARS-CoV-2 variants: the pandemics within the pandemic. Clin Microbiol Infect (2021) 27:1109–17. doi: 10.1016/j.cmi.2021.05.022
13. Lee TY, Lee H, Kim N, Jeon P, Kim JW, Lim HY, et al. Comparison of SARS-CoV-2 variant lethality in human angiotensin-converting enzyme 2 transgenic mice. Virus Res (2021) 305:198563. doi: 10.1016/j.virusres.2021.198563
14. Su CM, Wang L, Yoo D. Activation of NF-kappaB and induction of proinflammatory cytokine expressions mediated by ORF7a protein of SARS-CoV-2. Sci Rep (2021) 11:13464. doi: 10.1038/s41598-021-92941-2
15. Chen G, Wu D, Guo W, Cao Y, Huang D, Wang H, et al. Clinical and immunological features of severe and moderate coronavirus disease 2019. J Clin Invest (2020) 130:2620–9. doi: 10.1172/JCI137244
16. Huang C, Wang Y, Li X, Ren L, Zhao J, Hu Y, et al. Clinical features of patients infected with 2019 novel coronavirus in wuhan, China. Lancet (2020) 395:497–506. doi: 10.1016/s0140-6736(20)30183-5
17. Zhou Z, Ren L, Zhang L, Zhong J, Xiao Y, Jia Z, et al. Heightened innate immune responses in the respiratory tract of COVID-19 patients. Cell Host Microbe (2020) 27:883–890 e882. doi: 10.1016/j.chom.2020.04.017
18. Coperchini F, Chiovato L, Ricci G, Croce L, Magri F, Rotondi M, et al. The cytokine storm in COVID-19: Further advances in our understanding the role of specific chemokines involved. Cytokine Growth Factor Rev (2021) 58:82–91. doi: 10.1016/j.cytogfr.2020.12.005
19. Yang L, Xie X, Tu Z, Fu J, Xu D, Zhou Y, et al. The signal pathways and treatment of cytokine storm in COVID-19. Signal Transduct Target Ther (2021) 6:255. doi: 10.1038/s41392-021-00679-0
20. Channappanavar R, Fehr AR, Vijay R, Mack M, Zhao J, Meyerholz DK, et al. Dysregulated type I interferon and inflammatory monocyte-macrophage responses cause lethal pneumonia in SARS-CoV-Infected mice. Cell Host Microbe (2016) 19:181–93. doi: 10.1016/j.chom.2016.01.007
21. Halfmann PJ, Iida S, Iwatsuki-Horimoto K, Maemura T, Kiso M, Scheaffer SM, et al. SARS-CoV-2 omicron virus causes attenuated disease in mice and hamsters. Nature (2022) 603:687–92. doi: 10.1038/s41586-022-04441-6
22. Shrum B, Anantha RV, Xu SX, Donnelly M, Haeryfar SM, McCormick JK, et al. A robust scoring system to evaluate sepsis severity in an animal model. BMC Res Notes (2014) 7:1–11. doi: 10.1186/1756-0500-7-233
23. Li F, Han M, Dai P, Xu W, He J, Tao X, et al. Distinct mechanisms for TMPRSS2 expression explain organ-specific inhibition of SARS-CoV-2 infection by enzalutamide. Nat Commun (2021) 12:866. doi: 10.1038/s41467-021-21171-x
24. Leist SR, Dinnon KH 3rd, Schafer A, Tse LV, Okuda K, Hou YJ, et al. A mouse-adapted SARS-CoV-2 induces acute lung injury and mortality in standard laboratory mice. Cell (2020) 183:1070–85 e1012. doi: 10.1016/j.cell.2020.09.050
25. Huang da W, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc (2009) 4:44–57. doi: 10.1038/nprot.2008.211
26. Sherman BT, Hao M, Qiu J, Jiao X, Baseler MW, Lane HC, et al. DAVID: a web server for functional enrichment analysis and functional annotation of gene lists (2021 update). Nucl Acids Res (2022) 50:W216–21. doi: 10.1093/nar/gkac194
27. Surjit M, Lal SK. The SARS-CoV nucleocapsid protein: a protein with multifarious activities. Infect Genet Evol (2008) 8:397–405. doi: 10.1016/j.meegid.2007.07.004
28. Savastano A, Ibanez de Opakua A, Rankovic M, Zweckstetter M. Nucleocapsid protein of SARS-CoV-2 phase separates into RNA-rich polymerase-containing condensates. Nat Commun (2020) 11:6041. doi: 10.1038/s41467-020-19843-1
29. Lu L, Zhang H, Dauphars DJ. & he, y. w. a potential role of interleukin 10 in COVID-19 pathogenesis. Trends Immunol (2021) 42:3–5. doi: 10.1016/j.it.2020.10.012
30. Oladunni FS, Park JG, Pino PA, Gonzalez O, Akhter A, Allue-Guardia A, et al. Lethality of SARS-CoV-2 infection in K18 human angiotensin-converting enzyme 2 transgenic mice. Nat Commun (2020) 11:6122. doi: 10.1038/s41467-020-19891-7
31. Yinda CK, Port JR, Bushmaker T, Offei Owusu I, Purushotham JN, Avanzato VA, et al. K18-hACE2 mice develop respiratory disease resembling severe COVID-19. PloS Pathog (2021) 17:e1009195. doi: 10.1371/journal.ppat.1009195
32. Zheng J, Wong LR, Li K, Verma AK, Ortiz ME, Wohlford-Lenane C, et al. COVID-19 treatments and pathogenesis including anosmia in K18-hACE2 mice. Nature (2021) 589:603–7. doi: 10.1038/s41586-020-2943-z
33. Sun SH, Chen Q, Gu HJ, Yang G, Wang YX, Huang XY, et al. A mouse model of SARS-CoV-2 infection and pathogenesis. Cell Host Microbe (2020) 28:124–133 e124. doi: 10.1016/j.chom.2020.05.020
34. Dong W, Mead H, Tian L, Park JG, Garcia JI, Jaramillo S, et al. The K18-human ACE2 transgenic mouse model recapitulates non-severe and severe COVID-19 in response to an infectious dose of the SARS-CoV-2 virus. J Virol (2022) 96:e0096421. doi: 10.1128/JVI.00964-21
35. Yang S, Cao L, Xu W, Xu T, Zheng B, Ji Y, et al. Comparison of model-specific histopathology in mouse models of COVID-19. J Med Virol (2022) 94:3605–12. doi: 10.1002/jmv.27747
36. McCray PB Jr., Pewe L, Wohlford-Lenane C, Hickey M, Manzel L, Shi L, et al. Lethal infection of K18-hACE2 mice infected with severe acute respiratory syndrome coronavirus. J Virol (2007) 81:813–21. doi: 10.1128/JVI.02012-06
37. Gamage AM, Tan KS, Chan WOY, Liu J, Tan CW, Ong YK, et al. Infection of human nasal epithelial cells with SARS-CoV-2 and a 382-nt deletion isolate lacking ORF8 reveals similar viral kinetics and host transcriptional profiles. PloS Pathog (2020) 16:e1009130. doi: 10.1371/journal.ppat.1009130
38. Hui KPY, Ho JCW, Cheung MC, Ng KC, Ching RHH, Lai KL, et al. SARS-CoV-2 omicron variant replication in human bronchus and lung ex vivo. Nature (2022) 603:715–20. doi: 10.1038/s41586-022-04479-6
39. Shuai H, Chan JF, Hu B, Chai Y, Yuen TT, Yin F, et al. Attenuated replication and pathogenicity of SARS-CoV-2 B.1.1.529 omicron. Nature (2022) 603:693–9. doi: 10.1038/s41586-022-04442-5
40. Plante JA, Mitchell BM, Plante KS, Debbink K, Weaver SC, Menachery VD, et al. The variant gambit: COVID-19’s next move. Cell Host Microbe (2021) 29:508–15. doi: 10.1016/j.chom.2021.02.020
41. Bansal K, Kumar S. Mutational cascade of SARS-CoV-2 leading to evolution and emergence of omicron variant. Virus Res (2022) 315:198765. doi: 10.1016/j.virusres.2022.198765
42. Zhan Y, Yin H. & yin, j. y. B.1.617.2 (Delta) variant of SARS-CoV-2: features, transmission and potential strategies. Int J Biol Sci (2022) 18:1844–51. doi: 10.7150/ijbs.66881
43. Tareq AM, Emran TB, Dhama K, Dhawan M, Tallei TE. Impact of SARS-CoV-2 delta variant (B.1.617.2) in surging second wave of COVID-19 and efficacy of vaccines in tackling the ongoing pandemic. Hum Vaccin Immunother (2021) 17:4126–7. doi: 10.1080/21645515.2021.1963601
44. Paz M, Aldunate F, Arce R, Ferreiro I, Cristina J. An evolutionary insight into severe acute respiratory syndrome coronavirus 2 omicron variant of concern. Virus Res (2022) 314:198753. doi: 10.1016/j.virusres.2022.198753
45. Kupferschmidt K. Where did ‘weird’ omicron come from? Science (2021) 374:1179. doi: 10.1126/science.acx9738
46. Islam A, Ferdous J, Islam S, Sayeed MA, Dutta Choudhury S, Saha O, et al. Evolutionary dynamics and epidemiology of endemic and emerging coronaviruses in humans, domestic animals, and wildlife. Viruses (2021) 13:1908. doi: 10.3390/v13101908
47. Chen Q, Huang XY, Liu Y, Sun MX, Ji B, Zhou C, et al. Comparative characterization of SARS-CoV-2 variants of concern and mouse-adapted strains in mice. J Med Virol (2022) 94:3223–32. doi: 10.1002/jmv.27735
48. Koyama T, Miyakawa K, Tokumasu RSSJ, Kudo M, Ryo A. Evasion of vaccine-induced humoral immunity by emerging sub-variants of SARS-CoV-2. Future Microbiol (2022) 17:417–24. doi: 10.2217/fmb-2022-0025
49. Ryabkova VA, Churilov LP, Shoenfeld Y. Influenza infection, SARS, MERS and COVID-19: Cytokine storm - the common denominator and the lessons to be learned. Clin Immunol (2020) 223:108652. doi: 10.1016/j.clim.2020.108652
50. Peckham H, de Gruijter NM, Raine C, Radziszewska A, Ciurtin C, Wedderburn LR, et al. Male Sex identified by global COVID-19 meta-analysis as a risk factor for death and ITU admission. Nat Commun (2020) 11:6317. doi: 10.1038/s41467-020-19741-6
Keywords: SARS-CoV-2, variants of concern, immune response, cytokine storm, K18-hACE2 mice model
Citation: Aw ZQ, Mok CK, Wong YH, Chen H, Mak TM, Lin RTP, Lye DC, Tan KS and Chu JJH (2022) Early pathogenesis profiles across SARS-CoV-2 variants in K18-hACE2 mice revealed differential triggers of lung damages. Front. Immunol. 13:950666. doi: 10.3389/fimmu.2022.950666
Received: 23 May 2022; Accepted: 04 October 2022;
Published: 27 October 2022.
Edited by:
Kari Ann Shirey, University of Maryland, United StatesReviewed by:
Yusha Araf, Shahjalal University of Science and Technology, BangladeshJih-Jin Tsai, Kaohsiung Medical University Hospital, Taiwan
Copyright © 2022 Aw, Mok, Wong, Chen, Mak, Lin, Lye, Tan and Chu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Kai Sen Tan, kaistan@nus.edu.sg; Justin Jang Hann Chu, miccjh@nus.edu.sg
†These authors have contributed equally to this work