Dynamically evolving novel overlapping gene as a factor in the SARS-CoV-2 pandemic

Abstract

Understanding the emergence of novel viruses requires an accurate and comprehensive annotation of their genomes. Overlapping genes (OLGs) are common in viruses and have been associated with pandemics but are still widely overlooked. We identify and characterize ORF3d, a novel OLG in SARS-CoV-2 that is also present in Guangxi pangolin-CoVs but not other closely related pangolin-CoVs or bat-CoVs. We then document evidence of ORF3d translation, characterize its protein sequence, and conduct an evolutionary analysis at three levels: between taxa (21 members of Severe acute respiratory syndrome-related coronavirus), between human hosts (3978 SARS-CoV-2 consensus sequences), and within human hosts (401 deeply sequenced SARS-CoV-2 samples). ORF3d has been independently identified and shown to elicit a strong antibody response in COVID-19 patients. However, it has been misclassified as the unrelated gene ORF3b, leading to confusion. Our results liken ORF3d to other accessory genes in emerging viruses and highlight the importance of OLGs.

Keywords: ORF3d; SARS-CoV-2; evolutionary biology; genome annotation; infectious disease; microbiology; natural selection; overlapping genes; pandemic; virus.

Figure 1.. Gene repertoire and evolutionary relationships of Severe acute respiratory syndrome-related coronavirus species members.

Only genes downstream of ORF1ab are shown, beginning with the Spike gene S. (A) Four types of genes and their relative positions in the SARS-CoV-2 Wuhan-Hu-1 genome (NCBI: NC_045512.2). Genes are colored by type: novel overlapping genes (OLGs) (gold; ORF3d only); conserved OLGs (burgundy); accessory (green); and structural (blue). Note that ORF3b has been truncated relative to SARS-CoV genomes, whereas ORF8 remains intact (i.e. has not been split into ORF8a and ORF8b). (B) Genes with intact ORFs in each of 21 Severe acute respiratory syndrome-related coronavirus genomes. Gene positions are shown relative to each genome, i.e. homologous genes are not precisely aligned. Only the full-length isoforms of ORF3a and ORF3d are shown (for shorter isoforms, see Table 1). Note that the first 20 codons of ORF3d overlap the last codons of ORF3c (Supplementary file 1), such that the beginning of ORF3d involves a triple overlap (ORF3a/ORF3c/ORF3d). ORF3b is full-length in only three sequences (SARS-CoV TW11, SARS-CoV Tor2, and bat-CoV Rs7327), while the remaining sequences have premature STOP codons (Supplementary file 1). ORF8 is not novel in SARS-CoV-2 (contra Chan et al., 2020), but is intact in all but five sequences (split into ORF8a and ORF8b in SARS-CoVs TW11 and Tor2; deleted in bat-CoVs BtKY72, BM48-31, and JTMC15). ORF9b and ORF9c are found throughout this virus species, yet rarely annotated in genomes at NCBI.

Figure 1—figure supplement 1.. Codon permutation analysis to identify candidate overlapping genes in all three forward-strand reading frames of the SARS-CoV-2 genome.

Genome positions of known and hypothesized genes are indicated with horizontal red lines. Reading frames 1, 2, and 3 begin at positions 1, 2, or 3 of the Wuhan-Hu-1 reference genome, respectively, with genome coordinates shown at the bottom (Supplementary file 1). Yellow indicates low p-values (natural logarithm scale), while gray indicates absence of a STOP-to-STOP ORF longer than 30 codons or presence outside of an annotated non-overlapping gene (not tested). Note that low p-values are expected for the non-overlapping genes. Only the first of four SARS-CoV-2 ORFs homologous to ORF3b is shown (Figure 4; Supplementary file 1).

Figure 2.. SARS-CoV-2 gene expression in ribosome profiling and mass spectrometry datasets.

Ribosome profiling (Ribo-seq) data were obtained from Finkel et al., 2020; mass spectrometry data were obtained from Davidson et al., 2020 and Bezstarosti et al., 2020. Reads were trimmed to their first (5′) nucleotide to minimize statistical dependence while preserving reading frame. Results are shown after pooling samples treated with harringtonine (SRR11713360, SRR11713361) or lactimidomycin (SRR11713358, SRR11713359) at 5 hr post-infection. (A) Ribosome profiling coverage (read depth) near translation initiation sites, measured as mean reads per million mapped reads. Only reads of length 29–31 nucleotides were used, chosen for their enrichment at the start sites of highly expressed annotated genes (Figure 2—figure supplement 4). Light blue backgrounds denote fully (internal) overlapping genes. Annotated genes show an accumulation of 5′ ends of protected reads upstream of the gene’s start site (vertical gray dotted lines), peaking near the ribosome P-site offset of −12 nt (red = maximum depth). Distributions are largely consistent across individual samples (Figure 2—figure supplement 1). The ranges of the y-axes vary according to expression level, with the most highly expressed gene (N; Figure 2—figure supplement 2) having the largest range. (B) Correlation between protein expression as estimated by mass spectrometry and ribosomal profiling. ‘iBAQ prop’ refers to the relative (proportion of maximum) protein intensity-based absolute quantification (iBAQ) value (Schwanhäusser et al., 2011). Genes are denoted by shape and ordinally colored by iBAQ prop from high (red) to low (blue). ‘Upstream peak’ refers to the maximum read depth observed at the approximate P-site offset (red bars in A), while mean read depths were measured across each gene using 30 nt reads separately for in-frame (codon position 1) and out-of-frame (codon positions 2 and 3) sites (non-overlapping gene regions only, except for ORF9b). (C) Reading frame of ribosome profiling reads in the ORF3c/ORF3d region of ORF3a. Solid lines show the fraction of reads in each frame, summed across samples in sliding windows of 30 nt (step size = 1 nt; read length = 30 nt). Color denotes frame: green = reference frame (ORF3a); burgundy = ss13, the forward-strand frame encoding ORF3c, whose codon position 3 overlaps codon position 1 of ORF3a; and gold = ss12, the forward-strand frame encoding ORF3d, whose codon position 2 overlaps codon position 1 of ORF3a. Values are shown for the central nucleotide of each window, with shaded regions corresponding to 95% binomial confidence intervals. Alternative frame translation is suggested where a given frame (solid line) exceeds its average across the remainder of the gene (horizontal dashed line; non-OLG regions of ORF3a). Vertical gray dotted lines indicate gene start and end sites. Gray bars show read depth for each window, with a maximum of 2889 reads at genome position 25442.

Figure 2—figure supplement 1.. Ribosome profiling coverage near translation initiation sites for individual samples.

Details as in Figure 2A, except samples were not pooled and genes are arranged vertically by start site in the genome from 5′ to 3′ (top to bottom). Treatments are indicated as Harr = harringtonine and LTM = lactimidomycin at 5 hr post-infection (hpi).

Figure 2—figure supplement 2.. Correlation between protein expression as estimated by mass spectrometry and ribosomal profiling, for individual treatments and codon positions.

Details as in Figure 2B, except estimates are shown for all sites in a gene (All frames) and out-of-frame reads are shown separately for codon positions (Codon pos) 2 and 3. ORF9b is a fully (internal) overlapping gene occurring within N, with codon position 3 of ORF9b corresponding to codon position 1 of N (i.e. frame ss13 of N). As a result, codon position 3 of ORF9b is predominantly occupied by ribosomes translating N, the most highly expressed gene, leading to overestimation of ORF9b expression (Codon pos 3; right). Ribosome profiling read depth at codon position 1 exhibits a marginal Spearman’s rank correlation with mass spectrometry for cycloheximide (CHX), harringtonine (Harr), and lactimidomycin (LTM) treatments; however, as expected, all correlations for RNA-seq (mRNA), as well as reads mapping to other codon positions, were not significant.

Figure 2—figure supplement 3.. Reading frames occupied by the 5′ ends of ribosome profiling reads as a function of read length.

Processed, aligned reads from Finkel et al., 2020 were trimmed to their 5′ ends (first nt). Fractions of reads mapping to each codon position (x-axis) are shown as a function of read length (y-axis; 26–32 nt) after pooling samples according to treatment (CHX = cycloheximide, Harr = harringtonine, LTM = lactimidomycin, mRNA = RNA seq) and time (hours post-infection; hpi). Codon position is defined relative to gene start sites of annotated non-overlapping genes, excluding any overlapping gene regions. Gray bars show relative total counts for reads of different lengths, with different scales for 5 hpi (lower coverage; max = 284,608 reads at CHX/30nt) and 24 hpi (higher coverage; max = 3,149,251 reads at Harr/29nt). At 5 hpi, reads of length 29 and 30 nt tend to have 5′ ends mapping to the correct protein-coding reading frame (codon position 1) of the gene in which they occur, and also tend to be most abundant, as expected given the size of ribosome-protected fragments (~30 nt). Samples from 24 hpi have flatter (left-shifted) read length distributions and no clear frame signal. The horizontal dashed gray line corresponds to the expectation for RNA fragments not protected by ribosomes, that is 33% of reads mapping to each codon position (random).

Figure 2—figure supplement 4.. Ribosome profiling read accumulation at gene start sites as a function of read length.

Points indicate the percentage of all reads in the region shown whose first nucleotides (5′ ends) map to each genome position, colored according to read length. In each case, the longer read lengths have 5′ ends that predominate at approximately −12 nt from the gene start site (horizontal dashed lines, right), corresponding to the ribosome P-site offset (center of yellow highlighted region). A peak in total read depth (gray bars) is also observed at this position, as seen in Figure 2A. Highly expressed non-overlapping (N, M, S, and ORF7b) and overlapping (ORF9b) genes are shown on the left, while ORF3a and its hypothesized overlapping genes are shown on the right. Note that ORF9b begins shortly after N, such that the upstream peak of N occurs at position −22 nt relative to the start site of ORF9b.

Figure 2—figure supplement 5.. Reading frame of ribosome profiling reads in the ORF3c/ORF3d region of ORF3a for individual treatments.

Details as in Figure 2C, except results are shown separately after pooling the two biological replicates for each treatment (CHX = cycloheximide, Harr = harringtonine, LTM = lactimidomycin, mRNA = RNA seq) (Finkel et al., 2020). Genes and their reading frames are denoted by color according to Figure 1, with Frames 1, 2, and 3 corresponding to the frames beginning at codon positions 1, 2, and 3 of ORF3a, respectively. Clear evidence for translation in the appropriate reading frame is observed for ORF3c (burgundy) and ORF3d/ORF3d-2 (gold). The status of the first half of ORF3d is obscured by the signal of ORF3c, but the increase in Frame 1 (ORF3a) reads and drop in proportion of Frame 2 (ORF3c) reads at the start of ORF3d is suggestive of an otherwise-unexpected biological change at this locus.

Figure 2—figure supplement 6.. Reading frame of ribosome profiling reads in the ORF3c/ORF3d region of ORF3a for individual treatments using a smaller sliding window size of 19 nt.

Details as in Figure 2C and Figure 2—figure supplement 5, except results are shown for a sliding window size of 19 nt (step size = 1 nt).

Figure 2—figure supplement 7.. Reading frame of ribosome profiling reads in the N/ORF9b/ORF9c region.

Details as in Figure 2C and Figure 2—figure supplement 5, except Frames 1, 2, and 3 correspond to the frames beginning at positions 1, 2, and 3 of the Wuhan-Hu-1 reference genome, respectively. Expression of ORF9b is evidenced by a peak in the fraction of reads mapping to Frame 3 (burgundy) at the appropriate region, as well as a peak in read depth (gray bars) immediately upstream (Figure 2A). Translational status of ORF9c is ambiguous, with minimal enrichment of reads in the appropriate frame (burgundy) along the gene.

Figure 2—figure supplement 8.. Reading frame of ribosome profiling reads in the full ORF3a gene.

Details as in Figure 2C and Figure 2—figure supplement 5. No clear evidence is observed for the expression of any of the ORF3b ORFs present in SARS-CoV-2, but the translational status of these ORFs requires further investigation.

Figure 2—figure supplement 9.. Reading frame of ribosome profiling reads in the ORF3c/ORF3d region of ORF3a for reads of length 29 nt.

Details as in Figure 2C and Figure 2—figure supplement 5, except results are shown for reads of length 29 nt, which are also indicative of frame but have slightly lower coverage (Figure 2—figure supplement 3) and exhibit frame-related periodicity. Evidence for ORF3c (Frame 2; burgundy) is more pronounced using this read length, while the peak for ORF3d (Frame 3; gold) occurs slightly later, overlapping ORF3d-2.

Figure 3.. SARS-CoV-2 protein sequence properties.

(A) Predicted densities of MHC class I-bound CD8+ T cell 9 amino acid (aa) epitopes (top) and MHC class II-bound CD4+ T cell 15 aa epitopes (bottom). Results for proteins encoded downstream of ORF1ab are shown. Mean numbers of predicted epitopes per residue (blue bars) are calculated as the number of epitopes overlapping each amino acid position divided by protein length. Error bars show 95% confidence intervals. Two sets of negative controls were also tested: (1) n = 1000 randomized peptides generated from each protein by randomly sampling its amino acids with replacement (dark gray bars), representing the result expected given amino acid content alone; and (2) short unannotated ORFs, the peptides encoded by n = 103 putatively nonfunctional ORFs present in the SARS-CoV-2 genome, representing the result expected for ORFs that have been evolving without functional constraint. For the short unannotated ORFs, the horizontal gray dotted line shows the mean number of epitopes per residue, and the gray-shaded region shows the 95% confidence interval (i.e. 2.5% to 97.5% quantiles). ORF3d, N, and ORF8 have the lowest MHC class I epitope densities; ORF3d, ORF3b, ORF7b, and ORF10 have the lowest MHC class II epitope densities. (B) Hydrophobicity profiles of amino acid sequences encoded by the three forward-strand reading frames of the ORF3a/ORF3c/ORF3d gene region, as calculated by the VOLPES server, using the unitless ‘Factor 1’ consensus hydrophobicity scale. Frame is reported using ORF3a as the reference, for example ss12 refers to the frame encoding ORF3d, for which codon position 2 overlaps codon position 1 of ORF3a (Figure 2C).

Figure 3—figure supplement 1.. Structural prediction for the ORF3d protein.

Independent computational modeling predictions of α-helices at the secondary (left inset, carried out in NetSurf v2) and tertiary structure levels (right inset, carried out in I-TASSER). Folding concordance with the closest protein structure match is shown (rightmost inset, aligned with TM-align). For explanation of shown metrics, see https://zhanglab.ccmb.med.umich.edu/I-TASSER. Chan et al., 2020 also predict a fold with α-helices (Raven Kok, pers. comm.).

Figure 3—figure supplement 2.. Correlations between hydrophobicity profiles of amino acid sequences encoded by the three forward-strand reading frames of the ORF3a region.

Vertical dotted gray lines denote gene start sites. (A) Hydrophobicity profiles of the amino acid sequences encoded by all three forward-strand reading frames of ORF3a. Details as in Figure 3B, except the entire region encoding ORF3a is shown. (B) Correlations (Spearman’s rank) between hydrophobicity profiles for the three forward-strand reading frames. Color denotes significance.

Figure 3—figure supplement 3.. Correlation between hydrophobicity profiles in the amino acid sequences encoded by all three forward-strand frames of ORF3a, by gene subregion.

Details as in Figure 3—figure supplement 2. ORF3a (left) refers to those codons of ORF3a that encode residues not overlapping any hypothesized overlapping genes. Codons counts refer to ORF3a, that is 21 codons of ORF3a (20 codons of ORF3d) are involved in a triple overlap (at least one nucleotide) with ORF3c.

Figure 4.. Amino acid variation in proteins encoded by genes overlapping ORF3a in viruses of the species Severe acute respiratory syndrome-related coronavirus.

(A) Amino acid sequences of ORF3c, ORF3d, and ORF3b, as encoded by SARS-CoV-2 (Wuhan-Hu-1; NCBI: NC_045512.2). Note that, while ORF3b encodes a protein of 154 amino acids (aa) in SARS-CoV (NCBI: NC_004718.3), a premature STOP codon in SARS-CoV-2 has resulted in an ORF encoding only 22 aa. (B) Amino acid alignments of ORF3c, ORF3d, and ORF3b (3b) show their sequence conservation. Black lines indicate STOP codons in ORF3d and ORF3b, showing their restricted taxonomic ranges. Intact ORF3d is restricted to SARS-CoV-2 and pangolin-CoV GX/P5L; however, note that ORF3d-2 (denoted 3d-2), a shorter isoform of ORF3d, could have a slightly wider taxonomic range if TTG or GTG are permitted as translation initiation codons. Full-length ORF3b (ORF3b region) is found throughout members of this virus species, but truncated early in most genomes outside of SARS-CoV (Supplementary file 1), with the shortest isoform (denoted 3b) found in SARS-CoV-2 and closely related viruses.

Figure 5.. Natural selection analysis of viral nucleotide differences at three hierarchical evolutionary levels.

Nucleotide differences in each virus gene were analyzed at three host levels: between-taxa divergence (d) among Severe acute respiratory syndrome-related coronavirus genomes infecting bat, human, and pangolin species; between-host diversity (π) for SARS-CoV-2 infecting human individuals (consensus-level); and within-host diversity (π) for SARS-CoV-2 infecting human individuals (deep sequencing). Each gene/level is shaded according to the ratio of mean nonsynonymous to synonymous differences per site to indicate purifying selection (d_N/d_S < 1 or π_N/π_S < 1; blue) or positive selection (d_N/d_S > 1 or π_N/π_S > 1; red). The extremely low ratio for ORF3c was artificially adjusted to allow the display of other ratios, and a Jukes-Cantor correction was applied to d_N and d_S values. Values range from a minimum of π_N/π_S = 0.04 (ORF3c, between-host; p=0.0410) to a maximum of 21.0 (ORF9c, between-host; p=0.126), where significance was evaluated using Z-tests of the hypothesis that d_N-d_S = 0 or π_N-π_S = 0 (10,000 bootstrap replicates, codon unit). The mean of all pairwise comparisons is shown for sequenced genomes only, i.e. no ancestral sequences were reconstructed or inferred. For each gene, sequences were only included in the between-species analysis if a complete, intact ORF (no STOPs) was present. Genes containing an overlapping gene (OLG) in a different frame were analyzed separately for non-OLG and OLG regions using SNPGenie and OLGenie, respectively. For ORF3b, only the region corresponding to the first ORF in SARS-CoV-2 (Table 1) was analyzed. The short overlap between ORF1a and ORF1b (nsp11 and nsp12) was excluded from the analysis. Error bars represent the standard error of mean pairwise differences. See Materials and methods for further details.

Figure 5—figure supplement 1.. SARS-CoV-2 between-host nucleotide diversity and allele frequencies as a function of date during the initial period of the COVID-19 pandemic.

Analyses refer to human SARS-CoV-2 consensus sequences obtained from GISAID. Results show sliding windows of 14 days (step size = 7 days) representing 13 time points since the first GISAID sample was collected (EPI_ISL_402123 on 12/24/2019). Regions with overlapping genes were excluded from estimates of π. Shaded regions show standard error of mean π_N, π_S, or π_N/π_S (10,000 bootstrap replicates, codon unit). (A) Nonsynonymous (π_N) and synonymous (π_S) nucleotide diversity. (B) π_N/π_S, where the horizontal dashed gray line denotes the π_N/π_S ratio expected under neutrality (1.0). (C) Diversity of sampling locations measured as entropy, defined as -∑p*ln(p), where p is the number of distinct (unique) locations or countries reported for a given window (Ewens and Grant, 2001). (D) Allele frequency trajectories of the ORFd-LOF mutation, colored by continent and limited to sufficient sample sizes (Supplementary file 1).

Figure 5—figure supplement 2.. Correlation between natural selection and gene expression.

A weak negative Spearman correlation between d_N/d_S (π_N/π_S) and gene expression is observed across all three evolutionary levels (between-taxa, between-host, and within-host). The x-axes of different panels differ by method for inferring gene expression, which are those shown in Figure 2B: (1) ‘Mass spec’ refers to log₁₀(iBAQ prop), where iBAQ prop is the proportion of the maximum intensity-based absolute quantification (iBAQ) observed; (2) ‘Upstream peak’ is the peak in ribosome profiling read depth upstream of the gene start site, measured for reads of length 29-31nt (Figure 2A); and (3) ‘In-frame mean’ is the mean read depth at Codon pos 1, measured for reads of length 30nt (Figure 2B). A given gene (color) has only one selection value (y-axis) per evolutionary level, and only one gene expression value (x-axis) per method for determining expression. Selection was calculated as d_N/d_S (π_N/π_S) except for the overlapping gene ORF9b, where the OLG-appropriate measure d_NN/d_NS (π_NN/π_NS) was used (see Materials and methods).

Figure 6.. Between-taxa sliding window analysis of natural selection on overlapping frames of ORF3a.

Pairwise sliding window analysis (window size = 50 codons; step size = 1 codon) of selection across members of the species Severe acute respiratory syndrome-related coronavirus. OLG-appropriate d_N/d_S values were computed using OLGenie (Nelson et al., 2020a), a method that is conservative (non-conservative) for detecting purifying (positive) selection, and a Jukes-Cantor correction for multiple hits was employed. On the left-hand side, each genome is compared to SARS-CoV-2 in the ss12 reading frame of ORF3a, which contains ORF3d (Table 1 and Figure 2C). This frame shows evidence for purifying selection specific to the ORF3d region that is limited to the comparison with pangolin-CoV GX/P5L. On the right-hand side, this analysis is repeated for the ss13 reading frame of ORF3a, which contains ORF3c and ORF3b (Table 1 and Figure 2C), this time with respect to SARS-CoV, where full-length ORF3b is functional. This frame shows constraint across much of this gene and virus species, and ORF3c in particular is deeply conserved.

Figure 6—figure supplement 1.. Between-taxa sliding window of genes overlapping N.

Pairwise OLGenie analysis of the N gene’s ss13 reading frame across members of the species Severe acute respiratory syndrome-related coronavirus. Details as in Figure 6, except only the ss13 frame is shown.

Figure 7.. Pandemic spread of the EP+1 haplotype and the hitchhiking of ORF3d-LOF.

The mutational path leading to EP+1+LOF is shown in the upper panel. Cumulative frequencies of haplotypes in samples from Germany and five other countries with the most abundant sequence data are shown in the lower panel. Countries are grouped into early founder (left) and late founder (right) based on the presence or absence of SARS-CoV-2 samples from January, respectively. In the early founder group, EP–3 (gray) is observed earlier than other haplotypes in France and the US, and EP (red) is observed early in Germany, giving them the advantage of a founder effect. However, neither EP nor EP–3 dominate later spread; instead, EP+1 (yellow) and EP+1+LOF (blue) increase much faster despite their later appearance in these countries. In the late founder group, multiple haplotypes appear at approximately the same time, but EP–3 and EP spread more slowly. The green dashed line shows the combined frequencies of EP+1 and EP+1+LOF (yellow and blue, respectively). Note that EP–1 is never observed in our dataset.

Figure 8.. High-frequency within-host mutations.

For each site, mutations that occur in more than 2.5% of samples are shown, limiting to samples where the major or fixed allele matches Wuhan-Hu-1 at that site. The y-axis shows the proportion of such samples having the indicated minor (derived) allele. Each locus has up to three possible single-nucleotide derived alleles compared to the reference background. Open circles (black outlines) show the proportion of samples having any of the three possible derived alleles ('All'), while solid circles (color fill) show the proportion of samples having a specific derived allele (equivalent to 'All' if only one variant is observed). For most sites, only one derived mutation type (e.g. C→U) is observed across all samples. Precluding co-infection by multiple genotypes and sequencing errors, derived mutations occurring in more than one sample (y-axis) must be identical by state but not descent (i.e. recurrent). Genome positions are plotted on the x-axis, with distinct genes shown in different colors and overlapping genes shown as black blocks within reference genes. Nonsynonymous and nonsense mutations (together denoted 'Nonsyn') are indicated with a red dot. Source data are available in the supplementary material.

Figure 8—figure supplement 1.. Recurrent nonsynonymous mutations observed in multiple human hosts.

Histogram of the derived allele frequency (DAF) of the four most common recurrent protein-coding mutations across all within-host samples (bin width = 0.01). Mutation A404U introduces a premature STOP codon in nsp1, whereas the remainder are nonsynonymous.