RNA-Seq offers a significant improvement in measuring cross-species gene expression over microarrays because it allows for orthologous genes from multiple species to be measured and compared to one another without the significant complications of species or strain-specific probe biases. For this reason, it is likely that RNA-Seq will replace microarrays in cross-species experiments. However, many of the challenges that were present in measuring gene expression using microarrays will continue to be present in RNA-Seq experiments because they originate from biological rather than technical sources. Overcoming these challenges in cross-species experiments will be particularly difficult because measurement artifacts can both mimic and obscure the effects of evolution. For example, we found that shorter genes tended to have higher expression divergence between species and are more likely to be differentially expressed. This finding unlikely to be an artifact of RNA-Seq measurements as it is consistent with a finding in a microarray-based study of drosophila [6
], and was observed in the microarray measurements of yeast [4
]. However, shorter genes that are truly DE are less likely to be detected as such than longer genes with expression changes of the same fold change because their transcripts produce fewer nucleotide fragments, and thus are measured relatively higher Poisson noise.
Normalizing read counts between samples can also be problematic. Under a neutral model of gene expression evolution, the number of genes with increased expression in one species versus another should equal the number with decreased expression. Deviations from symmetry can be viewed as evidence of selection [26
]. Most methods for normalizing RNA-Seq data assume that the overall expression of the two species will be the same, forcing some degree of symmetry. Because reference assemblies will usually differ in quality in cross-species experiments, it is important to use normalization strategies that are robust to the effects of missing regions and errors in genomes. Missing regions will have a limited effect on our normalization method because only genes that are present in both species are included in the normalization equations. The effect of errors can be better ameliorated by using an alignment strategy such as ours that allows for some mismatches between the reads and the reference genome. In our qPCR validation, our results did not show any systematic, directional differences compared to our RNA-Seq results. Despite this strategy, we did see a tendency for genes with consistent expression changes in S. bayanus
to have decreased expression. There was, however, no evidence that these genes represented a functional class of genes, and across the entire phylogeny genes with changes in expression were equally likely to have increased and decreased expression.
To confirm that our statistical approach adequately accounted for both technical and environmental variance in our experiment, we validated our DE calls using qPCR with RNA isolated from biological replicates grown separately from the samples used for RNA-Seq. By contrast, running qPCR on the same RNA used for RNA-Seq would have only assessed whether the technical variance was correctly assessed. RNA-Seq measurements are already known to have a high level of technical reproducibility, in our study and others [12
], and therefore validations of technical reproducibility add little value. Because our qPCR was run on independent samples, and our results generally reproduced within our calculated confidence, we are confident that significant p-values correspond to the probability of consistent results if the experiment were repeated.
The sensitivity of experiments to detect differentially expressed genes is determined by the ability to distinguish true changes in gene expression from alternate sources of variance. Sensitivity in RNA-Seq experiments will be limited for genes that are measured with low read counts and high Poisson counting noise. In this experiment, our samples were sequenced deeply enough that we could measure the expression of greater than 95% of genes with greater than 10 reads. This sequencing depth is much more difficult to achieve in mammals and other species with large transcriptomes. We found that for genes measured with fewer than 10 reads, the signals from biological and technical variance were largely obscured by Poisson sampling noise, and only large fold changes could be detected as significant. While deeper sequencing will reduce this Poisson noise, using the same number of reads to sequence an additional replicate will reduce both Poisson and non-Poisson variance, and should theoretically provide a greater increase in sensitivity than deeper sequencing.
In a review of studies of gene expression evolution, Clarke et al. [27
] noted that the reported DE rates between evolutionarily diverged organisms shows a high degree of variation between studies, and suggested that this variation was caused by differences in the sensitivity of the experiments. Additionally, the way in which biological variance is measured and accounted for in the experimental design is a crucial factor in determining not just how many genes will be detected as DE, but also the reproducibility of findings. For example, Cassone et al. [28
] reported that the calculated rate of DE genes between two insipient species of mosquitoes was 1-2% of genes when they applied a model that correctly controlled for the biological variance observed between colonies of mosquitoes, but would be measured at ~
54% if less stringent controls were applied. Their finding of little differential gene expression in these closely related strains is consistent with our results.
Large deviations from our linear model of expression divergence would be possible if a single mutation altered a gene's expression dramatically, and this event had a cascading effect on the expression of other genes. However, we believe that in wild populations such changes are likely to be selected against. For example, the change in expression of RME1 demonstrates that a single base mutation can dramatically alter the expression of a gene. However, despite the high divergence within the intergenic regions of these species, there were only ten examples of changes in expression that were of the same magnitude as the change observed in RME1. A role for purifying selection is further supported by our finding that the more diverged copy of duplicated genes also dropped in expression 35% of the time while there were no such cases for the less diverged copy. This finding is consistent with, and further clarifies, microarray-based findings of expression divergence following gene duplication events [29
While there are clear signs of purifying selection in our study, several lines of evidence support the premise that the primary driver of changes in expression is neutral drift. Only a small fraction of the genes differentially expressed between individual species showed evidence of being fixed in any lineage. While two Gene Ontology categories were over-represented in genes exhibiting lineage-specific DE, the majority of such genes were not members of these categories. Furthermore, the linear correlation between the expression divergence and the intergenic substitution rate suggests a clock-like divergence. While properties such as TATA boxes and gene length were associated with increased gene expression divergence, divergence rates within these categories still maintained a linear relationship with intergenic substitution rate, consistent with a drift model with a faster diffusion coefficient.
A recent study of a S. cerevisiae - S. bayanus
hybrid identified several categories of genes as being directionally differentially regulated between S. cerevisiae
and S. bayanus
, and suggested that these categories were under evolutionary selection [31
]. We did not replicate differential expression in any of their Gene Ontology categories in our study. There were, however, large differences between our study designs and the statistical analyses used. A hybrid construct is used to identify genes that are expressed at different levels due to differences in transcriptional efficiency caused by cis regulatory differences. In Tirosh et al. [32
] it was shown that genes observed as DE between species are not necessarily observed as DE within the hybrid, and vice versa. Additionally, our strains were prototrophic, whereas the hybrid strain had several auxotrophies and therefore required amino acid supplements in the media, making the growth conditions in the two experiments different. The auxotrophies may be particularly important differences in light of the fact that two of the five networks that Bullard et al. identified as directionally regulated in the hybrid (histidine and lysine biosynthesis) contained pathways that were not functioning normally because the strain was auxotrophic for his3 and lys2.