Comparative expression profiling reveals widespread coordinated evolution of gene expression across eukaryotes

Abstract

Comparative studies of gene expression across species have revealed many important insights, but have also been limited by the number of species represented. Here we develop an approach to identify orthologs between highly diverged transcriptome assemblies, and apply this to 657 RNA-seq gene expression profiles from 309 diverse unicellular eukaryotes. We analyzed the resulting data for coevolutionary patterns, and identify several hundred protein complexes and pathways whose expression levels have evolved in a coordinated fashion across the trillions of generations separating these species, including many gene sets with little or no within-species co-expression across environmental or genetic perturbations. We also detect examples of adaptive evolution, for example of tRNA ligase levels to match genome-wide codon usage. In sum, we find that comparative studies from extremely diverse organisms can reveal new insights into the evolution of gene expression, including coordinated evolution of some of the most conserved protein complexes in eukaryotes.

Fig. 1

Overview of phylogenetic expression profiling approach and data. a 18S-based cladogram of the species in this study. b Overview of our approach for inferring ortholog groups from RNA-seq transcripts (see also Methods and Supplementary Figure 1). c Heatmap of expression levels of the 4219 genes and 657 samples analyzed in this study. Samples and genes are clustered hierarchically. Color bar near the top shows the phylum of each sample, matching the colors in a

Fig. 2

Phylogenetic expression profiling reveals coordinated evolution within gene sets. a Overview of traditional phylogenetic profiling (PP; top) and phylogenetic expression profiling (PEP; bottom). PEP uses the quantitative gene expression levels across species rather than the binary presence/absence of a gene. Patterns of coordinated evolution hidden to PP can be potentially uncovered using PEP. b For ciliary genes, pairwise PP correlations increase with PEP correlation strength. c The 662 gene sets with significant PEP scores are clustered by the pairwise correlations between gene sets. The color bar below the dendrogram shows the 15 unique gene set groups the dendrogram was divided into and gene ontology enrichments for each group are highlighted in the same color. Black bars highlight notable groups of gene sets within the larger groups. d Proteasome genes were found to be undergoing coordinated evolution and are shown as a heatmap

Fig. 3

Comparison of PEP with cross-environment and cross-genotype expression profiling methods. a Scatterplot of the gene-set correlation medians for the gene sets found as significant using the PEP method introduced here vs. the correlation medians for the same gene sets in a data set of yeast expression across environments. The gene sets boxed in purple and green are highlighted in c and d, respectively. b Scatterplot of gene-set correlation medians as in a, but comparing PEP to yeast expression across deletion strains (i.e., across genotypes). c Comparison of gene-gene correlations for the deadenylation dependent mRNA decay gene set highlighted in a for the PEP method (lower triangle) and classical yeast coexpression correlations across environments (upper triangle). d Comparison of gene–gene correlations for the ribosome gene set highlighted in a as described in c

Fig. 4

Associations between expression and other sample information. a Scatterplot of the association between gene expression and the absolute value of latitude for the DNA polymerase POLH. Each point represents one RNA-seq sample. b Scatterplot of the association between expression of the aspartate-tRNA ligase DARS and the abundance of aspartate codons in the coding regions of each sample’s transcriptome