The Missing Expression Level–Evolutionary Rate Anticorrelation in Viruses Does Not Support Protein Function as a Main Constraint on Sequence Evolution

Figures

• 

  Fig. 1. The virus is an excellent system to distinguish two hypotheses explaining the E–R anticorrelation in cellular organisms. (A) Selective constraints imposed on the evolution of a protein sequence are hypothesized to entail the maintenance of protein function, the avoidance of cytotoxicity caused by protein misfolding and misinteraction, or both. Protein sequence evolution of cellular organisms may be constrained by both mechanisms, but the selection against cytotoxicity is relaxed for viral proteins (red dashed arrow). (B) The E–R anticorrelation has been observed in all three domains of cellular organisms. (C) Viruses can be used to test the relative importance of the two hypotheses explaining the E–R anticorrelation. If the E–R anticorrelation remains strong among viral ORFs, the cytotoxicity avoidance hypothesis is not supported (left), because the avoidance of cytotoxicity appears to be unnecessary for creating an E–R anticorrelation. In contrast, if the E–R anticorrelation is missing in viruses, the function maintenance hypothesis is not supported (right), because the maintenance of protein function per se is insufficient to create an E–R anticorrelation.

• 

  Fig. 2. Detection of the E–R anticorrelation in SARS-CoV-2. (A) A schematic shows how the protein evolutionary rate is estimated for Vero cell and SARS-CoV-2. The red stars on the right panel denote amino acid alternations. (B) A schematic shows how expression level is estimated for Vero cell and SARS-CoV-2. On the left penal, the expression level for each Vero cell gene is estimated as the normalized reads abundance aligned to it. On the right penal, the genome and two out of the eight subgenomes of SARS-CoV-2 are shown as examples. The gray lines denote the regions skipped by the leader-to-body fusion during the discontinuous transcription in coronaviruses. Note that the coding sequence (CDS) of a subgenome could be the 3′-UTR of another subgenome. Therefore, the expression level for each ORF in SARS-CoV-2 was estimated as the number of lead-containing reads of the corresponding subgenome. (C) The E–R correlation in the Vero cell (in blue) and in SARS-CoV-2 (in red). For Vero cells, expressed protein-coding genes are split into ten bins according to their expression levels: (−∞, 0.5), [0.5, 1.5), [1.5, 2.5), [2.5, 3.5), [3.5, 4.5), [4.5, 5.5), [5.5, 6.5), [6.5, 7.5), [7.5, 8.5), [8.5,…

• 

  Fig. 3. The E–R anticorrelation remains absent in SARS-CoV-2 when nsps are considered individually. (A) A schematic shows nsps derived from the proteolysis of pp1ab, which is encoded by ORF1ab. The black arrows denote the cleavage sites. (B) The absence of the E–R anticorrelation among 23 proteins encoded by the SARS-CoV-2 genome. (C) Similar to (B), but three proteins (S, E, and M) known to be packaged for secretion at the endoplasmic reticulum are excluded.

• 

  Fig. 4. The E–R anticorrelation remains absent in SARS-CoV-2 when the protein evolutionary rates in animal hosts are used. (A) The phylogenetic tree (red branches) from which the protein evolutionary rate (dN/dS) of individual SARS-CoV-2 ORFs was estimated. The scale bar represents 0.01 substitutions per nucleotide. (B) The absence of an anticorrelation between expression level, estimated by the number of leader-containing reads, and the protein evolutionary rate, estimated as the dN/dS ratio among four coronaviruses shown in (A).

• 

  Fig. 5. The absence of E–R anticorrelation in nine other virus species. (A–I) The E–R correlation in individual virus species. The number of variants used for the estimation of protein evolutionary rate for each virus species is shown in (J). For Vaccinia virus (H), a constant number (one) has been added to all expression data prior to applying the logarithm transform, to avoid the logarithm of zero. (J) A meta-analysis of E–R correlation for ten virus species. Spearman’s correlation coefficient for individual viruses is shown. The combined correlation coefficient is given by meta-analysis of correlation coefficients (metacor function in the “meta” package of R, which combines correlation coefficients from ten individual viruses into one pooled correlation estimate). Here, the random-effects model with Sidik–Jonkman estimator is used for the estimation of between-study heterogeneity τ2, and z-transformed correlations are used for the meta-analysis. The results of the fixed-effect model is presented because the test of heterogeneity gave a P value = 0.16 (I2 = 21.5%). The effective sample size N is also shown for the meta-analysis.

• 

  Fig. 6. The absence of E–R anticorrelation when the protein evolutionary rates are estimated as the dN/dS ratio of orthologous genes. (A–C) The E–R correlation coefficients in individual virus groups. The expression level ORFs of SARS-CoV-2 (A), human cytomegalovirus (B), or Vaccinia virus (C) were used; the protein evolutionary rates were estimated as the dN/dS ratio between SARS-CoV-2 and SARS-CoV (A), between human and chimpanzee cytomegalovirus (B), or among five Orthopoxvirus species (C). The scale bar of the phylogenetic tree in (C) represents 0.01 substitutions per nucleotide. For Vaccinia virus, a constant number (one) has been added to all expression data prior to applying the logarithm transform, to avoid the logarithm of zero. (D) Similar to figure 5J, showing a meta-analysis of E–R correlation for three groups of viruses in figure 6A–C. The results of the fixed-effect model is presented because the test of heterogeneity gave a P value = 0.88 (I2 = 0.0%). The effective sample size N is also shown for the meta-analysis.

• 

  Fig. 7. The presence of E–R anticorrelation among HERV ORFs. (A) The selection against cytotoxicity is restored when a viral genome is integrated into a host genome. (B) The syntenic reciprocal best BLAST hits endogenous retrovirus pairs were identified as human–chimpanzee orthologous pairs. (C) A schematic shows the estimation of the protein evolutionary rate for each HERV ORF. (D) The E–R correlation in HERVs (ρ = −0.20, P = 1.1 × 10−13, N = 1,316, Spearman’s correlation). The expression levels of HERV ORFs are estimated from RNA-seq data for human peripheral blood mononuclear cells from six healthy individuals. The plot shows the expression level in healthy #1 (replicate #1). Replicates #2–#6 are shown in supplementary figure S3, Supplementary Material online. The HERV ORFs are split into six bins according to their expression levels: (−∞, −9.5), [−9.5, −8.5), [−8.5, −7.5), [−7.5, −6.5), [−6.5, −5.5), and [−5.5, +∞). The mean (±standard errors) of the evolutionary rate is shown for each bin. Spearman’s correlation coefficient is calculated from the unbinned data. (E) The E–R anticorrelation among HERVs (indicated by the red arrow) was not significantly different from that among human…