Single-cell immune profiling of mouse liver aging reveals Cxcl2+ macrophages recruit neutrophils to aggravate liver injury.

BACKGROUND AND AIMS: Immune cells play a crucial role in liver aging. However, the impact of dynamic changes in the local immune microenvironment on age-related liver injury remains poorly understood. We aimed to characterize intrahepatic immune cells at different ages to investigate key mechanisms associated with liver aging. APPROACH AND RESULTS: We carried out single-cell RNA sequencing on mouse liver tissues at 4 different ages, namely, the newborn, suckling, young, and aged stages. The transcriptomic landscape, cellular classification, and intercellular communication were analyzed. We confirmed the findings by multiplex immunofluorescence staining, flow cytometry, in vitro functional experiments, and chimeric animal models. Nine subsets of 89,542 immune cells with unique properties were identified, of which Cxcl2+ macrophages within the monocyte/macrophage subset were preferentially enriched in the aged liver. Cxcl2+ macrophages presented a senescence-associated secretory phenotype and recruited neutrophils to the aged liver through the CXCL2-CXCR2 axis. Through the secretion of IL-1ß and TNF-α, Cxcl2+ macrophages stimulated neutrophil extracellular traps formation. Targeting the CXCL2-CXCR2 axis limited the neutrophils migration toward the liver and attenuated age-related liver injury. Moreover, the relationship between Cxcl2+ macrophages and neutrophils in age-related liver injury was further validated by human liver transplantation samples. CONCLUSIONS: This in-depth study illustrates that the mechanism of Cxcl2+ macrophage-driven neutrophil activation involves the CXCL2-CXCR2 axis and provides a potential therapeutic strategy for age-related liver injury.

Phylogenetic Comparative Analysis of Single-Cell Transcriptomes Reveals Constrained Accumulation of Gene Expression Heterogeneity during Clonal Expansion.

In the same way that a phylogeny summarizes the evolutionary history of species, a cell lineage tree describes the process of clonal expansion, in which gene expression differences between cells naturally accrue as a result of stochastic partitioning and imperfect expression control. How is functional homeostasis, a key factor in the biological function of any population of cells, maintained in the face of such continuous accumulation of transcriptomic heterogeneity remains largely unresolved. To answer this question, we experimentally determined the single-cell transcriptomes and lineage relationships of up to 50% cells in single-HEK293-seeded colonies. Phylogenetic comparative analyses of the single-cell transcriptomes on the cell lineage tree revealed three lines of evidence for the constrained accumulation of transcriptome heterogeneity among cells, including rapid saturation of transcriptomic heterogeneity upon four cell divisions, reduced expression differences within subtrees closer to expression boundaries, and cofluctuations among genes. Our analyses showcased the applicability of phylogenetic comparative methods in cell lineage trees, demonstrated the constrained accumulation of transcriptomic heterogeneity, and provided novel insight into the functional homeostasis of cell populations.

A Prediction Model for Clinical Outcomes of COVID-19 Hospitalized Patients: Construction and Accuracy Assessment.

BACKGROUND: The coronavirus disease 2019 (COVID-19) pandemic spread rapidly with considerable morbidity nationwide since China's liberalization in December 2022. Our work has focused on identifying different predictive factors from the laboratory examination of critically ill patients, and forecasting the unfavorable outcome of critically ill patients with COVID-19 through a combined diagnosis of biological markers. METHODS: We conducted a retrospective study at the Department of First Affiliated Hospital of Wenzhou Medical University, China, from December 24, 2022, to January 10, 2023, where 434 critically ill patients who met the inclusion criteria were involved. Machine analysis was employed to search for the parameters with the highest predictive value to calculate COVID-19 mortality by exploiting 66 typical laboratory results. RESULTS: Combined diagnosis of serum albumin (ALB), lactate dehydrogenase (LDH), direct bilirubin (Dbil), ferritin, pulse oxygen saturation (SpO2), and neutrophil count (NEUT#) was evaluated, and the result with the highest predictive value (NEUT#) was selected as the predictor for COVID-19 mortality with a sensitivity of 89.2% and a specificity of 77.4%. CONCLUSIONS: The increased levels of LDH, Dbil, ferritin, and NEUT#, along with lowered ALB and SpO2 levels are the most decisive variables for forecasting the mortality for COVID-19 according to our machine-learning-based model. The combined diagnosis could be used to improve further diagnostic performance.

Prediction of Antibiotic Resistance Evolution by Growth Measurement of All Proximal Mutants of Beta-Lactamase.

The antibiotic resistance crisis continues to threaten human health. Better predictions of the evolution of antibiotic resistance genes could contribute to the design of more sustainable treatment strategies. However, comprehensive prediction of antibiotic resistance gene evolution via laboratory approaches remains challenging. By combining site-specific integration and high-throughput sequencing, we quantified relative growth under the respective selection of cefotaxime or ceftazidime selection in ∼23,000 Escherichia coli MG1655 strains that each carried a unique, single-copy variant of the extended-spectrum ß-lactamase gene blaCTX-M-14 at the chromosomal att HK022 site. Significant synergistic pleiotropy was observed within four subgenic regions, suggesting key regions for the evolution of resistance to both antibiotics. Moreover, we propose PEARP and PEARR, two deep-learning models with strong clinical correlations, for the prospective and retrospective prediction of blaCTX-M-14 evolution, respectively. Single to quintuple mutations of blaCTX-M-14 predicted to confer resistance by PEARP were significantly enriched among the clinical isolates harboring blaCTX-M-14 variants, and the PEARR scores matched the minimal inhibitory concentrations obtained for the 31 intermediates in all hypothetical trajectories. Altogether, we conclude that the measurement of local fitness landscape enables prediction of the evolutionary trajectories of antibiotic resistance genes, which could be useful for a broad range of clinical applications, from resistance prediction to designing novel treatment strategies.

Risk factors of non-sentinel lymph node metastasis in breast cancer with 1-2 sentinel lymph node macrometastases underwent total mastectomy: a case-control study.

BACKGROUND: The randomized trials which include ACOSOG Z0011 and IBCSG 23-01 had found that the survival rates were not different in patients with cT1/2N0 and 1-2 sentinel lymph node (SLN)-positive, macro/micrometastases who underwent breast-conserving therapy, and micrometastases who underwent total mastectomy (TM), when axillary lymph node dissection (ALND) was omitted. However, for patients with cT1/2N0 and 1-2 SLN macrometastases who underwent TM; there was still insufficient evidence from clinical studies to support whether ALND can be exempted. This study aimed to investigate the risk factors of non-sentinel lymph node (nSLN) metastasis in breast cancer patients with 1-2 SLN macrometastases undergoing TM. METHODS: The clinicopathological data of 1491 breast cancer patients who underwent TM and SLNB from January 2017 to February 2022 were retrospectively analyzed. Univariate and multivariate analyses were performed to analyze the risk factors for nSLN metastasis. RESULTS: A total of 273 patients with 1-2 SLN macrometastases who underwent TM were enrolled. Postoperative pathological data showed that 35.2% patients had nSLN metastasis. The results of multivariate analysis indicated that tumor size (TS) (P = 0.002; OR: 1.051; 95% CI: 1.019-1.084) and ratio of SLN macrometastases (P = 0.0001; OR: 12.597: 95% CI: 4.302-36.890) were the independent risk factors for nSLN metastasis in breast cancer patients with 1-2 SLN macrometastases that underwent TM. The ROC curve analysis suggested that when TS ≤22 mm and ratio of SLN macrometastases ≤0.33, the incidence of nSLN metastasis could be reduced to 17.1%. CONCLUSIONS: The breast cancer patients with cT1/2N0 stage, undergoing TM and 1-2 SLN macrometastases, when the TS ≤22 mm and macrometastatic SLN does not exceed 1/3 of the total number of detected SLN, the incidence of nSLN metastasis is significantly reduced, but whether ALND can be exempted needs further exploration.

Bidirectional Genetic Control of Phenotypic Heterogeneity and Its Implication for Cancer Drug Resistance.

Negative genetic regulators of phenotypic heterogeneity, or phenotypic capacitors/stabilizers, elevate population average fitness by limiting deviation from the optimal phenotype and increase the efficacy of natural selection by enhancing the phenotypic differences among genotypes. Stabilizers can presumably be switched off to release phenotypic heterogeneity in the face of extreme or fluctuating environments to ensure population survival. This task could, however, also be achieved by positive genetic regulators of phenotypic heterogeneity, or "phenotypic diversifiers," as shown by recently reported evidence that a bacterial divisome factor enhances antibiotic resistance. We hypothesized that such active creation of phenotypic heterogeneity by diversifiers, which is functionally independent of stabilizers, is more common than previously recognized. Using morphological phenotypic data from 4,718 single-gene knockout strains of Saccharomyces cerevisiae, we systematically identified 324 stabilizers and 160 diversifiers and constructed a bipartite network between these genes and the morphological traits they control. Further analyses showed that, compared with stabilizers, diversifiers tended to be weaker and more promiscuous (regulating more traits) regulators targeting traits unrelated to fitness. Moreover, there is a general division of labor between stabilizers and diversifiers. Finally, by incorporating NCI-60 human cancer cell line anticancer drug screening data, we found that human one-to-one orthologs of yeast diversifiers/stabilizers likely regulate the anticancer drug resistance of human cancer cell lines, suggesting that these orthologs are potential targets for auxiliary treatments. Our study therefore highlights stabilizers and diversifiers as the genetic regulators for the bidirectional control of phenotypic heterogeneity as well as their distinct evolutionary roles and functional independence.

Analysis of Risk Factors of Early Mortality After Pericardiectomy For Constrictive Pericarditis.

BACKGROUND: We aimed to investigate risk factors of early mortality following pericardiectomy. METHODS: This was a retrospective study of patients undergoing pericardiectomy between January 1994 and May 2021 at The People's Hospital of Guangxi Zhuang Autonomous Region, Ruikang Hospital Affiliated to Guangxi University of Chinese Medicine, and The People's Hospital of Liuzhou City. RESULTS: This study included 826 patients, who were divided into two groups: group with operative deaths (N = 66) and group without operative deaths (N = 760). There were 66 operative deaths (66/826, 8.0%). The causes of operative deaths were multiorgan failure (86/826, 10.4%). Preoperative CVP (P < 0.001), chest drainage (P < 0.001), surgical duration (P < 0.001), fluid balance postoperative day D2 (P < 0.001), and tuberculosis pericarditis (P = 0.001) in group with operative deaths were significantly higher than those in group without operative deaths. Univariate and multivariate analyses showed that factors associated with operative deaths include male (P < 0.001), age (P < 0.001), ICU retention time (P < 0.001), postoperative hospitalization time (P < 0.001), preoperative central venous pressure (P = 0.018), postoperative central venous pressure (P < 0.001), D0 fluid balance (P < 0.001), D2 fluid balance (P < 0.001), postoperative chest drainage (P = 0.029), surgical duration (P = 0.003), serum creatinine baseline (P = 0.002), serum creatinine 24h after surgery (P < 0.001), serum creatinine 48h after surgery (P < 0.001), blood lactate (P < 0.001), and tuberculosis pericarditis (P = 0.033). CONCLUSION: In our study, incomplete pericardial dissection, fluid overload, and tuberculosis pericarditis are associated with operative deaths following pericardiectomy.

Incomplete Pericardial Dissection, Fluid Overload, Delayed Diagnosis And Treatment, And Tuberculosis Pericarditis Are Associated With Low Cardiac Output Syndrome Following Pericardiectomy.

BACKGROUND: We aimed to investigate risk factors of LCOS following pericardiectomy. METHODS: This was a retrospective study of patients undergoing pericardiectomy at three hospitals between January 1994 and May 2021. RESULTS: A total of 826 patients were divided into two groups: group with LCOS (N = 126) and group without LCOS (N = 700). The incidence of postoperative LCOS was 15.3%. There were 66 operative deaths (8.0%). Univariable and multivariable analyses showed that factors are associated with LCOS, including postoperative LVEDD (P < 0.001), preoperative LVEDD (P < 0.001), time between symptoms and surgery (P < 0.001), thickness of pericardium (P < 0.001), intubation time (P = 0.002), hospitalized time postoperative (P < 0.001), preoperative central venous pressure (P = 0.016), postoperative central venous pressure (P = 0.034), D0 fluid balance (P = 0.019), D2 fluid balance (P = 0.017), postoperative chest drainage (P < 0.001), surgical duration (P < 0.001), bleeding during operation (P = 0.001), serum creatinine 24h after surgery (P < 0.001), serum creatinine 48h after surgery (P = 0.017), fresh-frozen plasma (P = 0.005), packed red cells (P = 0.006), and tuberculosis pericarditis (P = 0.026). CONCLUSION: In our study, incomplete pericardial dissection, fluid overload, delayed diagnosis and treatment, and tuberculosis pericarditis are associated with LCOS following pericardiectomy.

Determinants of the rate of protein sequence evolution.

The rate and mechanism of protein sequence evolution have been central questions in evolutionary biology since the 1960s. Although the rate of protein sequence evolution depends primarily on the level of functional constraint, exactly what determines functional constraint has remained unclear. The increasing availability of genomic data has enabled much needed empirical examinations on the nature of functional constraint. These studies found that the evolutionary rate of a protein is predominantly influenced by its expression level rather than functional importance. A combination of theoretical and empirical analyses has identified multiple mechanisms behind these observations and demonstrated a prominent role in protein evolution of selection against errors in molecular and cellular processes.

Persistent viral RNA positivity during the recovery period of a patient with SARS-CoV-2 infection.

As an emerging infectious disease, the clinical course and virological course of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) infection remain to be further investigated. In this case report, we described a case of SARS-CoV-2 infection with the clinical course for more than 2 months. This patient had recovered from pneumonia after treatment. The viral RNA of throat swabs became negative and the viral-specific antibodies were produced during the recovery period. However, the viral RNA reappeared and additionally persisted in throat swabs for more than 40 days. In addition, the viral RNA was detected in multiple types of specimens with extremely high titers in the saliva. In conclusion, these findings indicate that SARS-CoV-2 can cause a long clinical course. The coexistence of viral RNA and viral-specific antibodies may imply an immune evasion of SARS-CoV-2 from the host's immune system.

Risk stratification of cervical lesions using capture sequencing and machine learning method based on HPV and human integrated genomic profiles.

From initial human papillomavirus (HPV) infection and precursor stages, the development of cervical cancer takes decades. High-sensitivity HPV DNA testing is currently recommended as primary screening method for cervical cancer, whereas better triage methodologies are encouraged to provide accurate risk management for HPV-positive women. Given that virus-driven genomic variation accumulates during cervical carcinogenesis, we designed a 39 Mb custom capture panel targeting 17 HPV types and 522 mutant genes related to cervical cancer. Using capture-based next-generation sequencing, HPV integration status, somatic mutation and copy number variation were analyzed on 34 paired samples, including 10 cases of HPV infection (HPV+), 10 cases of cervical intraepithelial neoplasia (CIN) grade and 14 cases of CIN2+ (CIN2: n = 1; CIN2-3: n = 3; CIN3: n = 9; squamous cell carcinoma: n = 1). Finally, the machine learning algorithm (Random Forest) was applied to build the risk stratification model for cervical precursor lesions based on CIN2+ enriched biomarkers. Generally, HPV integration events (11 in HPV+, 25 in CIN1 and 56 in CIN2+), non-synonymous mutations (2 in CIN1, 12 in CIN2+) and copy number variations (19.1 in HPV+, 29.4 in CIN1 and 127 in CIN2+) increased from HPV+ to CIN2+. Interestingly, 'common' deletion of mitochondrial chromosome was significantly observed in CIN2+ (P = 0.009). Together, CIN2+ enriched biomarkers, classified as HPV information, mutation, amplification, deletion and mitochondrial change, successfully predicted CIN2+ with average accuracy probability score of 0.814, and amplification and deletion ranked as the most important features. Our custom capture sequencing combined with machine learning method effectively stratified the risk of cervical lesions and provided valuable integrated triage strategies.

Dosage sensitivity of X-linked genes in human embryonic single cells.

BACKGROUND: During the evolution of mammalian sex chromosomes, the degeneration of Y-linked homologs has led to a dosage imbalance between X-linked and autosomal genes. The evolutionary resolution to such dosage imbalance, as hypothesized by Susumu Ohno fifty years ago, should be doubling the expression of X-linked genes. Recent studies have nevertheless shown that the X to autosome expression ratio equals ~ 1 in haploid human parthenogenetic embryonic stem (pES) cells and ~ 0.5 in diploid pES cells, suggesting no doubled expression for X-linked genes and refuting Ohno's hypothesis. RESULTS: Here, by reanalyzing an RNA-seq-based single-cell transcriptome dataset of human embryos, we found that from the 8-cell stage until the time-point just prior to implantation, the expression levels of X-linked genes are not two-fold upregulated in male cells and gradually decrease from two-fold in female cells. Additional analyses of gene expression noise further suggest that the dosage sensitivity of X-linked genes is weaker than that of autosomal genes in differentiated female cells, which contradicts a key assumption in Ohno's hypothesis, that most X-linked genes are dosage sensitive. Moreover, the dosage-sensitive housekeeping genes are preferentially located on autosomes, implying selection against X-linkage for dosage-sensitive genes. CONCLUSIONS: We observed dosage imbalance between X-linked and autosomal genes, as well as relatively high expression noise from X-linked genes. These results collectively suggest that X-linked genes are less dosage sensitive than autosomal genes, putting one primary assumption of Ohno's hypothesis in question.

Nascent RNA folding mitigates transcription-associated mutagenesis.

Transcription is mutagenic, in part because the R-loop formed by the binding of the nascent RNA with its DNA template exposes the nontemplate DNA strand to mutagens and primes unscheduled error-prone DNA synthesis. We hypothesize that strong folding of nascent RNA weakens R-loops and hence decreases mutagenesis. By a yeast forward mutation assay, we show that strengthening RNA folding and reducing R-loop formation by synonymous changes in a reporter gene can lower mutation rate by >80%. This effect is diminished after the overexpression of the gene encoding RNase H1 that degrades the RNA in a DNA-RNA hybrid, indicating that the effect is R-loop-dependent. Analysis of genomic data of yeast mutation accumulation lines and human neutral polymorphisms confirms the generality of these findings. This mechanism for local protection of genome integrity is of special importance to highly expressed genes because of their frequent transcription and strong RNA folding, the latter also improves translational fidelity. As a result, strengthening RNA folding simultaneously curtails genotypic and phenotypic mutations.

Intra and Interspecific Variations of Gene Expression Levels in Yeast Are Largely Neutral: (Nei Lecture, SMBE 2016, Gold Coast).

It is commonly, although not universally, accepted that most intra and interspecific genome sequence variations are more or less neutral, whereas a large fraction of organism-level phenotypic variations are adaptive. Gene expression levels are molecular phenotypes that bridge the gap between genotypes and corresponding organism-level phenotypes. Yet, it is unknown whether natural variations in gene expression levels are mostly neutral or adaptive. Here we address this fundamental question by genome-wide profiling and comparison of gene expression levels in nine yeast strains belonging to three closely related Saccharomyces species and originating from five different ecological environments. We find that the transcriptome-based clustering of the nine strains approximates the genome sequence-based phylogeny irrespective of their ecological environments. Remarkably, only ∼0.5% of genes exhibit similar expression levels among strains from a common ecological environment, no greater than that among strains with comparable phylogenetic relationships but different environments. These and other observations strongly suggest that most intra and interspecific variations in yeast gene expression levels result from the accumulation of random mutations rather than environmental adaptations. This finding has profound implications for understanding the driving force of gene expression evolution, genetic basis of phenotypic adaptation, and general role of stochasticity in evolution.

Deciphering the Genic Basis of Yeast Fitness Variation by Simultaneous Forward and Reverse Genetics.

The budding yeast Saccharomyces cerevisiae is the best studied eukaryote in molecular and cell biology, but its utility for understanding the genetic basis of phenotypic variation in natural populations is limited by inefficient association mapping due to strong and complex population structure. To overcome this challenge, we generated genome sequences for 85 strains and performed a comprehensive population genomic survey of a total of 190 diverse strains. We identified considerable variation in population structure among chromosomes and identified 181 genes that are absent from the reference genome. Many of these nonreference genes are expressed and we functionally confirmed that two of these genes confer increased resistance to antifungals. Next, we simultaneously measured the growth rates of over 4,500 laboratory strains, each of which lacks a nonessential gene, and 81 natural strains across multiple environments using unique DNA barcode present in each strain. By combining the genome-wide reverse genetic information gained from the gene deletion strains with a genome-wide association analysis from the natural strains, we identified genomic regions associated with fitness variation in natural populations. To experimentally validate a subset of these associations, we used reciprocal hemizygosity tests, finding that while the combined forward and reverse genetic approaches can identify a single causal gene, the phenotypic consequences of natural genetic variation often follow a complicated pattern. The resources and approach provided outline an efficient and reliable route to association mapping in yeast and significantly enhance its value as a model for understanding the genetic mechanisms underlying phenotypic variation and evolution in natural populations.

Codon-by-codon modulation of translational speed and accuracy via mRNA folding.

Rapid cell growth demands fast protein translational elongation to alleviate ribosome shortage. However, speedy elongation undermines translational accuracy because of a mechanistic tradeoff. Here we provide genomic evidence in budding yeast and mouse embryonic stem cells that the efficiency-accuracy conflict is alleviated by slowing down the elongation at structurally or functionally important residues to ensure their translational accuracies while sacrificing the accuracy for speed at other residues. Our computational analysis in yeast with codon resolution suggests that mRNA secondary structures serve as elongation brakes to control the speed and hence the fidelity of protein translation. The position-specific effect of mRNA folding on translational accuracy is further demonstrated experimentally by swapping synonymous codons in a yeast transgene. Our findings explain why highly expressed genes tend to have strong mRNA folding, slow translational elongation, and conserved protein sequences. The exquisite codon-by-codon translational modulation uncovered here is a testament to the power of natural selection in mitigating efficiency-accuracy conflicts, which are prevalent in biology.

Determinative developmental cell lineages are robust to cell deaths.

All forms of life are confronted with environmental and genetic perturbations, making phenotypic robustness an important characteristic of life. Although development has long been viewed as a key component of phenotypic robustness, the underlying mechanism is unclear. Here we report that the determinative developmental cell lineages of two protostomes and one deuterostome are structured such that the resulting cellular compositions of the organisms are only modestly affected by cell deaths. Several features of the cell lineages, including their shallowness, topology, early ontogenic appearances of rare cells, and non-clonality of most cell types, underlie the robustness. Simple simulations of cell lineage evolution demonstrate the possibility that the observed robustness arose as an adaptation in the face of random cell deaths in development. These results reveal general organizing principles of determinative developmental cell lineages and a conceptually new mechanism of phenotypic robustness, both of which have important implications for development and evolution.

Serum viral duplex-linear DNA proportion increases with the progression of liver disease in patients infected with HBV.

OBJECTIVE: HBV has two forms of genomic DNA, relaxed-circular DNA (rcDNA) and duplex-linear DNA (dlDNA). Compared to rcDNA, dlDNA has been demonstrated to integrate more frequently into host cellular chromosomes, which may have oncogenic consequences. However, the dlDNA proportion relative to total HBV DNA and its clinical significance in patients remain to be investigated. DESIGN: Based on the structural difference between rcDNA and dlDNA, we developed a peptide nucleic acid (PNA)-mediated quantitative real-time PCR (qPCR) clamping assay to measure the proportions of dlDNA in total HBV DNA in sera obtained from patients with chronic hepatitis B (CHB), liver cirrhosis (LC) or LC-developed hepatocellular carcinoma (HCC). The factors that influence the proportion of dlDNA were also investigated. RESULTS: The average dlDNA proportion was approximately 7% in the sera of chronic HBV-infected patients and was elevated in CHB patients with abnormal levels of alanine aminotransferase. The sera dlDNA proportions increased to approximately 14% and 20% in the patients with LC and HCC, respectively. Interferon-α treatment slightly increased the dlDNA proportion in the responders; and nucleotide analogue therapy spuriously elevated the proportion. Moreover, treatment of human hepatoma cells supporting HBV replication with inflammatory cytokines significantly altered the dlDNA proportion in vitro. CONCLUSIONS: Using a novel PNA-mediated qPCR clamping assay, we first showed that serum dlDNA proportions progressively increased during the development of HBV-related liver diseases. The dlDNA proportion can be regulated by inflammatory cytokines, suggesting an association among inflammation, increased production of HBV dlDNA and development of HCC.

Human long noncoding RNAs are substantially less folded than messenger RNAs.

Long noncoding RNAs (lncRNAs) do not code for proteins but function as RNAs. Because the functions of an RNA rely on either its sequence or secondary structure, lncRNAs should be folded at least as strongly as messenger RNAs (mRNAs), which serve as messengers for translation and are generally thought to lack secondary structure-dependent RNA-level functions. Contrary to this prediction, analysis of genome-wide experimental data of human RNA folding reveals that lncRNAs are substantially less folded than mRNAs even after the control of expression level and GC% (percentage of guanines and cytosines), although both lncRNAs and mRNAs are more strongly folded than expected by chance. In contrast to mRNAs, lncRNAs show neither the positive correlation between folding strength and expression level nor the negative correlation between folding strength and evolutionary rate. These and other results support that although RNA folding undoubtedly plays a role in RNA biology it is also important in translation and/or protein biology.

Differential requirements for mRNA folding partially explain why highly expressed proteins evolve slowly.

The cause of the tremendous among-protein variation in the rate of sequence evolution is a central subject of molecular evolution. Expression level has been identified as a leading determinant of this variation among genes encoded in the same genome, but the underlying mechanisms are not fully understood. We here propose and demonstrate that a requirement for stronger folding of more abundant mRNAs results in slower evolution of more highly expressed genes and proteins. Specifically, we show that: (i) the higher the expression level of a gene, the greater the selective pressure for its mRNA to fold; (ii) random mutations are more likely to decrease mRNA folding when occurring in highly expressed genes than in lowly expressed genes; and (iii) amino acid substitution rate is negatively correlated with mRNA folding strength, with or without the control of expression level. Furthermore, synonymous (d(S)) and nonsynonymous (d(N)) nucleotide substitution rates are both negatively correlated with mRNA folding strength. However, counterintuitively, d(S) and d(N) are differentially constrained by selection for mRNA folding, resulting in a significant correlation between mRNA folding strength and d(N)/d(S), even when gene expression level is controlled. The direction and magnitude of this correlation is determined primarily by the G+C frequency at third codon positions. Together, these findings explain why highly expressed genes evolve slowly, demonstrate a major role of natural selection at the mRNA level in constraining protein evolution, and reveal a previously unrecognized and unexpected form of nonprotein-level selection that impacts d(N)/d(S).