Extreme Genomic CpG Deficiency in SARS-CoV-2 and Evasion of Host Antiviral Defense

1Coronaviruses (CoVs) evolve in mammalian hosts and carry genomic signatures of their host-specific environment, especially the host-specific antiviral and RNA modification activities. Many pathogenic single-stranded RNA viruses, including CoVs, exhibit strong CpG deficiency. Two mammalian enzymes are inferred to contribute to the observed CpG deficiency. The zinc finger antiviral protein (ZAP, known as ZC3HAV1 in mammals or hZAP in human), a key component in mammalian interferon-mediated immune response, binds specifically to CpG dinucleotides in viral RNA genomes via its RNA-binding domain. ZAP inhibits viral replication and mediates viral genome degradation. ZAP has two isoforms (ZAP-L and ZAP-S); both participate in initiating antiviral activities but only ZAP-S mediates the return to homeostasis after the antiviral response. ZAP acts against not only retroviruses, such as HIV-1, but also Echovirus 7 and Zika virus, both being positive-sense single-stranded RNA viruses like CoVs. In particular, selection against CpG in viral RNA disappears in ZAP-deficient cells, suggesting that ZAP may be the only cellular agent targeting CpG in viral RNA genomes.

2Experimental evidence is consistent with the interpretation that CpG deficiency in RNA viruses has evolved in response to these cytoplasmic CpG-specific antiviral activities. During natural evolution of HIV-1 within individual patients, viral fitness decreased with increasing CpG dinucleotides. Experimental increase of CpG dinucleotides in CpG-deficient viral genomes consistently leads to strong decrease in viral replication and virulence, prompting the proposal of vaccine-development strategies involving increasing CpG to attenuate pathogenic RNA viruses.

3Another antiviral enzyme is APOBEC3G, found in innate immune cells. APOBEC3G was originally thought specific to single-stranded DNA, such as reverse-transcribed HIV-1, but is now known to modify a variety of RNA viruses, deaminating C to U. This would be effective against RNA viruses if the deaminated sites are functionally important. APOBEC3G copurifies with highly edited mRNA substrates and therefore could act on CoV genomes which are positive-sense single-stranded RNA. Although APOBEC3G is not strongly CpG-specific, it could contribute to CpG deficiency when coupled with ZAP-mediated antiviral activities targeting CpG. Modification of CpG to UpG in nonfunctional regions could reduce viral susceptibility to CpG-mediated attack by ZAP relative to viruses with unmodified CpG dinucleotides.

4Both ZAP and APOBEC3G exhibit tissue-specific expression patterns in human. Both are expressed in lungs, but ZAP is the most highly expressed where lymphocytes are the most abundant (bone marrow, lymph node, appendix, and spleen), whereas APOBEC3G is the most highly expressed in lymph node, spleen, and testis. A severely CpG deficient virus may indicate an evolutionary history in ZAP-abundant tissues, such as strongly CpG-deficient HIV-1, infecting host T cells in lymph organs where ZAP is abundant. The presence of such viruses indicates that they have found ways to evade ZAP-mediated cellular antiviral defense.

5The differential expression of ZAP and APOBEC3G in different host or host tissues is expected to leave different genomic signatures on viral RNA genomes. We may use the conventional index of CpG deficiency implemented in DAMBE:

6The index is expected to be 1, with no deficiency or excess; <1, if deficient; and >1, if excess. The 1252 betacoronavirus (BetaCoV) full-length genomes deposited in GenBank (of which 1,127 are unique) have mean±SE value of 0.516 ± 0.0017 for ICpG, which is significantly (P < 0.0001) smaller than their null expectation of 1.

7If a CoV infects a different host tissue with different ZAP abundance, then its RNA genome will experience different selection pressure against its CpG. This difference in cellular antiviral activity would result in differences in ICpG during viral genomic evolution. In contrast, a CoV infecting a specific host tissue for a long time would experience the same cellular antiviral and RNA modification environment and is consequently expected to have similar and stable ICpG.

8SARS-CoV-2 and its most closely related known relative (BatCoV RaTG13) have the lowest ICpG among its close relatives, both being outliers in a plot of viral genomic ICpG vs. GC% (fig. 1). Three groups of BetaCoV most closely related to SARS-CoV-2 are represented in figure 1. Group 1 consists of genomes of human CoV-HKU1 which is found only in human and circulates among human populations without any dependence on other mammalian species as intermediate or reservoir species. Group 2 includes BetaCoV 1 genomes with two types of hosts: 1) ungulates (with bovine and equine CoV as well as porcine hemagglutinating encephalomyelitis virus); and 2) human, with CoV-OC43 being a recent derivative of bovine CoVs. Group 3 are all SARS-related CoVs from three types of hosts: 1) Rhinolophus bats which serve as a natural reservoir of SARS-related CoVs and the new SARS-CoV-2, 2) civets (from which CoV genomes with 99.6% identity to SARS virus genomes were identified; ), and 3) human patients infected by SARS-CoV-2. Figure 1 shows that genomic GC% and ICpG can differ among different viral lineages in the same host or among different hosts for the same viral lineage.

9The most striking pattern in figure 1 is an isolated but dramatic shift in the lineage leading to BatCoV RaTG13 which was reported to be sampled from a bat (Rhinolophus affinis) in Yunnan Province in 2013 but only sequenced by Wuhan Institute of Virology after the outbreak of SARS-CoV-2 infection in late 2019. This bat CoV genome is the closest phylogenetic relative of SARS-CoV-2, sharing 96% sequence similarity. Many studies have shown an association between decreased CpG (low ICpG) in viral RNA genomes and increased virulence, not only in HIV evolving within individual patients but also in experimentally CpG dinucleotide-enriched viral genomes. The association between decreased CpG and increased virulence in RNA viruses is mainly due to interferon-induced ZAP protein which binds to CpG dinucleotides in viral RNA genomes by its RNA-binding domain, inhibits viral replication, and facilitates viral genome degradation. Thus, a decreased ICpG in a viral pathogen suggests an increased threat to public health, but an increased ICpG decreases the threat because such viral pathogens, with increased ICpG and reduced virulence, would be akin to natural vaccines. Many viral researchers have in fact proposed vaccine development by increasing CpG in viral RNA genomes.

10In this context, it is unfortunate that BatCoV RaTG13 was not sequenced in 2013; otherwise, the downshifting in ICpG might have served as a warning due to two highly significant implications. First, the virus likely evolved in a tissue with high ZAP expression which favors viral genomes with a low ICpG. Second and more importantly, survival of the virus indicates that it has successfully evaded ZAP-mediated antiviral defense. In other words, the virus has become stealthy and dangerous to humans.

11The ICpG value for BatCoV RaTG13 genome is 0.40875, much lower than ICpG values observed in all other BetaCoV genomes sampled from bat species in the genus Rhinolophus. There are 56 BetaCoV genomes sampled from Rhinolophus bats inhabiting south and southeastern Asia (but mostly from central and southern China). Nature had essentially inoculated BetaCoVs into various Rhinolophus lineages and allowed genomic evolution to happen (supplementary fig. 1, Supplementary Material online). Although genomic ICpG values have fluctuated in different viral lineages, only BatCoV RaTG13 has been observed to possess an extraordinarily low ICpG (supplementary fig. 1, Supplementary Material online). This suggests that the ancestor of BatCoV RaTG13 and SARS-CoV-2 may have evolved in a mammalian tissue with high expression of ZAP and emerged with an unusually low ICpG. This mammalian tissue likely is not in Rhinolophus bats because low ICpG has not been observed in other BatCoV lineages (supplementary fig. 1, Supplementary Material online). Identifying a virus with comparably low ICpG would suggest candidate host species possessing tissues with cellular environments that select strongly against CpG in viral genomes.

12Among all BetaCoVs available in GenBank on February 3, 2020, there are 1,127 unique genomes of which 927 genomes have explicit host designations (supplementary file Betacoronavirus_CpG.xlsx, Supplementary Material online). Surprisingly, no available BetaCoV genome from diverse natural hosts has a genomic ICpG and GC% combination close to that observed in SARS-CoV-2 and BatCoV RaTG13 (fig. 2). BetaCoV lineages parasitizing Rhinolophus bats overall have relatively low ICpG values (fig. 2).

13BetaCoV infecting dromedary camels offers a weak hint that camel digestive system may select more strongly against CpG in viral genomes than camel respiratory system. Camel CoVs form two clusters. One cluster overlaps completely with MERS viruses (fig. 2) that infect mammalian respiratory system. The other cluster includes camel CoV HKU23 strains positioned close to bovine CoV (grouped under “Ungulate_CoV” in figs. 1 and 2), both belonging to Embecovirus and infecting mainly mammalian digestive system but also respiratory systems. Those viruses infecting camel digestive system have lower genomic ICpG and GC% than those infecting camel respiratory system (fig. 2).

14To search for a mammalian host with the potential to select viral lineages with low ICpG values, I expanded the search to include all complete alphacoronavirus (AlphaCoV) genomes (supplementary file alphacoronavirus_CpG.xlsx, Supplementary Material online). All complete AlphaCoV genomes (>27,000 nt) with explicit host information are plotted in ICpG and GC content in figure 3. Five points are worth highlighting. First, only genomes from canine coronaviruses (CCoVs), which had caused a highly contagious intestinal disease worldwide in dogs, have genomic ICpG and GC% values similar to those observed in SARS-CoV-2 and BatCoV RaTG13 (fig. 3A). The genome (accession no. KP981644) is from the most virulent pantropic CCoV invading multiple canine organs. It belongs to a clade with the lowest observed ICpG values (fig. 3B).

15Second, canids, like camels, also have CoVs infecting their respiratory system (canine respiratory CoV or CRCoV belonging to BetaCoV). There are two genomes sequenced for CRCoV (accession nos. JX860640 and KX432213). Their genomic ICpG values are 0.4756 and 0.4684, respectively, substantially higher than those for CCoVs infecting the digestive system (fig. 3A). Thus, similar to the pattern observed in CoVs infecting camels, CCoVs infecting canine digestive system have ICpG much lower than CRCoVs infecting canine respiratory system.

16Third, none of the available AlphaCoV genomes from bats or other mammalian host species possess genomic ICpG and GC% values similar to those observed in SARS-CoV-2 and BatCoV RaTG13 (fig. 3). Thus, although AlphaCoV infects a diverse array of bat lineages, these bat tissues do not seem to generate AlphaCoV strains with low ICpG values comparable with SARS-CoV-2 and BatCoV RaTG13.

17Fourth, I want to highlight one data point involving a CCoV genome represented as a green dot in figure 3 (highlighted by a green arrow in fig. 3A, genome accession KC175339). The CCoV has a genomic GC% of 38.17% and ICpG of 0.4986, much higher than the rest. The virus was originally isolated from a dog but had been propagated extensively in cell culture before being sequenced (Whittaker GR, personal communication). Viruses are propagated in cells that expresses the right cellular receptor for viral entry, but do not mount an immune response to kill the virus or get killed by the virus. The consequent relaxation of selection against the virus (and against CpG in the CCoV genome) in cell culture would allow CpG in the viral RNA genome to rebound through mutation, which would explain the increased ICpG (KC175339 in the phylogeny in fig. 3B). This process of regaining CpG is reminiscent of CpG-specific methylation in Mycoplasma species where CpG was regained when some lineages lost CpG-specific methyltransferases, with a fast-evolving lineage (M. pneumoniae) regaining CpG faster than a slow-evolving lineage (M. genitalium; ). This rapid change in ICpG with environmental change as shown in figure 3B has two important implications. First, it suggests the feasibility of tracking certain host-switching or tissue-switching events (which would be impossible if it takes hundreds of years for a virus to change ICpG). Second, many experimental studies have demonstrated attenuated virulence in RNA viruses with increasing CpG in the viral RNA genome and suggested this as an efficient means of vaccine development. The observed increase in ICpG in the CCoV genome through cell culture propagation shows a simple way of increasing CpG by simply propagating the virus without selection against CpG dinucleotides in the viral genome.

18Fifth, the cellular receptor for SARS-CoV-2 entry into the cell is angiotensin I converting enzyme 2 (ACE2; ). ACE2 is pervasively expressed in human digestive system, at the highest levels in small intestine and duodenum (fig. 3C), with relatively low expression in lung. This suggests that mammalian digestive system is likely to be infected by CoVs. This is consistent with the interpretation that the low ICpG in SARS-CoV-2 was acquired by the ancestor of SARS-CoV-2 evolving in mammalian digestive system. The interpretation is further corroborated by a recent report that a high proportion of COVID-19 patients also suffer from digestive discomfort. In fact, 48.5% presented with digestive symptoms as their chief complaint.

19Humans are the only other host species observed to produce CoV genomes with low-genomic ICpG values, as shown by the cluster of human AlphaCoV NL63 genomes (fig. 3). This virus mainly not only infects the respiratory system but also causes digestive problems in 33% of the patients reporting respiratory problems. In a comprehensive study of the first 12 COVID-19 patients in the United States (Midgley and ), one patient reported diarrhea as the initial symptom before developing fever and cough (Midgley and ). Stool samples from 7 out of 10 patients tested positive for SARS-CoV-2, including 3 patients with diarrhea (Midgley and ), corroborating a previous report of SARS-CoV-2 detection in stool. In particular, live SARS-CoV-2 virus was isolated from stool of a COVID-19 patient. In this context, it is significant that BatCoV RaTG13, as documented in its genomic sequence in GenBank (MN996532), was isolated from a fecal swab. These observations are consistent with the hypothesis that SARS-CoV-2 has evolved in mammalian intestine or tissues associated with intestine.

20Figures 1 and 2 do not include all BetaCoVs from all hosts, so BetaCoVs from other mammalian species may also possess low ICpG values. One example is viruses isolated from pangolins. Nine SARS-CoV-2-like genomes have recently been isolated and sequenced from pangolin and deposited in GISAID database (gisaid.org). The one with the highest sequence coverage (GISAID ID: EPI_ISL_410721) has an ICpG value of 0.3929, close to the extreme low end of ICpG values observed among available SARS-CoV-2 genomes. Thus, SARS-CoV-2, BatCoV RaTG13, and those from pangolin may either have a common ancestor with a low ICpG or have convergently evolved low ICpG values.

21Other than ZAP and ABOBEC3G, the enigmatic DNA methyltransferase2 (Dnmt2;,; ), originally thought to be a Dnmt, may also contribute to viral RNA modification. However, Dnmt2 appears to methylate only small RNA. For this reason, it may not be important in shaping ICpG in large CoV RNA genomes, although it has been observed to relocate from the nucleus to cytoplasmic stress granules, where it may participate in the methylation of mRNA.

22These observations allow formation of a hypothesis for the origin and initial transmission of SARS-CoV-2. First, the ancestor of SARS-CoV-2 and BatCoV RaTG13 infected the intestine of a mammalian species (e.g., canids or human ingesting bat meat). Second, the presumably strong selection against CpG in the viral RNA genome in canid intestine resulted in rapid evolution of the virus, with many CpG → UpG mutations leading to reduced genomic ICpG and GC%. The licking of anal regions in canids during mating and other circumstances facilitate viral transmission from the digestive system to the respiratory system. Finally, the reduced viral genomic ICpG allowed the virus to evade human ZAP-mediated immune response and became a severe human pathogen. Because SARS-CoV-2, TG13 and the related pangolin-derived coronaviruses all have a low-GC genome, the simplest hypothesis is that the low-GC genome was gained in their common ancestor. However, it is also possible that the viral lineage gained the low-GC genome only recently in the digestive system of a canid and spread to other species. This suggests the importance of monitoring SARS-like CoVs in feral dogs in the fight against SARS-CoV-2.

23Although the specific origins of SARS-CoV-2 are of vital interest in the current world health environment, this study more broadly suggests that important evidence of viral evolution can be revealed by consideration of the interaction of host defense with viral genomes, including selective pressure exerted by host tissues on viral genome composition.

24Supplementary data are available at Molecular Biology and Evolution online.