Downloaded from genome.cshlp.org on January 20, 2021 - Published by Cold Spring Harbor Laboratory Press Resource A comparison of gene expression and DNA methylation patterns across tissues and species Lauren E. Blake,1,9 Julien Roux,1,2,3,9 Irene Hernando-Herraez,4 Nicholas E. Banovich,1 Raquel Garcia Perez,4 Chiaowen Joyce Hsiao,1 Ittai Eres,1 Claudia Cuevas,1 Tomas Marques-Bonet,4,5,6,7 and Yoav Gilad1,8 1 Department of Human Genetics, University of Chicago, Chicago, Illinois 60637, USA; 2Department of Biomedicine, University of Basel, 4031 Basel, Switzerland; 3Swiss Institute of Bioinformatics, 4031 Basel, Switzerland; 4Universitat Pompeu Fabra, Institute of Evolutionary Biology, 88 08003 Barcelona, Spain; 5Passeig de Lluís Companys, Catalan Institution of Research and Advanced Studies, 23 08010 Barcelona, Spain; 6Barcelona Institute of Science and Technology, Centre for Genomic Regulation, 88 08003 Barcelona, Spain; 7Universitat Autònoma de Barcelona, Institut Català de Paleontologia Miquel Crusafont, 08193 Barcelona, Spain; 8 Department of Medicine, University of Chicago, Chicago, Illinois 60637, USA Previously published comparative functional genomic data sets from primates using frozen tissue samples, including many data sets from our own group, were often collected and analyzed using nonoptimal study designs and analysis approaches. In addition, when samples from multiple tissues were studied in a comparative framework, individuals and tissues were con- founded. We designed a multitissue comparative study of gene expression and DNA methylation in primates that minimizes confounding effects by using a balanced design with respect to species, tissues, and individuals. We also developed a com- parative analysis pipeline that minimizes biases attributable to sequence divergence. Thus, we present the most comprehen- sive catalog of similarities and differences in gene expression and DNA methylation levels between livers, kidneys, hearts, and lungs, in humans, chimpanzees, and rhesus macaques. We estimate that overall, interspecies and inter-tissue differences in gene expression levels can only modestly be accounted for by corresponding differences in promoter DNA methylation. However, the expression pattern of genes with conserved inter-tissue expression differences can be explained by corre- sponding interspecies methylation changes more often. Finally, we show that genes whose tissue-specific regulatory patterns are consistent with the action of natural selection are highly connected in both gene regulatory and protein–protein inter- action networks. [Supplemental material is available for this article.] Gene regulatory differences between humans and other primates Khan et al. 2013) and ancestral or derived phenotypes that may are hypothesized to underlie human-specific traits (King and be relevant to human diseases (Cooper and Shendure 2011; Wilson 1975). Over the last decade, dozens of comparative geno- Romero et al. 2012). mic studies focused on characterizing mRNA expression level dif- One caveat that is shared among practically all comparative ferences between primates in a large number of tissues (e.g., studies in primates is related to difficulty in obtaining multiple tis- Enard et al. 2002; Khaitovich et al. 2005; Blekhman et al. 2008; sue samples from the same individual. To date, there have been no Brawand et al. 2011; Barbash and Sakmar 2017), typically focusing published comparative studies in primates that have analyzed on differences between humans and other primates. A few studies multiple tissues sampled from the same individuals across multi- have also characterized inter-primate differences in regulatory ple species in a balanced design (Romero et al. 2012). As a result, mechanisms and phenotypes other than gene expression levels, regulatory differences between tissues are always confounded such as DNA methylation levels, chromatin modifications and ac- with regulatory differences between individuals (Blekhman et al. cessibility, and protein expression levels (Cain et al. 2011; Pai et al. 2008; Brawand et al. 2011; Pai et al. 2011; Chen et al. 2019). 2011; Hernando-Herraez et al. 2013, 2015a,b; Ward et al. 2013; In turn, catalogs from these studies cannot be used to compare Stergachis et al. 2014; Zhou et al. 2014; Villar et al. 2015). These tissue-specific regulatory differences between species to inter- studies often construct catalogs of gene expression levels and other tissue differences in regulatory variation within species (see mechanisms. These catalogs have been useful to better understand “Discussion” in Pai et al. 2011). the evolutionary processes that led to adaptations in humans Our group and others often use previously published catalogs (Enard et al. 2002; Cáceres et al. 2003; Karaman et al. 2003; of comparative data from primates in the different studies. Khaitovich et al. 2004a, 2005; Gilad et al. 2005; Lemos et al. Although we do not expect previously observed patterns to 2005; Blekhman et al. 2008, 2010; Nowick et al. 2009; Babbitt be erroneous, we are aware that data on gene-specific interspecies et al. 2010; Pai et al. 2011; Shibata et al. 2012; Capra et al. 2013; regulatory differences, and especially data that pertain to com- parisons of divergence across tissues, may be inaccurate for the reasons we discussed above. We thus designed the current 9 These authors contributed equally to this work. Corresponding author: gilad@uchicago.edu Article published online before print. Article, supplemental material, and publi- © 2020 Blake et al. This article, published in Genome Research, is available cation date are at http://www.genome.org/cgi/doi/10.1101/gr.254904.119. under a Creative Commons License (Attribution-NonCommercial 4.0 Interna- Freely available online through the Genome Research Open Access option. tional), as described at http://creativecommons.org/licenses/by-nc/4.0/. 250 Genome Research 30:250–262 Published by Cold Spring Harbor Laboratory Press; ISSN 1088-9051/20; www.genome.org www.genome.org Downloaded from genome.cshlp.org on January 20, 2021 - Published by Cold Spring Harbor Laboratory Press Gene regulatory patterns across primate tissues A B liver, and lung tissue samples from four human, four chimpanzee, and four rhe- sus macaque individuals (Fig. 1A; Sup- plemental Table S1A). From these 48 samples, we harvested RNA and DNA in parallel (Methods). After confirming that the RNA from all samples was of ac- ceptable quality (Supplemental Fig. S1A; Supplemental Table S1B), we performed RNA-sequencing to obtain estimates of gene expression levels. Additional details C D about the donors, tissue samples, sample processing, and sequencing information can be found in “Methods” and Supple- mental Table S1. We estimated gene expression levels using an approach designed to prevent biases driven by sequence divergence across species, similar to the approach of Blake et al. (2018b). Briefly, we map- ped RNA-sequencing reads to each species’ respective genome. To compare E F gene expression levels across species, we only calculated the number of reads mapping to exons that can be classified as clear orthologs across all three species (Supplemental Table S1B). We excluded data from genes that were lowly ex- pressed in more than half of the samples as well as data from one human heart sample that was an obvious outlier, prob- ably owing to a sample swap (Supple- Figure 1. Surveying gene expression and DNA methylation in diverse tissues across primates. (A) Study mental Fig. S2A,B). We normalized the design. (B) Principal components analysis (PCA) of gene expression levels in 47 samples. (C) Normalized distribution of gene expression levels to gene expression (quantile-transformed RPKMs) from four donors in the GTEx heart collection (“heart remove systematic expression differences same individuals”) compared to the lungs from different (“lung different individuals”) and the same do- between species (maximizing the num- nors’ lungs (“lung same individuals”) in AC020922.1. (D) PCA of average methylation levels 250 bp up- ber of genes with invariant expression stream and downstream in 47 samples. (E) Density function of the correlation between gene expression and DNA methylation levels in human–chimpanzee orthologous genes. (F) Density function of the cor- levels across species corresponds to our relation between gene expression and DNA methylation levels in genes orthologous across humans and null hypothesis) (Methods). Through rhesus macaques. this process, we obtained TMM- and cyclic loess-normalized log2 counts per million (CPM) values for 12,184 ortholo- study to produce a new comprehensive catalog of comparative gous genes to be used in downstream analyses (Supplemental Ta- gene expression and DNA methylation data from humans, chim- ble S2). panzees, and rhesus macaques, attempting to minimize possible Elements of study design, including sample processing, have confounders. been shown to impact gene expression data (Gilad and Mizrahi- The goal of our study is not to challenge previous conclusions Man 2015). Consequently, we tested the relationship between a or document specific differences between the current and previous large number of technical factors recorded throughout our exper- data. Instead, we aim to provide a new and more accurate com- iments and the biological variables of interest in our study, namely parative catalog of inter-tissue and interspecies differences in tissue and species (Methods; Supplemental Material; Supplemen- gene regulation between humans and other primates, with sub- tal Table S3A,B). We found that there were no technical confound- stantial sample and study design documentation. Overall, we be- ers with tissue but two technical factors were confounded with lieve that this catalog can be useful for many future applications species: time postmortem until collection and RNA extraction and can serve as a new benchmark for regulatory divergence in date (Supplemental Fig. 1B,C). Because of the opportunistic nature primates. of sample collection, these confounders are practically impossible to avoid in comparative studies of primates (especially apes). We discuss possible implications of these confounders throughout Results the paper. Study design and data collection Gene expression varies more across tissues than across species To comparatively study gene expression levels and DNA methyla- We first examined broad patterns in the gene expression data. A tion patterns in primates, we collected primary heart, kidney, principal component analysis (PCA) and a separate clustering Genome Research 251 www.genome.org Downloaded from genome.cshlp.org on January 20, 2021 - Published by Cold Spring Harbor Laboratory Press Blake et al. Table 1. Pairwise DE genes at FDR 1% DE between tissues (within species) Heart-kidney Heart-liver Heart-lung Kidney-liver Kidney-lung Liver-lung Human 4224 4776 4037 4248 3695 4701 Chimpanzee 4971 6295 5365 4625 4623 6247 Rhesus macaque 5980 6933 6814 5126 5721 7027 DE between species Heart Kidney Liver Lung Human versus chimpanzee 2195 799 1364 805 Human versus rhesus macaque 4098 2347 2868 2833 Chimpanzee versus rhesus macaque 2781 2286 3139 1917 analysis indicated that, as expected (Brawand et al. 2011; Barbosa- trasted liver and kidney, and the largest number of DE genes be- Morais et al. 2012; Merkin et al. 2012), the primary sources of tween liver and either heart or lung (Table 1; Supplemental Table gene expression variation are tissue (regression of PC1 by tissue = S4B). However, because our data were produced from bulk RNA-se- 0.81; P < 10−14; regression of PC2 by tissue = 0.70; P < 10−10) (Fig. quencing, we were unable to determine the impact of cell composi- 1B; Supplemental Tables S1A,B, S3A,B; Supplemental Fig. S3), fol- tion on the number of inter-tissue DE genes. lowed by species (regression of PC2 by species = 0.27; P < 10−3) We used the same framework of linear modeling to identify (Supplemental Tables S1A,B, S3A,B). We then confirmed that, gene expression differences between species, within each tissue globally, gene expression levels across tissues from the same indi- (Supplemental Table S4A). Depending on the tissue and species vidual are more highly correlated than gene expression levels we considered, we identified between 799 and 4098 interspecies across tissues from different individuals (Supplemental Fig. S2C). DE genes (at FDR = 1%) (Table 1). As expected, given the known This observation supports the intuitive notion that collecting phylogeny of the three species, within each tissue we classified and analyzing multiple tissues from the same individual is highly far fewer DE genes between humans and chimpanzees than be- desirable in functional genomics studies. tween either of these species and rhesus macaques (Supplemental We sought further explicit evidence that incorporating ba- Table S4B). lanced collection of multiple tissues from the same individuals is It is a common notion that genes with tissue-specific expres- an effective study design. To do so, we used lung and heart data sion patterns may underlie tissue-specific functions. Previous cat- from the GTEx Consortium (The GTEx Consortium 2017). We first alogs of such patterns in primates were always confounded by identified differentially expressed (DE) genes between lung and the effect of individual variation (because each tissue was sampled heart; we designated these classifications, which are based on hun- from a different individual). To classify tissue-specific genes using dreds of samples, as the “truth” (Supplemental Material; our data, we focused on genes that are either up-regulated or down- Supplemental Table S3E). Next, we repeatedly identified DE genes regulated in a single tissue relative to the other three tissues (within between lung and heart using GTEx data from randomly chosen one or more species). We define such genes as having a “tissue-spe- sets of just four samples from each tissue. We then compared the cific” expression pattern, acknowledging that this definition may results to DE genes identified from an equivalent analysis of sets only be relevant in the context of the four tissues we considered of four samples from each tissue, in which the tissue samples orig- here. inated from the same donor. Compared to the designated “truth,” Using this approach and considering the human data across DE analyses using data from tissue samples that are matched for all tissue comparisons, we identified 5284 genes with tissue-specif- donors result in a higher ratio of true positives to false positives ic gene expression patterns (FDR 1%) (Fig. 2A–D). By performing than analyses using tissue samples that are unmatched for donors similar analyses using the chimpanzee and rhesus macaque data, (P = 0.03) (Supplemental Table S3F). Given the small number of we found that the degree of conservation of tissue-specific expres- false positives in both data sets, study design is unlikely to impact sion patterns is higher than expected by chance (P < 10−16) (Fig. large-scale, highly robust trends across species. However, this study 2A–D). This observation is robust with respect to the statistical cut- design choice is particularly important if one is interested in indi- offs we used (Supplemental Table S6), indicating that many of vidual genes (as shown by an example in Fig. 1C). these conserved tissue-specific regulatory patterns are likely of functional significance. To broadly analyze the biological function of genes with con- Putatively functional tissue-specific gene expression patterns served tissue-specific expression, we performed a Gene Ontology To analyze the pairwise regulatory differences across tissues and (GO) enrichment analysis (Supplemental Material). We found species, we used the framework of a linear model (Methods). We these genes are indeed highly enriched with functional annota- first identified (at FDR <1%) 3695–7027 (depending on the compar- tions that are relevant to the corresponding tissue (Supplemental ison we considered) differences in gene expression levels between Tables S7A–D, S8). For example, genes with conserved heart-specif- tissues, within each species (Table 1; Supplemental Table S4A–C). ic expression patterns were enriched in GO categories related to Overall, the patterns of inter-tissue differences in gene expression muscle filament sliding (e.g., ACTA1, MYL2) and cardiac muscle levels are similar in the three species, significantly more so than ex- contraction (e.g., MYBPC3, TNNI3). pected by chance alone (P < 10−16, hypergeometric distribution) (Supplemental Material; Supplemental Table S5). A range of 17%– 26% of inter-tissue DE genes have conserved inter-tissue expression Functional analysis of gene regulatory differences patterns in all three species (Supplemental Table S5). Regardless of We sought further evidence that the classification of genes with species, we found the fewest inter-tissue DE genes when we con- conserved tissue-specific expression patterns is meaningful. To 252 Genome Research www.genome.org Downloaded from genome.cshlp.org on January 20, 2021 - Published by Cold Spring Harbor Laboratory Press Gene regulatory patterns across primate tissues A B sion patterns or genes whose tissue- specific expression patterns are not conserved (P < 10−5). When we only con- sidered genes that appear as nodes in the network, we found that genes with con- served tissue-specific expression patterns are more likely to be classified as hubs in the networks than genes without tissue- specific expression patterns or genes whose tissue-specific expression patterns are not conserved (P < 0.007). C D Motivated by these findings, we focused on gene expression patterns that are consistent with the action of nat- ural selection (as described in Blekhman et al. 2008; Supplemental Material; Sup- plemental Table S7E). We found that genes whose expression patterns are con- sistent with the action of either stabiliz- ing or directional selection (top 10%) (Supplemental Table S7F) have more interactions with other genes in the network than genes whose expression E patterns are not consistent with the action of natural selection (bottom 10%; P < 0.05 for all comparisons) (Fig. 2E). This observation is fairly robust with respect to percentile cutoff (Supple- mental Table 7F). We repeated a similar analysis using protein–protein interaction data from the Human Protein Atlas (Uhlén et al. 2015; Yu et al. 2015; Lindskog 2016; Thul and Lindskog 2018) in all four tis- sues. We again found that genes whose expression patterns are consistent with F selection have more annotated protein– protein interactions (P < 0.05 in all eight comparisons) (Fig. 2F; Supplemental Table 7G). These interaction results sug- gest that functionally important genes are carefully regulated. Furthermore, this tight regulation occurs at both the gene expression and protein levels in primates. Variation in DNA methylation across tissues and species We used low-coverage whole-genome bisulfite sequencing (BS-seq) to study Figure 2. Tissue-specific DE genes (FDR = 0.01). The number of conserved tissue-specific DE genes across all three species is greater than the number expected by chance: (A) heart; (B) kidney; (C) liver; DNA methylation patterns in each sam- (D) lung. In each tissue, genes with tissue-specific regulatory patterns that are consistent with the action ple. The bisulfite conversion reaction ef- of natural selection (top) are more likely to appear in gene coexpression networks (E) and have an in- ficiency was higher than 99.4% for all creased number of protein–protein interactions (F ) than those that are less consistent with the action samples (Supplemental Table S1C). of natural selection (bottom). (∗ ) P < 0.05; (∗∗∗ ) P < 0.001. The x-axes of E and F are cut off at 80 interac- tions for readability. In both cases, <5% of the data points are beyond this cutoff. Following sequencing, we mapped the high-quality BS-seq reads to in silico bisulfite-converted genomes of the corre- do so, we considered transcription coexpression networks (Stuart sponding species. We measured DNA methylation levels in 12.5 et al. 2003; Zhang and Horvath 2005) based on GTEx data from million to 22.9 million CpG sites per sample, with a minimum cov- heart and lung (Pierson et al. 2015). We found that genes with con- erage of two sequencing reads per site (Supplemental Table S1C). served tissue-specific expression patterns are more likely to appear We estimated local methylation levels by smoothing the data as nodes in the networks than genes without tissue-specific expres- across nearby CpG sites (Supplemental Material; Supplemental Genome Research 253 www.genome.org Downloaded from genome.cshlp.org on January 20, 2021 - Published by Cold Spring Harbor Laboratory Press Blake et al. Table 2. Pairwise differentially methylated regions (DMRs) in autosomal chromosomes Tissue DMRs (within species) Heart-kidney Heart-liver Heart-lung Kidney-liver Kidney-lung Liver-lung Human 30,291 (1.7%) 22,561 (1.3%) 14,208 (0.8%) 17,910 (1.1%) 16,521 (1.0%) 12,842 (0.8%) Chimpanzee 17,699 (1.1%) 28,767 (1.7%) 8524 (0.5%) 23,847 (1.4%) 12,076 (0.7%) 22,107 (1.3%) Rhesus macaque 23,023 (1.5%) 41,280 (2.7%) 7026 (0.5%) 35,910 (2.4%) 15,889 (1.1%) 32,636 (2.1%) Species DMRs Heart Kidney Liver Lung Human versus chimpanzee 14,504 (1.1%) 10,667 (0.8%) 9476 (0.8%) 8617 (0.7%) Human versus rhesus macaque 25,539 (5.3%) 21,292 (4.5%) 21,639 (4.5%) 17,696 (3.9%) Cutoff recommended by Hansen et al. (2012) (percentage out of total pairwise methylated regions). Figs. S4–S6; Hansen et al. 2012). To facilitate a comparison of meth- species tend to be conserved, consistent with the observations of ylation levels across species, we annotated 10.5 million ortholo- previous studies (Martin et al. 2011; Molaro et al. 2011; Pai et al. gous CpGs in the human and chimpanzee genomes, as well as a 2011; Hernando-Herraez et al. 2013, 2015b). smaller set of 2.4 million orthologous CpGs in all three primate ge- We next focused specifically on tissue-specific DMRs, as these nomes (Supplemental Table S1C–E). To identify differences in DNA may contribute to tissue-specific function. In contrast to differenc- methylation levels between tissues and species we again used a lin- es in DNA methylation between any pair of tissues, a tissue-specific ear model framework (Methods; Fig. 1D). Focusing on DNA meth- DMR is defined as having a similar methylation level in three of the ylation patterns across tissues within species, we identified tissues we considered, but a significantly different DNA methyla- between 7026 and 41,280 differentially methylated regions be- tion level in the remaining tissue. We found that there were tween tissues, within species (T-DMRs) (Table 2; Supplemental Ta- more DMRs specific to liver (3278 to 11,433 DMRs depending on ble S9A; Blake et al. 2018a). Pairwise comparisons between hearts the species) than to kidney (2300 to 3957 DMRs), heart (1597 to and lungs showed the lowest number of T-DMRs, regardless of spe- 2969 DMRs), or lung (453 to 5018 DMRs) (Fig. 3A–D; Supplemen- cies (7026 in rhesus macaques, 8524 in chimpanzees, 14,208 in hu- tal Table S10B). Tissue-specific DMRs are highly conserved regard- mans), whereas comparisons involving heart and liver showed the less of the comparisons we made (P < 10−13 for all comparisons, at largest number of T-DMRs (22,561 in humans, 28,767 in chimpan- least 25% bp overlap was required to be considered shared). zee, and 41,280 in rhesus macaques) (Ta- ble 2). We found that human T-DMRs overlapped genic and regulatory features A B significantly more than expected by chance. In particular, there is an enrich- ment of T-DMRs in intergenic regions, in- trons, 5′ UTRs, 3′ UTRs, and active enhancers (P < 0.04 for all tests) (as de- fined by Andersson et al. 2014; Supple- mental Table S9B). We found strong evidence for C D T-DMR conservation across all three spe- cies (P < 10−16 across all comparisons) (Supplemental Table S10A). Although this level of conservation is higher than expected by chance, we recognize that in each tissue comparison we performed, we had incomplete power to identify T-DMRs and so the true conservation of T-DMRs is expected to be even higher. E F To compare T-DMRs across species more effectively, we considered DNA methyla- tion data from all T-DMR orthologous re- gions that were classified as such in at least one species. When we performed hierarchical clustering using ortholog- ous DNA methylation data from these T-DMRs, the data clustered first by tissue, then by species (Supplemental Fig. S7). This trend is robust with respect to the species used to initially locate T-DMRs (Supplemental Figs. S8, S9). Thus, our re- Figure 3. Tissue-specific DMRs (FDR = 0.01). The number of conserved tissue-specific DMRs is greater sults suggest that in general, inter-tissue than expected by chance: (A) heart; (B) kidney; (C ) liver; (D) lung. Genes with the closest TSSs to con- DNA methylation differences within a served tissue-specific DMRs are enriched for relevant functional annotations in hearts (E) and livers (F). 254 Genome Research www.genome.org Downloaded from genome.cshlp.org on January 20, 2021 - Published by Cold Spring Harbor Laboratory Press Gene regulatory patterns across primate tissues In all four tissues, >59% of conserved DMRs are hypometh- levels, P < 0.003 in all tissues, t-test comparing the absolute values ylated in a tissue-specific manner. We evaluated the overlap be- of the effect sizes for both groups of DE genes). tween tissue-specific DMRs and genomic regions marked with Next, we examined the genes in which DNA methylation dif- H3K27ac, a mark often associated with active gene expression ferences may underlie inter-tissue gene expression differences (ex- (The ENCODE Project Consortium 2012). We found that con- ample in Fig. 4A–C). Using adaptive shrinkage, we found that 9%– served hypomethylated tissue-specific DMRs were annotated 17% of inter-tissue gene expression differences could potentially with H3K27ac more frequently than tissue-specific DMRs identi- be explained by DNA methylation differences across tissues (FSR fied only in humans (P < 0.001, difference of proportions test) 5%). When we performed the control analysis and considered (Supplemental Material; Supplemental Table S10C). We then only data from genes that were not DE between tissues, <1% of ef- asked about the potential impact of these conserved hypometh- fect sizes differed once we accounted for the DNA methylation ylated tissue-specific DMRs on the expression of nearby genes. We data (Fig. 4F; Supplemental Table S11C–E). found that genes with the closest TSSs to conserved tissue-specific Finally, we focused on regulatory patterns that are most likely DMRs are highly enriched with relevant functional annotations in to be functional—namely, conserved inter-tissue gene regulatory hearts and livers (the tissues with the largest numbers of conserved differences. These differences were more likely to be explained hypomethylated tissue-specific DMRs) (Fig. 3E,F; Supplemental by variation in DNA methylation levels than nonconserved in- Table S10D; Supek et al. 2011). For example, conserved heart-spe- ter-tissue gene expression differences (minimum difference is cific DMRs are closest to genes in cardiovascular-related pathways, 7%, P < 0.005 for all comparisons; at FDR <5% and FSR <5%) including ventricular cardiac muscle cell development, canonical (Fig. 4D,E,H; Supplemental Table S11C–E). This observation is ro- Wnt signaling pathway, and the MAPK7 cascade. Overall, these bust with respect to the FDR and FSR cutoff used (Supplemental observations suggest that conserved tissue-specific DMRs are likely Table S11C–E). Indeed, the correlation between DNA methylation to underlie tissue-specific gene regulation in primates. and gene expression data is higher for genes with conserved inter- tissue expression patterns compared to genes whose expression patterns are not conserved (Fig. 1E,F). Interspecies differences in gene expression and DNA methylation One way to maintain conserved inter-tissue gene expression levels differences could be through DNA methylation level differences. Our comparative catalog can be used to identify DNA methylation We compared the genes whose variation in inter-tissue gene ex- differences that could potentially explain gene expression differ- pression can potentially be explained by variation in DNA methyl- ences across species and tissues. To do so, we first identified ation levels (assuming no independent effect of an unobserved 7,725 orthologous genes with expression data and corresponding factor) to all genes with conserved inter-tissue expression differ- promoter DNA methylation data in humans and chimpanzees, ences. We found that these genes are enriched for “essential tissue and 4155 orthologous genes with the same information for all functions” (Supplemental Table S11F). For example, the heart three species. We then determined to what extent divergence in genes are enriched for cardiac and smooth muscle contraction, DNA methylation levels could potentially underlie interspecies whereas those in liver are enriched for regulation of cholesterol differences in gene expression by comparing the gene expression transport and hormone secretion (Fig. 4G; Supplemental Table effect size associated with “species” before and after accounting S11F). These observations suggest that DNA methylation levels for methylation levels. To determine significant effect size differ- may mark or even drive differences in the expression levels in func- ences, we applied adaptive shrinkage (Stephens 2017)—a flexible tionally relevant genes. empirical Bayes approach for estimating false discovery rate (Methods). We note that this mediation approach does not consid- er the possibility that a third, unobserved event may be causally re- Discussion sponsible for both the DNA methylation and expression patterns. We designed a comparative study of gene regulation in humans, Considering DE genes between humans and chimpanzees (in chimpanzees, and rhesus macaques that minimized confounding at least one tissue), we found that between 11% and 25% of genes effects and bias. Consistent with previous studies, we found a high (depending on tissue) showed a difference in the effect of species degree of conservation in gene expression levels when we consid- on gene expression levels once average promoter methylation lev- ered the same tissue across species (Barlow 1993; Brawand et al. els were accounted for (significant difference in effect size classi- 2011; Sharp et al. 2011; Lin et al. 2014; Gallego Romero et al. fied at FSR 5%, represented by red in Supplemental Fig. S10; 2015). We also found evidence for conservation of tissue-specific Supplemental Table S11A). As a control analysis, we considered DMRs. Our observations are qualitatively consistent with those only the genes that were not originally classified as DE between of previous studies that mostly used microarrays to measure DNA humans and chimpanzees, and we found that the difference in methylation levels (Pai et al. 2011; Hernando-Herraez et al. the effect size of species on gene expression levels was reduced in 2015b; Tsankov et al. 2015), however, the high resolution of our <1% of genes once DNA methylation was accounted for (FSR BS-seq data allowed us to examine a much larger number of CpG 5%) (Supplemental Fig. S10; Supplemental Table S11A); thus, our sites. Thus, we were able to show that although DNA methylation approach is well calibrated. can potentially explain a modest proportion of expression differ- We applied the same approach to the human and rhesus ma- ences between tissues (Pai et al. 2011), it is more likely to impact caque data and found that the percentage of genes for which gene conserved tissue-specific gene expression levels. expression differences could potentially be explained by DNA We created and made available the most comprehensive, and methylation differences ranges from 21% in the lung to 40% in likely most accurate comparative catalog of gene expression and the liver (Supplemental Fig. S11; Supplemental Table S11B). This methylation levels in humans, chimpanzees, and rhesus ma- observation may reflect the more extreme gene expression differ- caques. Comparative functional genomic studies in primates, in- ences between humans and rhesus macaques than between hu- cluding from our own laboratory, often are not designed to test mans and chimpanzees (before accounting for DNA methylation specific hypotheses. Rather, many of these comparative genome- Genome Research 255 www.genome.org Downloaded from genome.cshlp.org on January 20, 2021 - Published by Cold Spring Harbor Laboratory Press Blake et al. A B C D E F G H Figure 4. Inter-tissue DNA methylation and gene expression levels (FDR = 0.05 and FSR = 0.05). (A–C) A representative example of the PRKACA gene in which the variation of methylation levels (A) may explain the differences in gene expression levels (B) between human heart and kidney. (C) The residuals of normalized gene expression levels after regressing out methylation levels. (D–G) Next, we compared the tissue effect sizes before and after controlling for DNA methylation levels in inter-tissue DE and non-DE genes, separately. Genes in red are significant at S-value <0.05. Effect size differences in conserved DE genes in human heart relative to human kidney (D), nonconserved DE genes in human heart relative to human kidney (E), and non-DE genes in human heart and human kidney (F ). (G) The conserved DE genes are enriched for heart-related function. (H) Variation in DNA methylation is more likely to explain variation in conserved DE genes than nonconserved DE genes (DE in the human tissues listed, but not in the same tissues in chimpanzees). scale studies aim to build catalogs of similarities and differences in help us better understand the evolutionary processes that led to gene regulation between humans and other primates. These cata- adaptations in humans. These catalogs are also used to establish in- logs have been shown to be quite useful; for example, they can be formed models of the relative importance of changes in different used to identify interspecies regulatory changes that have likely molecular mechanisms to regulatory evolution (Khaitovich et al. evolved under natural selection (Enard et al. 2002; Cáceres et al. 2004b; Warnefors and Eyre-Walker 2012) and to inform us about 2003; Karaman et al. 2003; Khaitovich et al. 2004a, 2005; Gilad ancestral or derived phenotypes that may be relevant to human et al. 2005; Lemos et al. 2005; Blekhman et al. 2008, 2010; diseases (Cooper and Shendure 2011; Gallego Romero et al. Nowick et al. 2009; Babbitt et al. 2010; Pai et al. 2011; Shibata 2015). Ultimately, comparative catalogs of gene regulatory pheno- et al. 2012; Capra et al. 2013; Khan et al. 2013), and thereby types are used to develop and test specific hypotheses regarding 256 Genome Research www.genome.org Downloaded from genome.cshlp.org on January 20, 2021 - Published by Cold Spring Harbor Laboratory Press Gene regulatory patterns across primate tissues the connection between interspecies regulatory changes and phys- Methods iological, anatomical, and cognitive phenotypic differences be- tween species. Sample description In this study, we used a comparative catalog to identify spe- We collected heart (left ventricle), kidney (cortex), liver, and lung cies-specific and, in particular, tissue-specific regulatory patterns, tissues from four individual donors in human (Homo sapiens, all of because these genes are often drug targets (Dezső et al. 2008) and reported Caucasian ethnicity), chimpanzee (Pan troglodytes), and are likely important for the evolution of human traits (Blekhman Indian rhesus macaque (Macaca mulatta), for a total of 48 samples et al. 2008). We showed that genes with conserved tissue-specific (3 species × 4 tissues × 4 individuals) (Fig. 1A). The choice of these regulatory patterns have more regulatory interactions and pro- tissues was guided by their relative homogeneity with respect to tein–protein interactions than genes whose regulatory patterns cellular composition (e.g., Balashova and Abdulkadyrov 1984), are not conserved or are not tissue-specific. These patterns became which does not change substantially across primate species. In even more pronounced when we focused on genes whose expres- contrast, other tissues, such as brain subparts, differ substantially sion patterns are consistent with the action of natural selection. in cellular composition across primates (Brodal 1983), which could Together, these observations consistently support the inference potentially confound the analyses. that when genes perform an important function that needs to be Human samples were obtained from the National Disease carefully regulated, evolution can act on multiple levels of the reg- Research Interchange (IRB protocol 14378B). Non-human samples were obtained from several sources, including the Yerkes primate ulatory cascade in primates. center and the Southwest Foundation for Biomedical Research, un- Focusing on species-specific patterns of tissue-specific gene der IACUC protocol 71619. When possible, samples were collected regulation, our observations can help formulate specific functional from adult individuals whose cause of death was unrelated to the hypotheses regarding human-specific adaptations. For example, tissues studied. genes with tissue-specific gene regulation identified only in hu- mans are enriched for GO pathways that may contribute to hu- man-specific features, including sodium ion import across the RNA library preparation and sequencing plasma membrane in kidneys (e.g., SLC9A3 and TRPM4), the gly- In total, we prepared 48 unstranded RNA-sequencing libraries as cogen biosynthetic process in livers (e.g., PGM1 and AKT1), and previously described (Marioni et al. 2008; Blekhman et al. 2010). paraxial mesoderm morphogenesis in lungs (e.g., MST1R and Twenty-four barcoded adapters were used to multiplex different MAP9). samples on two pools of libraries. RNA-sequencing libraries were sequenced on 26 lanes on four different flow-cells on an Illumina HiSeq 2500 sequencer in either the Gilad laboratory or Consideration of study design and record keeping at the University of Chicago Genomics Facility (50-bp single-end Regardless of the model system used and the types of data that are reads) (Supplemental Material; Supplemental Table S1B). collected, study design is critical. Perhaps because comparative studies in primates typically rely on opportunistic sample collec- Quantifying the number of RNA-seq reads from orthologous tion, there are not recognized standards for study design that are genes kept and consistently reported in most existing studies (including We used FastQC (version 0.10.0; https://www.bioinformatics many earlier studies from our group). We thus believe that it is .babraham.ac.uk/projects/fastqc/) to generate a read quality re- worthwhile to explicitly discuss a few important considerations re- port and TrimGalore (version 0.2.8; https://www.bioinformatics garding study design and the recording of metadata. .babraham.ac.uk/projects/trim_galore/), a wrapper based on Without a balanced study design, it would have been impos- Cutadapt (version 1.2.1) (Martin 2011), to trim adapter sequences sible to independently estimate the effects of individual, tissue, from RNA-seq reads. We trimmed using a stringency of 3. To cut and species on our data. Because the sources of confounding fac- the low-quality ends of reads, we used a quality threshold (Phred tors are difficult to predict in advance, we strongly recommend score) of 20. Reads shorter than 20 nt after trimming were eliminat- that samples are collected using a balanced design with respect ed before mapping (Supplemental Table S1B). to as many parameters as possible. These include the distribution For each sample, we used TopHat2 (version 2.0.8b) (Kim et al. of tissue samples per individual, the number of individuals from 2013) to map the reads to the correct species’ genome: human reads each species, sex, age range, cause of death and collection time to the hg19 genome, chimpanzee reads to the panTro3 genome, (in the case of postmortem tissues), or sample collection and cell and rhesus macaque reads to the rheMac2 genome (Supplemental culturing (in the case of iPSC-based models). All steps of sample Material). Expression level estimates may be biased across the processing (RNA extraction, library preparation, and so forth) species owing to factors such as mRNA transcript size and different genome annotation qualities. To circumvent these issues, we should be done in batches that are randomized or balanced with only retained reads that mapped to a set of 30,030 Ensembl respect to species, tissue, and any other variables of interest. gene orthologous meta-exons available for each of the three ge- Most important, all sample processing steps should be record- nomes, as described and used previously (Blekhman et al. 2010; ed in a sample history file that includes anything that happened to Blekhman 2012; Gallego Romero et al. 2015). We defined the num- any sample. We have documented many of these steps in ber of reads mapped to orthologous genes as the sum of the reads Supplemental Table S1A–E. This documentation can help provide mapped to the orthologous meta-exons of each gene. We quanti- evidence that a phenomenon is driven by biological rather than fied gene expression levels using the program coverageBed from technical factors. It may also benefit future studies by facilitating the BEDTools suite (Quinlan and Hall 2010) and then performed effective meta-analysis of multiple data sets, which would help TMM and cyclic loess normalization (Supplemental Material). to address the problems of tissue availability and small sample siz- Because all three releases of GTEx use hg19 and only a fraction es. We believe that, moving forward, it should be a requirement of the 1000 Genomes Project data are available in GRCh38 coordi- that these metadata are available with every published compara- nates, we also used hg19. However, to show that the results we re- tive genomic data set. port would not change much if we used the GRCh38 build, we Genome Research 257 www.genome.org Downloaded from genome.cshlp.org on January 20, 2021 - Published by Cold Spring Harbor Laboratory Press Blake et al. leveraged the fact that differential expression analysis compares limma + voom framework to identify significant differences in gene expression levels from groups of samples (e.g., human liver one tissue versus the group of the other tissues (Law et al. 2014; samples to human lung samples). Therefore, we compared the Ritchie et al. 2015). ranks of the normalized gene expression levels in the 15 human To identify interspecies differences in gene expression pat- samples mapped using hg19 to the same samples mapped to terns across tissues within species (tissue-by-species interactions), GRCh38. The correlations of these ranks were extremely high (me- we used the limma + voom framework and looked for the signifi- dian Pearson’s correlation = 0.96). These strong correlations sug- cance of tissue-by-group interactions. In one analysis, the groups gest that our general conclusions—and indeed, many genes we were great ape versus rhesus macaque and in another analysis, identified as DE—would remain if we had used the GRCh38 build. the groups were human versus non-human primates. To minimize the number of interactions, we compared one tissue relative to a group of the other three tissues (e.g., great ape vs. rhesus macaque, Analysis of technical variables heart vs. non-heart). Significant tissue-specific interactions were To assess whether the study’s biological variables of interest—tis- detected using the adaptive shrinkage method, ashr (Stephens sue and species—were confounded with the study’s recorded sam- 2017). Specifically, for each test, we input the regression estimates ple and technical variables, we used a previously described from limma to ashr: regression coefficients, posterior standard er- approach (Blake et al. 2018b). rors, and posterior degrees of freedom. We used the default settings For the 12 RNA-seq related technical variables that were the in ashr to calculate the shrunken regression coefficients (called most highly correlated with tissue or species, we assessed which the “posterior mean” in ashr), false discovery rate (FDR, also technical variables constitute the “best set” of independent vari- known as Q-value), and false sign rate (FSR, also known as S-value: ables to be included in a linear model for gene expression levels. the probability that sign of the estimated effect size is wrong in ei- Because of the partial correlations between the variables, we ap- ther direction). We assigned directionality based on the sign of the plied lasso regression using the package “glmnet” (Friedman posterior mean and determined significance based on the FSR. et al. 2010). Before performing the analysis, we also protected our variables of interest, tissue and species, in the model for each gene. We summarized each technical variable’s influence across The impact of matched tissue samples on DE results the genes by counting the number of times each technical variable To determine the impact of matched tissue samples on DE results, was included in the best set of the gene models. We found that we compared inter-tissue DE analysis results when using tissues none of the technical variables appeared in more than 25% of from the same or different individuals in GTEx v7 data (The the best sets (i.e., >25% of the gene models). Therefore, we chose GTEx Consortium 2017). We first subset the GTEx raw gene ex- not to include these technical variables in our model for testing pression count data to only individuals for which there was gene differential expression. count information in the heart and lung tissues, for genes included Finally, during our analysis of technical factors, we discovered in all three tissues. (The most GTEx samples were in heart and that RNA extraction date was confounded with species. In 2012, lung, so we decided to focus on these samples.) Furthermore, to we extracted RNA from the chimpanzee samples on March 8, minimize the number of covariates needed in the linear model, from the human samples on 3 d between March 12 and 29, and we decided to only analyze individuals of the same sex and whose from the rhesus samples on March 6. To test the relationship be- samples were sequenced on the same platform (sex = 1 and plat- tween the date of RNA extraction and gene expression PCs in hu- form = 1 from the GTEx documentation). We then normalized mans, we performed individual linear models on PCs 1–5 using the data and performed DE analysis using a voom + limma pipe- RNA extraction date as a predictor. None of the models were stat- line. In the linear model, we included tissue and three GTEx-pro- istically significant at FDR 10%, suggesting that tissue type is vided covariates (covariates 1 and 2 and inferred covariate 1 from more highly associated with gene expression levels than RNA ex- the covariate file for each tissue) as fixed effects and individual traction date. as a random effect. We chose to renormalize the raw counts data rather than use the normalized counts from GTEx because the voom + limma pipeline requires raw counts to assign voom Differential expression analysis using a linear model-based weights. We considered the output of the inter-tissue DE analysis framework for all individuals (DE vs. non-DE genes at FDR 5%) as the “ground To perform differential expression analysis, we used the same ap- truth.” To evaluate the impact of our study design, we then subset proach as Blake et al. (2018b). We applied a linear model-based em- the gene expression information to the individuals for which there pirical Bayes method (Smyth 2004; Smyth et al. 2005) that is information in all three tissues. We obtained gene expression accounts for the mean-variance relationship of the RNA-sequenc- level information from four randomly selected individuals and ing read counts, using weights specific to both genes and samples used the voom + limma pipeline to identify inter-tissue DE genes. (Law et al. 2014). Next, we compared the list of DE genes from this analysis to the To be considered a “tissue-specific DE gene” under our strin- ground truth list. We performed this downsampling procedure gent definition, the gene must be in the same direction and statisti- for tissues from the same four individuals as well as from four dif- cally significant in all pairwise comparisons including the given ferent individuals 10 times each and compared the number of true tissue but not significant in any comparison without that tissue. and false positives from the tests. For the analysis with eight differ- For example, for a gene to be classified as having heart-specific ent individuals, there were no repeated individuals, so we did not up-regulation in a given species, the gene needed to be up-regulat- use the “duplicateCorrelation” function in voom. ed (a significant, positive effect size) in heart versus liver, heart ver- sus lung, heart versus kidney, but not significantly different between the liver versus lung, liver versus kidney, or kidney versus BS-seq library preparation, sequencing, and mapping lung, in the same species. Under the more lenient definition We prepared a total of 48 whole-genome BS-seq libraries from ex- of tissue-specific DE genes, we compared the gene expression tracted DNA as previously described (Tung et al. 2012; Banovich level of one tissue to the mean of the other three tissues. To do et al. 2014). We aligned the trimmed reads to the human (hg19, so, we grouped the three tissues together and again used the February 2009), chimpanzee (panTro3, October 2010), or rhesus 258 Genome Research www.genome.org Downloaded from genome.cshlp.org on January 20, 2021 - Published by Cold Spring Harbor Laboratory Press Gene regulatory patterns across primate tissues macaque (rheMac2, January 2006) genomes, and to the lambda first had to determine the orthologous TSSs. We began with the phage genome using the Bismark aligner (version 0.8.1) (Krueger 12,184 orthologous genes in our RNA-sequencing analysis. Of and Andrews 2011). these, we found that 11,131 of these orthologous genes had an We estimated the percentage of methylation at an individual hg19 RefSeq TSS annotation (https://sourceforge.net/projects/ cytosine site by the ratio of the number of cytosines (unconverted) seqminer/files/Reference%20coordinate/refGene_hg19_TSS.bed/ found in mapped reads at this position, to the total number of download). We used liftOver to find orthologous sites in the chim- reads covering this position (sequenced as cytosine or thymine, panzees and rhesus macaque genomes in 9682 of those 11,131 i.e., converted or unconverted) using the methylation extractor genes. We then determined which of the hg19 RefSeq TSS annota- tool from Bismark (version 0.8.1). We additionally collapsed infor- tions were closest to the first hg19 orthologous exon, and we mation from both DNA strands (because CpG methylation status is repeated this process with the other two species and their respec- highly symmetrical on opposite strand) (Lister et al. 2009) to tive genomes. We found that 9604 of 9682 of the closest TSS anno- achieve better precision in methylation estimates across the tations in humans had the same liftOver coordinates in the other genome. two species. We then calculated the distance between the first To obtain CpGs that mapped to multiple species, the chim- orthologous exon to the TSS site in all three species individually. panzee and rhesus macaque CpG sites were mapped to human To minimize this difference between the three species, we filtered coordinates (hg19) using chain files from ftp://hgdownload all genes with a maximum distance difference across the species of .cse.ucsc.edu/goldenPath/hg19/liftOver/ and the liftOver tool larger than 2500 bp (for reference, the 75th percentile of the max- from the UCSC Genome Browser (Karolchik et al. 2014). These files imum difference in distance was 2078 bp). After this filtering step, had previously been filtered for paralogous regions and repeats, 7263 autosomal genes remained, and 4155 genes had at least two but we also removed positions that were not remapped to their orthologous CpGs 250 bp upstream of and 250 bp downstream original position when we mapped from human back to their orig- from the orthologous TSS. We chose a 250-bp window around inal genome. Chimpanzee and rhesus macaque CpG sites were the TSS based on DNA methylation levels around the promoter mapped to human, even if their orthologous positions were not in Banovich et al. (2014) and calculated the average of orthologous a CpG site in human. CpGs within this window for the 4155 genes. Using the same method but in humans and chimpanzees only, we found and cal- culated the average of orthologous CpGs within this window for Identifying differentially methylated regions (DMRs) 7725 genes. We were interested in identifying regions showing consistent dif- ferences between pairs of tissues or pairs of species, taking biolog- ical variation into account. To identify DMRs we used the linear model-based framework in the Bioconductor package bsseq (ver- Joint analysis of promoter DNA methylation and gene expression sion 0.10.0) (Hansen et al. 2012). For a given pairwise comparison levels (e.g., human liver vs. human heart), the bsseq package produces a To determine whether DNA methylation may underlie interspe- signal-to-noise statistic for each CpG site similar to a t-test statistic, cies differences in gene expression levels, we used a joint analysis assuming that methylation levels in each condition have equal method as described below. For each gene, we analyzed the gene variance. As recommended by Hansen et al. (2012), we used a expression levels, along with the accompanying average methyla- low-frequency mean correction to improve the marginal distribu- tion levels 250 bp upstream of and downstream from the TSS tion of the t-statistics. Similar to previous studies using this meth- (found above). For a given tissue, we first determined the effect odology, a t-statistic cutoff of −4.6, 4.6 was used for significance of species on gene expression levels using a linear model, with spe- (Hansen et al. 2011, 2014). DMRs were defined as regions contain- cies and RIN score as fixed effects (Model 1). Next, we parameter- ing at least three consecutive significant CpGs, an average methyl- ized a linear model attempting to predict expression levels ation difference of 10% between conditions, and at least one CpG exclusively from methylation levels. We refer to these residuals every 300 bp (Hansen et al. 2012). We used BEDTools (version as “methylation-corrected” gene expression values. We then 2.26.0) (Quinlan and Hall 2010) to calculate the number of over- used these values to again determine the effect of species, this lapping DMRs across tissues and/or species (Blake et al. 2018a). time on gene expression levels “corrected” for methylation, using We required overlapping DMRs to have a minimum base pair over- a linear model with species and RIN score as fixed effects (Model 2). lap of at least 25%, unless otherwise stated. To determine the contribution of DNA methylation levels to inter- To be considered a tissue-specific DMR, the region was re- species differences in gene expression, we computed the difference quired to be a significant tissue DMR (T-DMR) in one tissue com- in the species effect size between Model 1 and Model 2 for each pared to the other three tissues pairwise (in the same direction) gene, as well as the standard error of the difference. Large effect but not a significant T-DMR across any of the other three tissues size differences between Models 1 and 2 for a given gene suggest in pairwise comparisons. We again used BEDTools to ensure a that methylation status may be a significant driver of DE. To assess minimum base pair overlap of 25%. Once the tissue-specific the significance of this difference, we used adaptive shrinkage DMRs were identified within a species, we then classified them (ashr) (Stephens 2017) to compute the posterior mean of the dif- as species-specific or conserved. To be considered conserved ferences in the effect sizes, using vashr (Lu and Stephens 2016), (across humans and chimpanzees or all three species), the tissue- with the degrees of freedom equal to the number of samples in specific DMR had to be significant in all species in the comparison the linear model minus 2. The shrunken variances from vashr and have a minimum base pair overlap of the T-DMR of at were used in the ashr posterior mean computation. From this least 25%. procedure, we obtained the number of genes where species has a significant difference in effect sizes before and after regressing out methylation. We assessed significance using the S-value Calculating the average methylation levels of conserved statistic (FS) (Stephens 2017). Using the S-values, rather than the promoters Q-values, not only takes significance into account but also has To calculate the DNA methylation levels of orthologous CpGs the added benefit of assessing our confidence in the direction of around the transcription start site (TSS) of orthologous genes, we the effect. Genome Research 259 www.genome.org Downloaded from genome.cshlp.org on January 20, 2021 - Published by Cold Spring Harbor Laboratory Press Blake et al. We performed the above analysis separately for interspecies References DE genes and non-DE genes, and in each tissue individually. We Andersson R, Gebhard C, Miguel-Escalada I, Hoof I, Bornholdt J, Boyd M, identified interspecies DE genes in our tissue of interest as those Chen Y, Zhao X, Schmidl C, Suzuki T, et al. 2014. An atlas of active en- with a significant species term in the model of species and RIN hancers across human cell types and tissues. Nature 507: 455–461. score as fixed effects. We assessed significance of DE genes at doi:10.1038/nature12787 FDR 5%, unless otherwise noted. Babbitt CC, Fedrigo O, Pfefferle AD, Boyle AP, Horvath JE, Furey TS, Wray GA. 2010. Both noncoding and protein-coding RNAs contribute to We also applied the same analysis framework to determine gene expression evolution in the primate brain. Genome Biol Evol 2: whether DNA methylation may underlie inter-tissue differences 67–79. doi:10.1093/gbe/evq002 in gene expression levels. For the inter-tissue DE genes and non- Balashova V, Abdulkadyrov K. 1984. Cellular composition of hemopoietic DE genes, we replaced “species” with “tissue” as a fixed effect in tissue of the liver and spleen in the human fetus. Arkh Anat 86: 80. models 1 and 2. We assessed significance with various FDR and Banovich NE, Lan X, McVicker G, van de Geijn B, Degner JF, Blischak JD, Roux J, Pritchard JK, Gilad Y. 2014. Methylation QTLs are associated FSR thresholds, as specified in the text. with coordinated changes in transcription factor binding, histone mod- ifications, and gene expression levels. PLoS Genet 10: e1004663. doi:10 .1371/journal/pgen.1004663 Barbash S, Sakmar TP. 2017. Brain gene expression signature on primate ge- Data access nomic sequence evolution. Sci Rep 7: 17329. doi:10.1038/s41598-017- All raw and processed sequencing data generated in this study have 17462-3 Barbosa-Morais NL, Irimia M, Pan Q, Xiong HY, Gueroussov S, Lee LJ, been submitted to the NCBI Gene Expression Omnibus (GEO; Slobodeniuc V, Kutter C, Watt S, Colak R, et al. 2012. The evolutionary https://www.ncbi.nlm.nih.gov/geo/) under accession number landscape of alternative splicing in vertebrate species. Science 338: GSE112356. Data and scripts used in this paper are available at 1587–1593. doi:10.1126/science.1230612 https://github.com/Lauren-Blake/Reg_Evo_Primates and in the Barlow DP. 1993. Methylation and imprinting: from host defense to gene regulation? Science 260: 309–310. doi:10.1126/science.8469984 Supplemental Code. The results of our scripts can be viewed at Blake LE, Roux J, Gilad Y. 2018a. Pairwise DMRs for, “A comparison of gene https://lauren-blake.github.io/Regulatory_Evol/analysis/. expression and DNA methylation patterns across tissues and species” (Version 1). doi:10.5281/zenodo.1495829 Blake LE, Thomas SM, Blischak JD, Hsiao CJ, Chavarria C, Myrthil M, Gilad Y, Pavlovic BJ. 2018b. A comparative study of endoderm differentiation Competing interest statement in humans and chimpanzees. Genome Biol 19: 162. doi:10.1186/s13059- 018-1490-5 The authors declare no competing interests. Blekhman R. 2012. A database of orthologous exons in primates for compar- ative analysis of RNA-seq data. Nat Preced doi:hdl:10101/npre.2012 .7054.1 Acknowledgments Blekhman R, Oshlack A, Chabot AE, Smyth GK, Gilad Y. 2008. Gene regu- lation in primates evolves under tissue-specific selection pressures. Members of the Gilad, Marques-Bonet, Robinson-Rechavi, Ste- PLoS Genet 4: e1000271. doi:10.1371/journal.pgen.1000271 Blekhman R, Marioni JC, Zumbo P, Stephens M, Gilad Y. 2010. Sex-specific phens, and Pritchard laboratories provided helpful discussions and lineage-specific alternative splicing in primates. Genome Res 20: and comments on the manuscript. In particular, Matthew Ste- 180–189. doi:10.1101/gr.099226.109 phens provided guidance on integrating the DNA methylation Brawand D, Soumillon M, Necsulea A, Julien P, Csárdi G, Harrigan P, Weier and gene expression data, and Michelle Ward and Natalia Gonza- M, Liechti A, Aximu-Petri A, Kircher M, et al. 2011. The evolution of gene expression levels in mammalian organs. Nature 478: 343–348. les provided helpful comments on a draft of the manuscript. We doi:10.1038/nature10532 thank Kasper Hansen, Jenny Tung, and Luis Barreiro for discus- Brodal A. 1983. The perihypoglossal nuclei in the macaque monkey and the sions regarding multispecies methylation analysis. We acknowl- chimpanzee. J Comp Neurol 218: 257–269. doi:10.1002/cne.902180303 edge Athma Pai for sharing a methylation protocol, Bryce van de Cáceres M, Lachuer J, Zapala MA, Redmond JC, Kudo L, Geschwind DH, Geijn for help with the WASP pipeline, and Charlotte Soneson, Lockhart DJ, Preuss TM, Barlow C. 2003. Elevated gene expression levels distinguish human from non-human primate brains. Proc Natl Acad Sci Jacob Degner, and Unjin Lee. We also thank the Yerkes Primate 100: 13030–13035. doi:10.1073/pnas.2135499100 Center and Southwest Foundation for Biomedical Research, Cain CE, Blekhman R, Marioni JC, Gilad Y. 2011. Gene expression differenc- Anne Stone, and Jéssica Hernández Rodríguez for providing and/ es among primates are associated with changes in a histone epigenetic or helping to ship the tissue samples. L.E.B. was supported by modification. Genetics 187: 1225–1234. doi:10.1534/genetics.110 .126177 the National Science Foundation Graduate Research Fellowship Capra JA, Erwin GD, McKinsey G, Rubenstein JL, Pollard KS. 2013. Many (DGE-1144082). Additionally, L.E.B. and I.E. were funded by the human accelerated regions are developmental enhancers. Philos Trans National Institutes of Health (NIH) Genetics and Gene Regulation R Soc Lond B Biol Sci 368: 20130025. doi:10.1098/rstb.2013.0025 Training Grant (T32 GM07197). J.R. was funded by a Swiss Nation- Chen J, Swofford R, Johnson J, Cummings BB, Rogel N, Lindblad-Toh K, Haerty W, Palma FD, Regev A. 2019. A quantitative framework for char- al Science Foundation postdoc mobility fellowship (PBLAP3- acterizing the evolutionary history of mammalian gene expression. 134342) and the FP7 Marie Curie International Outgoing Fellow- Genome Res 29: 53–63. doi:10.1101/gr.237636.118 ship PRIMATE_REG_ EVOL. This project was funded in part by Cooper GM, Shendure J. 2011. Needles in stacks of needles: finding disease- the Office of Research Infrastructure/Office of the Director causal variants in a wealth of genomic data. Nat Rev Genet 12: 628–640. doi:10.1038/nrg3046 (ORIP/OD) P51OD011132 grant. T.M.-B. is supported by Dezső Z, Nikolsky Y, Sviridov E, Shi W, Serebriyskaya T, Dosymbekov D, BFU2017-86471-P (Ministry of Economy and Competitiveness/ Bugrim A, Rakhmatulin E, Brennan RJ, Guryanov A, et al. 2008. A com- European Regional Development Fund, European Union), NIH prehensive functional analysis of tissue specificity of human gene ex- U01 MH106874 grant, Howard Hughes Medical Institute: Interna- pression. BMC Biol 6: 49. doi:10.1186/1741-7007-6-49 tional Early Career, Obra Social “La Caixa” and Secretaria d’Univer- Enard W, Khaitovich P, Klose J, Zöllner S, Heissig F, Giavalisco P, Nieselt- Struwe K, Muchmore E, Varki A, Ravid R, et al. 2002. Intra- and interspe- sitats i Recerca and Centres de Recerca de Catalunya Programme cific variation in primate gene expression patterns. Science 296: 340– del Departament d’Economia i Coneixement de la Generalitat de 343. doi:10.1126/science.1068996 Catalunya. The content presented in this article is solely the re- The ENCODE Project Consortium. 2012. An integrated encyclopedia of sponsibility of the investigators and does not necessarily reflect DNA elements in the human genome. Nature 489: 57–74. doi:10 .1038/nature11247 the official views of the funders. The funders had no role in study Friedman J, Hastie T, Tibshirani R. 2010. Regularization paths for general- design, data collection and analysis, decision to publish, or prepa- ized linear models via coordinate descent. J Stat Softw 33: 1–22. doi:10 ration of the manuscript. .18637/jss.v033.i01 260 Genome Research www.genome.org Downloaded from genome.cshlp.org on January 20, 2021 - Published by Cold Spring Harbor Laboratory Press Gene regulatory patterns across primate tissues Gallego Romero I, Pavlovic BJ, Hernando-Herraez I, Zhou X, Ward MC, Lindskog C. 2016. The Human Protein Atlas—an important resource for ba- Banovich NE, Kagan CL, Burnett JE, Huang CH, Mitrano A, et al. sic and clinical research. Expert Rev Proteomics 13: 627–629. doi:10.1080/ 2015. A panel of induced pluripotent stem cells from chimpanzees: a re- 14789450.2016.1199280 source for comparative functional genomics. eLife 4: e07103. doi:10 Lister R, Pelizzola M, Dowen RH, Hawkins RD, Hon G, Tonti-Filippini J, Nery .7554/eLife.07103 JR, Lee L, Ye Z, Ngo QM, et al. 2009. Human DNA methylomes at base Gilad Y, Mizrahi-Man O. 2015. A reanalysis of mouse ENCODE comparative resolution show widespread epigenomic differences. Nature 462: 315– gene expression data. F1000 Res 4: 121. doi:10.12688/f11000research 322. doi:10.1038/nature08514 .16536.12681 Lu M, Stephens M. 2016. Variance adaptive shrinkage (vash): flexible empir- Gilad Y, Rifkin SA, Bertone P, Gerstein M, White KP. 2005. Multi-species mi- ical Bayes estimation of variances. Bioinformatics 32: 3428–3434. doi:10 croarrays reveal the effect of sequence divergence on gene expression .1093/bioinformatics/btw483 profiles. Genome Res 15: 674–680. doi:10.1101/gr.3335705 Marioni JC, Mason CE, Mane SM, Stephens M, Gilad Y. 2008. RNA-seq: an The GTEx Consortium. 2017. Genetic effects on gene expression across hu- assessment of technical reproducibility and comparison with gene ex- man tissues. Nature 550: 204–213. doi:10.1038/nature24277 pression arrays. Genome Res 18: 1509–1517. doi:10.1101/gr.079558.108 Hansen KD, Timp W, Bravo HC, Sabunciyan S, Langmead B, McDonald OG, Martin M. 2011. Cutadapt removes adapter sequences from high- Wen B, Wu H, Liu Y, Diep D, et al. 2011. Increased methylation varia- throughput sequencing reads. EMBnet J 17: 10. doi:10.14806/ej.17.1 tion in epigenetic domains across cancer types. Nat Genet 43: 768– .200 775. doi:10.1038/ng.865 Martin DI, Singer M, Dhahbi J, Mao G, Zhang L, Schroth GP, Pachter L, Boffelli D. 2011. Phyloepigenomic comparison of great apes reveals a Hansen KD, Langmead B, Irizarry RA. 2012. BSmooth: from whole genome correlation between somatic and germline methylation states. Genome bisulfite sequencing reads to differentially methylated regions. Genome Res 21: 2049–2057. doi:10.1101/gr.122721.111 Biol 13: R83. doi:10.1186/gb-2012-13-10-r83 Merkin J, Russell C, Chen P, Burge CB. 2012. Evolutionary dynamics of gene Hansen KD, Sabunciyan S, Langmead B, Nagy N, Curley R, Klein G, Klein E, and isoform regulation in mammalian tissues. Science 338: 1593–1599. Salamon D, Feinberg AP. 2014. Large-scale hypomethylated blocks asso- doi:10.1126/science.1228186 ciated with Epstein-Barr virus–induced B-cell immortalization. Genome Molaro A, Hodges E, Fang F, Song Q, McCombie WR, Hannon GJ, Smith AD. Res 24: 177–184. doi:10.1101/gr.157743.113 2011. Sperm methylation profiles reveal features of epigenetic inheri- Hernando-Herraez I, Prado-Martinez J, Garg P, Fernandez-Callejo M, Heyn tance and evolution in primates. Cell 146: 1029–1041. doi:10.1016/j H, Hvilsom C, Navarro A, Esteller M, Sharp AJ, Marques-Bonet T. .cell.2011.08.016 2013. Dynamics of DNA methylation in recent human and great ape Nowick K, Gernat T, Almaas E, Stubbs L. 2009. Differences in human and evolution. PLoS Genet 9: e1003763. doi:10.1371/journal.pgen.1003763 chimpanzee gene expression patterns define an evolving network of Hernando-Herraez I, Garcia-Perez R, Sharp AJ, Marques-Bonet T. 2015a. transcription factors in brain. Proc Natl Acad Sci 106: 22358–22363. DNA methylation: insights into human evolution. PLoS Genet 11: doi:10.1073/pnas.0911376106 e1005661. doi:10.1371/journal.pgen.1005661 Pai AA, Bell JT, Marioni JC, Pritchard JK, Gilad Y. 2011. A genome-wide Hernando-Herraez I, Heyn H, Fernandez-Callejo M, Vidal E, Fernandez- study of DNA methylation patterns and gene expression levels in mul- Bellon H, Prado-Martinez J, Sharp AJ, Esteller M, Marques-Bonet T. tiple human and chimpanzee tissues. PLoS Genet 7: e1001316. doi:10 2015b. The interplay between DNA methylation and sequence diver- .1371/journal.pgen.1001316 gence in recent human evolution. Nucleic Acids Res 43: 8204–8214. Pierson E, GTEx Consortium, Koller D, Battle A, Mostafavi S, Ardlie KG, Getz doi:10.1093/nar/gkv693 G, Wright FA, Kellis M, Volpi S, et al. 2015. Sharing and specificity of co- Karaman MW, Houck ML, Chemnick LG, Nagpal S, Chawannakul D, expression networks across 35 human tissues. PLoS Comput Biol 11: Sudano D, Pike BL, Ho VV, Ryder OA, Hacia JG. 2003. Comparative anal- e1004220. doi:10.1371/journal.pcbi.1004220 ysis of gene-expression patterns in human and African great ape cul- Quinlan AR, Hall IM. 2010. BEDTools: a flexible suite of utilities for compar- tured fibroblasts. Genome Res 13: 1619–1630. doi:10.1101/gr.1289803 ing genomic features. Bioinformatics 26: 841–842. doi:10.1093/bioinfor Karolchik D, Barber GP, Casper J, Clawson H, Cline MS, Diekhans M, matics/btq033 Dreszer TR, Fujita PA, Guruvadoo L, Haeussler M, et al. 2014. The Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. 2015. UCSC Genome Browser database: 2014 update. Nucleic Acids Res 42: limma powers differential expression analyses for RNA-sequencing D764–D770. doi:10.1093/nar/gkt1168 and microarray studies. Nucleic Acids Res 43: e47. doi:10.1093/nar/ Khaitovich P, Muetzel B, She X, Lachmann M, Hellmann I, Dietzsch J, gkv007 Steigele S, Do HH, Weiss G, Enard W, et al. 2004a. Regional patterns Romero IG, Ruvinsky I, Gilad Y. 2012. Comparative studies of gene expres- of gene expression in human and chimpanzee brains. Genome Res 14: sion and the evolution of gene regulation. Nat Rev Genet 13: 505–516. 1462–1473. doi:10.1101/gr.2538704 doi:10.1038/nrg3229 Khaitovich P, Weiss G, Lachmann M, Hellmann I, Enard W, Muetzel B, Sharp AJ, Stathaki E, Migliavacca E, Brahmachary M, Montgomery SB, Wirkner U, Ansorge W, Pääbo S. 2004b. A neutral model of transcrip- Dupre Y, Antonarakis SE. 2011. DNA methylation profiles of human ac- tome evolution. PLoS Biol 2: e132. doi:10.1371/journal.pbio.0020132 tive and inactive X chromosomes. Genome Res 21: 1592–1600. doi:10 Khaitovich P, Hellmann I, Enard W, Nowick K, Leinweber M, Franz H, Weiss .1101/gr.112680.110 G, Lachmann M, Pääbo S. 2005. Parallel patterns of evolution in the ge- Shibata Y, Sheffield NC, Fedrigo O, Babbitt CC, Wortham M, Tewari AK, nomes and transcriptomes of humans and chimpanzees. Science 309: London D, Song L, Lee BK, Iyer VR, et al. 2012. Extensive evolutionary 1850–1854. doi:10.1126/science.1108296 changes in regulatory element activity during human origins are associ- Khan Z, Ford MJ, Cusanovich DA, Mitrano A, Pritchard JK, Gilad Y. 2013. ated with altered gene expression and positive selection. PLoS Genet 8: Primate transcript and protein expression levels evolve under compen- e1002789. doi:10.1371/journal.pgen.1002789 Smyth GK. 2004. Linear models and empirical Bayes methods for assessing satory selection pressures. Science 342: 1100–1104. doi:10.1126/science differential expression in microarray experiments. Stat Appl Genet Mol .1242379 Biol 3: Article3. doi:10.2202/1544-6115.1027 Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. 2013. Smyth GK, Michaud J, Scott HS. 2005. Use of within-array replicate spots for TopHat2: accurate alignment of transcriptomes in the presence of inser- assessing differential expression in microarray experiments. tions, deletions and gene fusions. Genome Biol 14: R36. doi:10.1186/gb- Bioinformatics 21: 2067–2075. doi:10.1093/bioinformatics/bti270 2013-14-4-r36 Stephens M. 2017. False discovery rates: a new deal. Biostatistics 18: 275– King MC, Wilson AC. 1975. Evolution at two levels in humans and chim- 294. doi:10.1093/biostatistics/kxw041 panzees. Science 188: 107–116. doi:10.1126/science.1090005 Stergachis AB, Neph S, Sandstrom R, Haugen E, Reynolds AP, Zhang M, Krueger F, Andrews SR. 2011. Bismark: a flexible aligner and methylation Byron R, Canfield T, Stelhing-Sun S, Lee K, et al. 2014. Conservation caller for Bisulfite-Seq applications. Bioinformatics 27: 1571–1572. of trans-acting circuitry during mammalian regulatory evolution. doi:10.1093/bioinformatics/btr167 Nature 515: 365–370. doi:10.1038/nature13972 Law CW, Chen Y, Shi W, Smyth GK. 2014. voom: precision weights unlock Stuart JM, Segal E, Koller D, Kim SK. 2003. A gene-coexpression network for linear model analysis tools for RNA-seq read counts. Genome Biol 15: global discovery of conserved genetic modules. Science 302: 249–255. R29. doi:10.1186/gb-2014-15-2-r29 doi:10.1126/science.1087447 Lemos B, Meiklejohn CD, Cáceres M, Hartl DL. 2005. Rates of divergence in Supek F, Bošnjak M, Škunca N, Šmuc T. 2011. REVIGO summarizes and vis- gene expression profiles of primates, mice, and flies: stabilizing selection ualizes long lists of gene ontology terms. PLoS One 6: e21800. doi:10 and variability among functional categories. Evolution 59: 126–137. .1371/journal.pone.0021800 doi:10.1554/04-251 Thul PJ, Lindskog C. 2018. The human protein atlas: a spatial map of the hu- Lin S, Lin Y, Nery JR, Urich MA, Breschi A, Davis CA, Dobin A, Zaleski C, Beer man proteome. Protein Sci 27: 233–244. doi:10.1002/pro.3307 MA, Chapman WC, et al. 2014. Comparison of the transcriptional land- Tsankov AM, Gu H, Akopian V, Ziller MJ, Donaghey J, Amit I, Gnirke A, scapes between human and mouse tissues. Proc Natl Acad Sci 111: Meissner A. 2015. Transcription factor binding dynamics during human 17224–17229. doi:10.1073/pnas.1413624111 ES cell differentiation. Nature 518: 344–349. doi:10.1038/nature14233 Genome Research 261 www.genome.org Downloaded from genome.cshlp.org on January 20, 2021 - Published by Cold Spring Harbor Laboratory Press Blake et al. Tung J, Barreiro LB, Johnson ZP, Hansen KD, Michopoulos V, Toufexis D, and chimpanzees. PLoS One 7: e34935. doi:10.1371/journal.pone Michelini K, Wilson ME, Gilad Y. 2012. Social environment is associated .0034935 with gene regulatory variation in the rhesus macaque immune system. Yu NY, Hallström BM, Fagerberg L, Ponten F, Kawaji H, Carninci P, Proc Natl Acad 109: 6490–6495. doi:10.1073/pnas.1202734109 Forrest AR, The FANTOM Consortium, Hayashizaki Y, Uhlén M, et al. Uhlén M, Fagerberg L, Hallström BM, Lindskog C, Oksvold P, Mardinoglu A, 2015. Complementing tissue characterization by integrating Sivertsson Å, Kampf C, Sjöstedt E, Asplund A, et al. 2015. Proteomics. transcriptome profiling from the Human Protein Atlas and from the Tissue-based map of the human proteome. Science 347: 1260419. FANTOM5 consortium. Nucleic Acids Res 43: 6787–6798. doi:10.1093/ doi:10.1126/science.1260419 nar/gkv608 Villar D, Berthelot C, Aldridge S, Rayner TF, Lukk M, Pignatelli M, Park TJ, Zhang B, Horvath S. 2005. A general framework for weighted gene co-ex- Deaville R, Erichsen JT, Jasinska AJ, et al. 2015. Enhancer evolution pression network analysis. Stat Appl Genet Mol Biol 4: Article17. doi:10 across 20 mammalian species. Cell 160: 554–566. doi:10.1016/j.cell .2202/1544-6115.1128 .2015.01.006 Zhou X, Cain CE, Myrthil M, Lewellen N, Michelini K, Davenport ER, Ward MC, Wilson MD, Barbosa-Morais NL, Schmidt D, Stark R, Pan Q, Stephens M, Pritchard JK, Gilad Y. 2014. Epigenetic modifications are Schwalie PC, Menon S, Lukk M, Watt S, et al. 2013. Latent regulatory po- associated with inter-species gene expression variation in primates. tential of human-specific repetitive elements. Mol Cell 49: 262–272. Genome Biol 15: 547. doi:10.1186/s13059-014-0547-3 doi:10.1016/j.molcel.2012.11.013 Warnefors M, Eyre-Walker A. 2012. A selection index for gene expression evolution and its application to the divergence between humans Received July 18, 2019; accepted in revised form January 2, 2020. 262 Genome Research www.genome.org Downloaded from genome.cshlp.org on January 20, 2021 - Published by Cold Spring Harbor Laboratory Press A comparison of gene expression and DNA methylation patterns across tissues and species Lauren E. Blake, Julien Roux, Irene Hernando-Herraez, et al. Genome Res. 2020 30: 250-262 originally published online January 17, 2020 Access the most recent version at doi:10.1101/gr.254904.119 Supplemental http://genome.cshlp.org/content/suppl/2020/02/03/gr.254904.119.DC1 Material References This article cites 77 articles, 23 of which can be accessed free at: http://genome.cshlp.org/content/30/2/250.full.html#ref-list-1 Open Access Freely available online through the Genome Research Open Access option. Creative This article, published in Genome Research, is available under a Creative Commons Commons License (Attribution-NonCommercial 4.0 International), as described at License http://creativecommons.org/licenses/by-nc/4.0/. Email Alerting Receive free email alerts when new articles cite this article - sign up in the box at the Service top right corner of the article or click here. To subscribe to Genome Research go to: http://genome.cshlp.org/subscriptions © 2020 Blake et al.; Published by Cold Spring Harbor Laboratory Press