Functional pangenome analysis suggests inhibition of the protein E as a readily available therapy for COVID-2019

The spread of the novel coronavirus (SARS-CoV-2) has triggered a global emergency, that demands urgent solutions for detection and therapy to prevent escalating health, social and economic impacts. The spike protein (S) of this virus enables binding to the human receptor ACE2, and hence presents a prime target for vaccines preventing viral entry into host cells1. The S proteins from SARS-CoV-1 and SARS-CoV-2 are similar2, but structural differences in the receptor binding domain (RBD) preclude the use of SARS-CoV-1–specific neutralizing antibodies to inhibit SARS-CoV-23. Here we used comparative pangenomic analysis of all sequenced Betacoronaviruses to reveal that, among all core gene clusters present in these viruses, the envelope protein E shows a variant shared by SARS and SARS-Cov2 with two completely-conserved key functional features, an ion-channel and a PDZ-binding Motif (PBM). These features trigger a cytokine storm that activates the inflammasome, leading to increased edema in lungs causing the acute respiratory distress syndrome (ARDS)4-6, the leading cause of death in SARS-CoV-1 and SARS-CoV-2 infection7,8. However, three drugs approved for human use may inhibit SARS-CoV-1 and SARS-CoV-2 Protein E, either acting upon the ion channel (Amantadine and Hexamethylene amiloride9,10) or the PBM (SB2035805), thereby potentially increasing the survival of the host, as already demonstrated for SARS-CoV-1in animal models. Hence, blocking the SARS protein E inhibits development of ARDS in vivo. Given that our results demonstrate that the protein E subcluster for the SARS clade is quasi-identical for the key functional regions of SARS-CoV-1 and SARS-CoV-2, we conclude that use of approved drugs shown to act as SARS E protein inhibitors can help prevent further casualties from COVID-2019 while vaccines and other preventive measures are being developed.


Introduction
The current epidemic of the novel coronavirus (SARS-CoV-2) has triggered a rising global health emergency (WHO https://www.who.int/emergencies/diseases/novel-coronavirus-2019) with disruptive health, social and economic repercussions. Three major actions, confinement, detection and therapy are required to contain this pandemic. Detection and therapy can be guided by virus sequence data, which was first made available in 17 December, 2019 at gsaid.org, soon after detection of the first cases. The SARS-CoV-2 genome sequence data are a pivotal resource on its own, but can gain greater value when embedded within those of other Betacoronavirus, allowing a comparative pangenomic analysis. This approach can help identify the core genome of Betacoronaviruses and extract accessory genomic features shared by a subset of these viruses or unique to SARS-CoV-2. Whereas core genomic features are required for the virus to be functional, accessory features are candidates to provide insights into the drivers of the unique capacities of SARS-CoV-2 explaining its spread and virulence. Genome annotation and structural modelling can then be used to assess the possible functions of these accessory features and guide approaches to detection and treatment.
Here we apply a comparative pangenomic approach of all Betacoronavirus genomes sequenced thus far (Table 1), to detect the core and accessory gene cluster of this genus, and then annotate the functions, further assessed through structural analysis. Furthermore, we provide insights into importance, key features and differences in the Envelope protein, E, conserved separately in SARS, MERS and coronaviruses from other animals.

Results and Discussion
We used a pangenome approach to explore the genome of SARS-CoV-2 (MN985325 and related isolates) in comparison to other Betacornaviruses, including SARS and MERS. Previous approaches 11 proceeded by extracting individual genes, aligning them, and then establishing a phylogeny. Our approach differs in that we first cluster all sequences to then calculate the alignment. This approach, based on using panX 12 , allows us to achieve a higher sensibility in the detection of gene clusters. The resulting pangenome including all core and accessory gene clusters, interactive trees and alignments from 24 Betacoronavirus genomes, including 4 isolates from SARS-CoV-2, is available at https://pangenomedb.cbrc.kaust.edu.sa/.
There are five essential genes in coronaviruses, the Spike protein (S), Membrane glycoprotein (M), Nucleocapsid (N), envelope protein (E) and the Orf1ab (a large polyprotein known as replicase/protease), all required to produce a structurally complete viral particle 13 . These five genes define the core pan-genome. However, panX 12 only retrieved four of these (S, M, N and the ORF1ab) as components of the core pan-genome ( Figure 1). However, the fifth gene in the Betacoronavirus core genome, the envelope protein (E), varied among genomes, defining three distinct subclusters within the envelope protein E of Betacoronaviruses ( Figure 2). One of these E clusters comprises only SARS-CoV-1 (AY274119), SARS-CoV-2 and two Bat Coronaviruses (MG772933 and KF63652). Among the other two protein E gene clusters, one includes MERS (JX869059 and KC164505) and several coronaviruses from different animals, while the third cluster includes coronaviruses from two Rousettus bats. To validate the conservation and specificity of the SARS E cluster, we compared the SARS-CoV-2 E sequence with all known sequences in NCBI's non redundant (NR) protein database. This comparison showed, consistent with our pan-genome analysis, that the subcluster of the protein E gene appearing in SARS-CoV-2, SARS and two bat coronaviruses, is defined by the same essential functional features, which are exactly conserved (100% similarity) between bats and SARS-CoV-2, but differ slightly (95% similarity with one deletion and three substitutions) from SARS ( Figure 2B and supplementary figure 7). A recent study compared genomes and gene families of all alpha, beta, delta and gamma coronaviruses, but did not explore the variability present within the protein E gene sequences, which we did through the clustering approach reported above.
We found five accessory gene clusters (15% of the accessory genome of genus Betacoronavirus) common to SARS-like coronaviruses, including SARS-CoV-2, out of 33 total accessory gene clusters. Four of these gene clusters (ORF3a, ORF6, ORF7a, and ORF8, Fig. 1A), appear only in SARS-CoV-2, bat-SL-CoVZC45 and the SARS virus. A gene cluster annotated as ORF10 appeared to be unique to SARS-CoV-2 ( Fig. 1). However, a blast-based search against the DNA of all sequences in NCBI's NT database showed that the ORF10 sequence matches with 100% coverage in DNA of other SARS-like genomes, but there is no open-reading frame predicted for those genomes in the matching region ( Fig 1A). Hence, the apparent uniqueness of ORF10 might simply be an artifact of the GeneBank annotation pipeline that did not predict an ORF for this gene in other matching genomes. However, for completeness, we report a functional analysis of ORF10 as Supplementary material (Supplementary material 1) since at present there is no experimental evidence available against or in favor of this being a real gene.
Compared to the essential genes S, M, N, E and ORF1ab, other genes such as ORF3a and ORF7a are poorly characterized, according to GenBank annotations (https://www.ncbi.nlm.nih.gov/genome/proteins/86693?genome_assembly_id=760344). Our bioinformatic protein structural analysis confirmed that SARS-CoV-2 ORF3a (cluster 12), also known as viroporin, and ORF7a (cluster 17), that is same as ORF8a 14 in SARS, retain the structural features observed in other SARS viruses, namely a multi-pass transmembrane domain and a cytoplasmic b-barrel or b-sandwich fold (ORF3a) and an N-terminal secpathway signal peptide, cleaved after residue 15, an immunoglobulin-like b-sandwich fold stabilized by two cysteine di-sulfide bonds, and a C-terminal single-pass transmembrane helix (ORF7a; Supplementary Figures 5,6). Hence these SARS-CoV-2 proteins are expected to act in the same way as they do in other SARS, namely as accessory proteins mostly localized in the endoplasmic reticulum-Golgi intermediate compartment, but also occurring on the cell membrane where they enhance viral pathogenicity and mortality through protein-protein interactions [15][16][17] .
The remaining two accessory gene clusters, ORF6, ORF8, present in SARS-CoV-2, SARS and bat-SL-CoVZC45 are functionally uncharacterized. All are very short polypeptides (61, 121 and 38 residues for ORF6 and 8, respectively) with a large percentage of hydrophobic residues (62%, 56%, respectively). None of the proteins have trans-membrane regions, but ORF8 has an N-terminal sec-pathway signal peptide with a cleavage site after residues 15, suggesting that it is secreted into the extracellular space ( Supplementary Figures 5,6). Following signal peptide cleavage, the ORF8 protein core is predicted to consist largely of b-strands and features 7 cysteines. We predict that this protein adopts a cysteine disulfide-bond stabilized b-sandwich structure similar to the soluble domain of ORF7a (Supplementary Figures 5,6), inferring that ORF8 also functions as ligand binding module. ORF6 consists of a long amphipathic helical region, followed by an acidic tail. Our deep-learning-based annotation (DeepGOPlus 18 ) suggested that ORFs 6, 7, 8 are involved in the regulation of molecular functions, either in response to stimuli (ORFs 8), or in cellular component biogenesis and organization (ORF6) (Suppl. Figure 2-4). Collectively, our analysis suggested no major function or gene cluster to be unique to SARS-CoV-2. Moreover, our analysis is in agreement with ORFs 3a, 6, 7a, 8 and The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/2020.02.17.952895 doi: bioRxiv preprint being accessory non-essential proteins, that would be inefficient targets for COVID-2019 therapy.
The S protein binds to the host receptor ACE2, and hence presents a prime target for preventing viral entry into host cells 1 . The S proteins from SARS-CoV-1 and SARS-CoV-2 are similar 2 , but structural differences in the receptor binding domain (RBD) preclude the use of SARS-CoV-1-specific neutralizing antibodies to inhibit SARS-CoV-2 3 . We therefore focused our attention on the protein E, which is highly similar in SARS-CoV. The E protein, also a viroporin 19 , was previously confirmed as a determinant of pathogenicity 4,5 in SARS-CoV-1. Protein E was used as target for SARS antivirals 6 , and studies using SARS-CoV with lacking or mutated protein E as vaccine candidates showed promising results [20][21][22][23] . Furthermore, the E gene is one of the key genes used in identification of SARS-CoV-2 through RT-PCR 7,24 .
The SARS-Cov-1 protein E has ion channel (IC) activity and features a postsynaptic density-95/discs large/zona occludens-1 (PDZ)-binding motif (PBM). Both IC and PBM were required for SARS-CoV to induce virulence in mice 14 . The protein E from SARS-CoV-2 differs from that of SARS-CoV-1 only by three substitutions and one deletion (Figure 2A Studies on SARS-CoV-1 demonstrated that the E protein uses the IC and PBM to trigger a cytokine storm that activates the inflammasome, leading to increased edema in lungs. Ultimately these events result in ARDS 4-6 , a leading cause of death in SARS-CoV-1 and SARS-CoV-2 infection 7,8 . However, the drugs Amantadine, Hexamethylene amiloride and also BIT225 (BIT225 is in clinical trials) completely block the IC activity of SARS-CoV-1 9,10 and restrict its reproduction, leading to better survival of the animal host. The PBM enables the E protein to interact with the cellular PDZ domains. Its association with the tandem PDZ domains of syntenin activates the MAP kinase p38, triggering the overexpression of inflammatory cytokines. Consequently, inhibition of p38 by SB203580 increases the survival of the host 5 .
The high level of conservation of the IC and PBM within clade 2, including SARS-CoV-1 and SARS-CoV-2, strongly suggests that the available inhibitors for the E protein of SARS-CoV-1 should also be active on SARS-CoV-2.

Conclusion
Therapies to combat the spread of SARS-CoV-2 and the lethality caused by the resulting COVID-2019 are currently focusing primarily on S, the viral spike protein 3,25,26 . However, despite high similarities between proteins S between SARS-CoV-1 and SARS-CoV-2 viruses, existing neutralizing antibodies are ineffective against SARS-CoV-2 3 . Hence, new antibodies that bind specifically to the SARS-CoV-2 spike protein need to be developed, tested and approved for human use, which would be a time-consuming process.
Our pangenomic analysis suggests that the protein E of all SARS viruses preserves its critical motifs used for pathogenesis. Hence, we predict that the readily available and approved inhibitors (Amantadine, Hexamethylene amiloride and SB203580), that proved efficient for alleviating ARDS in SARS-CoV-1-infected animal models should also be effective against SARS-CoV-2. While not preventing future spreading of the virus, these inhibitors might reduce the mortality rate while effective vaccines are being developed.

Methods
Collection of Betacoronavirus strains. The NCBI genome assemblies page was searched, on January 26, for taxonid of genus Betacoronavirus, resulting in a list of 22 genomes, including four isolates of SARS-CoV-2. We downloaded in addition to this initial list two more assemblies one for a bat coronavirus (MG772933) and another for MERS (KF600630), see Table 1. For easy identification of genes and annotations from GenBank, we included unique locus_tag identifiers containing locus id and an index number. Metadata was collected to define user friendly strain names and to understand phylogenetic trees in the context of virus host, country, taxon id and data release date.
Creating a pangenome for Betacoronavirus GenBank annotation and sequences were used as a starting point for this pangenomic analysis study. These GenBank annotations contain heterogeneous identifiers for genes so we included unique locus tags containing the strain accession number and gene index number for computational processing. In addition, we collected metadata for these coronaviruses for displaying on the interactive phylogenetic trees for analysis later on. Gene sequences were compared to each other for sequence similarity, followed by gene clustering. Phylogenetic analysis was carried out on the core gene clusters to produce a species tree and similarly for each gene cluster to produce gene trees. For this pangenome analysis we used panX 12 (Ding et al., 2018) that takes GenBank files as input and extracts gene coordinates, annotations from GenBank and sequences from these genomes on both nucleotide and protein levels. For interactive visualization of pangenome clusters, alignments and phylogenetic trees we used panX visualization module obtain from github, https://github.com/neherlab/pan-genomevisualization/. Other databases showing greater detail about many viruses, such as Virus Pathogen Resource (VIPR, https://www.viprbrc.org) and GSAID (https://www.gisaid.org) are very useful but these resources do not provide gene clustering or pan-genome analysis.
Re-annotation of proteins from Betacoronavirus genomes using KAUST Metagenomic Analysis Platform (KMAP) annotation pipeline. To annotate uncharacterized genes in Betacoronavirus, we utilized all available protein sequences available from the clustering analysis as an input to our Automatic Annotation of Microbial Genomes (AAMG) pipeline 27 , available via KAUST Metagenomic Analysis Platform (http://www.cbrc.kaust.edu.sa/kmap). AAMG uses sequence-based BLAST to UniProtKB (www.uniprot.org), KEGG (www.kegg.jp) and also Protein Family (PFam) domain detection using InterPro (www.ebi.ac.uk/interpro).
Prediction of function using DeepGOPlus. We used all genes from Betacoronavirus dataset in order to explore Gene Ontology predictions through our in-house DeepGOPlus 18 tool that combines deep convolutional neural network (CNN) model with sequence similarity-based predictions to derive relevant functional classes from Gene Ontology alongside a confidence score. Results were included with a score 0.1 and filtered for class specificity. Predicted Gene Ontology (GO) terms with a confidence score from DeepGOPlus were visualized using QuickGO (https://www.ebi.ac.uk/QuickGO/) for SARS-CoV-2 specific gene clusters. DeepGOPlus is available on github, https://github.com/bioontology-research-group/DeepGOPlus Structure based prediction of function ProtParam (https://web.expasy.org/cgi-bin/protparam/protparam) was used for calculating protein features. Phobius (http://phobius.sbc.su.se/) and SignalP-5.0 were used for prediction of transmembrane regions and signal peptides 28 . Jpred4 was used to calculate secondary structure features (http://www.compbio.dundee.ac.uk/jpred4/). 3D modelling was carried out using SwissModel (homology modelling) 29 or QUARK (ab initio structure predictions) 30 .

Author contributions
IA, TG and CMD conceived and designed the research, IA led the study and conducted the data analysis, A.K. helped developed the web-based resource and led computational components of the study, MK and SA contributed to the functional structural analysis, and AP contributed with inferences on virulence. IA and CMD developed the first draft of the manuscript and all authors contributed to writing and improving the manuscript, and approved the submission.

Supplementary Information 1, ORF10 gene in SARS-CoV-2:
A short gene predicted in in SARS-CoV-2, known as ORF10, appears unique only to SARS-CoV-2. A sequence comparison of this gene to DNA of all Betacoronavirus genomes in NCBI shows matches in SARS and SARS-like genomes with 100% coverage, although no gene has been predicted in matching location of ORF10 in SARS or SARS-like genomes. This raises a question whether ORF10 is an artifact of annotation in SARS-CoV-2 where it has been predicted or an artifact in SARS and SARS-like genomes where it has not been predicted. There is no experimental data available to approve or disprove the reality of this gene. For completeness purposes we included ORF10 in our functional and structural analysis of SARS-CoV-2 but not included in the main text.
Supplementary Figure 4-6 shows predicted Gene Ontology and structural information about ORF10 using deep learning Gene Ontology Prediction tool (DeepGOplus). It shows that, similar to ORF8, this gene regulates molecular functions and contributes to response to stimulus. It is predicted to be localized only in the extracellular region of the host cells. OF10 is predicted to harbor a long helix followed by a b-strand and appears to be localized only in the extracellular region of the host cells. It has been concluded that SARS-CoV-2 accessory genes, including ORF10, carry a helper function and do not serve as a major target for detection or therapy of COVID-2019.
Supplementary Figure 1. DeepGOPlus predicted Gene Ontology Graph for SARS-CoV-2 envelope protein, E. According to predictions, the protein can be localized in the host cell Golgi apparatus, vacuoles, endoplasmic reticulum, and membrane cell periphery of the host cell. It contributes to the process of disruption of the host cell and viral budding from Golgi membrane. Figure 7. Alignment of SARS-CoV-2 E protein comparison to hits available in NCBI's NR database. The first two appear from SARS-CoV-2 and rest of the hits appear from SARS and SARS-like Betacoronaviruses. This result shows SARS-CoV-2 E protein is conserved only in SARS clade and E from MERS or other animal coronaviruses clusters in different clusters (see Figure 2).