DNA methylation QTL mapping across diverse human tissues provides molecular links between genetic variation and complex traits

Studies of DNA methylation (DNAm) in solid human tissues are relatively scarce; tissue-specific characterization of DNAm is needed to understand its role in gene regulation and its relevance to complex traits. We generated array-based DNAm profiles for 987 human samples from the Genotype-Tissue Expression (GTEx) project, representing 9 tissue types and 424 subjects. We characterized methylome and transcriptome correlations (eQTMs), genetic regulation in cis (mQTLs and eQTLs) across tissues and e/mQTLs links to complex traits. We identified mQTLs for 286,152 CpG sites, many of which (>5%) show tissue specificity, and mQTL colocalizations with 2,254 distinct GWAS hits across 83 traits. For 91% of these loci, a candidate gene link was identified by integration of functional maps, including eQTMs, and/or eQTL colocalization, but only 33% of loci involved an eQTL and mQTL present in the same tissue type. With this DNAm-focused integrative analysis, we contribute to the understanding of molecular regulatory mechanisms in human tissues and their impact on complex traits. As part of the enhanced GTEx (eGTEx) project, 987 human samples from 9 tissue types and 424 donors are assayed using DNA methylation microarrays. Colocalization of GWAS variants, eQTLs and mQTLs shows diverse links between genetic variation, molecular phenotypes and complex traits.

Studies of DNA methylation (DNAm) in solid human tissues are relatively scarce; tissue-specific characterization of DNAm is needed to understand its role in gene regulation and its relevance to complex traits. We generated array-based DNAm profiles for 987 human samples from the Genotype-Tissue Expression (GTEx) project, representing 9 tissue types and 424 subjects. We characterized methylome and transcriptome correlations (eQTMs), genetic regulation in cis (mQTLs and eQTLs) across tissues and e/ mQTLs links to complex traits. We identified mQTLs for 286,152 CpG sites, many of which (>5%) show tissue specificity, and mQTL colocalizations with 2,254 distinct GWAS hits across 83 traits. For 91% of these loci, a candidate gene link was identified by integration of functional maps, including eQTMs, and/or eQTL colocalization, but only 33% of loci involved an eQTL and mQTL present in the same tissue type. With this DNAm-focused integrative analysis, we contribute to the understanding of molecular regulatory mechanisms in human tissues and their impact on complex traits.
The majority of common genetic variants that impact human traits are believed to exert their effects through the regulation of gene expression 1,2 . Genetic regulation of gene expression has been comprehensively characterized across many human tissue types by the GTEx project 3,4 , and expression quantitative trait loci (eQTLs) appear to underlie a substantial fraction of variant-trait associations [5][6][7] . However, our understanding of the regulatory mechanisms by which variants influence human traits is far from complete.

Genetic regulation of DNAm in cis exhibits tissue specificity
To characterize the genetic regulation of DNAm across tissues, we mapped genetic variants that affect DNAm levels of proximal (±500 kb) CpG sites in cis (mQTLs). To optimize mQTL discovery, we accounted for non-genetic DNAm variability by fitting per-tissue linear models including surrogate variables that capture technical, environmental and biological factors like age, cellular heterogeneity and smoking habits (Extended Data Figs. 2 and 3). We identified significant (false discovery rate (FDR) < 0.05) mQTL CpG sites in cis (mCpGs) and modeled mQTL effects jointly across tissues to increase QTL-mapping power by leveraging tissue-shared effects 42 . Additionally, we mapped secondary cis-QTL signals using a stepwise regression procedure 43 . Although mQTLs are substantially more abundant than eQTLs 20,44 , their distinctive genomic properties are not fully characterized. To facilitate comparisons of mQTL to cis-eQTL patterns, we used an analogous QTL-mapping approach to identify eQTLs.
We detected a total of 286,152 mQTL CpG sites, of which 45,543 (16%) have secondary signals in at least one tissue. Modeling mQTL effects jointly across tissues resulted in a substantial (more than fivefold) increase in detectable mCpGs per tissue, ranging from 108,844 in testis to 206,802 in lung (Fig. 2a and Supplementary Table 3). We quantified the mQTL replication rate in external muscle 20 and brain 23 datasets (Supplementary Table 3) and observed high replication rates (average π1 = 0.85).
For a particular CpG, genetic regulation of DNAm tends to be either highly tissue specific or highly shared across tissue types (Extended Data Fig. 4a,b), similar to patterns observed for eQTLs and mQTL catalogs derived from a variety of healthy, solid tissues are needed to contribute to the characterization of the etiology of complex traits.
Gene expression with subject-matched genotype data is available for many healthy tissue sources from hundreds of individuals 4,37 , but comparable data sources for DNAm are limited to fewer tissue types 38 . The enhancing GTEx (eGTEx) project 39 seeks to complement existing gene expression data from human tissues with additional molecular traits, including DNAm. Here, we profiled human DNAm for >750,000 genomic locations (CpGs or CpG sites) using the Illumina EPIC methylation array. The dataset includes 987 samples from 424 donors representing 9 tissue types: breast, kidney, colon, lung, muscle, ovary, prostate, testis and whole blood (Supplementary Table 1). Gene expression is also available for 3,872 samples derived from these 9 tissue types 4 , and 495 DNAm-profiled samples have gene expression and genotype data available. To further characterize the relationship between DNAm and gene expression in a tissue-specific manner, we identified CpGs for which DNAm was associated with local gene expression, that is expression quantitative trait methylation sites (eQTMs). To characterize the cis-genetic regulation of methylation across tissues and compare it to its gene expression counterpart, we mapped mQTLs and eQTLs. We contrasted mQTL to eQTL profiles and characterized their interdependency as well as differential tissue specificities, functional mechanisms and regulatory pleiotropy signatures. To evaluate the impact of mQTLs on human traits and assess their relative contribution compared to eQTLs, we integrated mQTLs and eQTLs with functional maps, including eQTMs, multi-context eQTLs and GWAS summary statistics of 87 GWASs. Overview of analyses performed is shown in Fig. 1.

DNAm and its correlation with local gene expression differs across tissues and by regulatory features
To investigate the contribution of tissue type to the similarity of DNAm-profiled samples, we applied dimension-reduction and clustering approaches (Methods), which revealed sample similarity by tissue of origin and comparable clusters to transcriptome-derived ones (Extended Data Fig. 1a,b). To characterize the interdependence  Fig. 4b). Compared to GTEx eQTLs, detected mQTLs, and corresponding effect sizes, are more likely to be shared across tissue types (Fig. 2b), as observed in blood cells 26 . On average, for a given mCpG, 46% of tissues show a similar (within a factor of 2) effect relative to the tissue with the strongest mQTL effect size, and effects are more similar across tissues for mCpGs located in proximal, compared to distal, gene regulatory regions (Extended Data Fig. 4c). On the other hand, a substantial fraction (5% on average) of the identified mCpGs in a given tissue were not detected in other tissues (Extended Data Fig. 4b).

Functional mechanisms that drive genetic regulation of DNAm differ from gene expression
Molecular mechanisms underlying interindividual DNAm variation are known to exist [45][46][47][48] , but their role in mQTLs, and overlap with eQTLs, is not fully characterized. To compare the molecular mechanisms of detected mQTLs relative to eQTLs, we integrated annotation of CGIs, genomic functional elements and chromatin states (Supplementary Table 4) and performed within-and meta-tissue enrichment analyses (Methods). We observe that, in contrast with eQTLs, the detected mQTLs are more likely to reside in regulatory regions that appear inactive in the mQTL-mapping context analyzed herein (Extended Data Fig. 6), in line with whole-genome mQTL studies 44 reporting that most mQTLs are located in quiescent genomic regions. Although both eQTLs and mQTLs are enriched in gene regulatory regions (Fig. 2c), only eQTLs are enriched in CGIs and gene transcripts, particularly in splicing and untranslated exon regions (UTRs), as previously described 4 .
Conversely, detected mQTLs are depleted from CGIs and genes, as previously observed in blood 21,31,34 , but are strongly enriched in distal enhancers and putative insulators (Fig. 2c). QTLs often originate from transcription factor (TF) binding alterations to TF binding sites (TFBSs) 49 . mQTLs can impact methylation of nearby CpGs by directly altering TF occupancy or by altering TF binding and mediating recruitment of TET1 and DNMT3A enzymes 8,25,28,32,50,51 , and several TFs have been associated with mQTLs 8,50 . Yet, the contribution of mQTL-relative to eQTL-associated TFBSs is not well characterized. To further characterize putative mQTL and eQTL TFBS, we integrated chromatin immunoprecipitation with sequencing-derived annotation of TFBSs corresponding to 339 TFs. We identified 126 TFBSs significantly enriched (FDR < 0.01, OR > 1) in eQTLs or mQTLs (Supplementary Table 4). We observed a distinctive TFBS enrichment profile for mQTLs compared to eQTLs. The eQTL enrichments with the smallest P values across tissues correspond to TFs involved in basal transcription (for example, RNA Polymerase II genes). In contrast, mQTLs were enriched, among other TFBSs, in binding sites of steroid receptors (for example, ESR1 and NR2F2) and other proteins known to be involved in 3D organization of the genome (for example, ZNF143) (Supplementary Table 3 and Fig. 2d).
Altogether, these results suggest that mQTLs and eQTLs, in part, differ in their underlying biological mechanisms, driven by distinct sets of TFs. Although eQTLs tend to result from variants residing in gene bodies and regulatory elements, mQTLs result in part from variants altering nongenic, distal regulatory elements, including insulators and elements bound by proteins involved in chromatin spatial conformation and long-range interactions.

Genetic co-regulation of DNAm and gene expression is relatively scarce
Given their divergent genomic enrichment profiles, it is expected that mQTLs and eQTLs are driven, in part, by different causal variants. To identify eQTL-mQTL pairs (e/mQTL) likely to share a common putative causal variant, we performed e/mQTL colocalization, and observed that the proportion of detectable e/mQTL colocalizations is moderate (Fig. 2e), as only 21% of mQTL loci are suggestively colocalized (PP4 > 0.5) with at least one eQTL. Despite limitations in accurately estimating this fraction, our results indicate that a considerable fraction of detected mQTLs do not show clear tissue-matching associations with local gene expression, as previously observed 20,23,32 .

Genetic regulation of DNAm is characterized by molecular regulatory pleiotropy
Given that correlated methylated CpGs tend to be close together 52 , a genetic variant influencing DNAm in cis is expected to display molecular pleiotropy as it often impacts multiple CpGs 53 . To characterize the molecular pleiotropic nature of mQTLs, we quantified the number of mCpGs and eGenes involved in mQTL-eQTL colocalizations (Methods). Overall, we observe pervasive pleiotropy; the majority (78%) of colocalized eQTLs-mQTLs impact multiple mCpGs, and a considerable minority (28%) impact multiple eGenes (Extended Data Fig. 7). The largest pleiotropic set, identified in ovary, is led by variant rs6433571 and involves 114 mCpGs and 8 eGenes in the HOXD gene cluster region (Fig. 3a), associated with epithelial ovarian cancer 54 . This pleiotropic effect is not driven by ovary-specific gene expression (Fig. 3b) but by ovary-specific genetic regulation of DNAm and expression (Fig. 3c,d). By means of QTL-GWAS colocalization, we identified 112/114 mCpGs and 7/8 eGenes as significantly (PP4 > 0.5) colocalized with ovarian cancer risk (Supplementary Table 5), including the HOXD1 and HOXD3 genes which have a suspected role in genetic risk of ovarian cancer 55,56 , as well as less characterized genes, for example HAGLR. Together, these findings illustrate how genetic variants can modify the DNAm landscape in a long-range and tissue-specific manner, with coordinated changes in gene expression and implications for disease risk.

Genetic regulation of DNAm impacts complex trait associations extensively
Although most mQTLs do not clearly impact human traits 33 , multiple studies have provided evidence that a small fraction of mQTLs are associated with human phenotypes and can point to causal disease-relevant pathways and contexts [18][19][20][21]33,57,58 . To evaluate the impact of mQTLs on traits in a systematic manner, and compare their effects to those of eQTLs, we performed QTL-GWAS colocalization by integrating FDR < 0.05 QTLs with 83 GWAS datasets that had at least one QTL-overlapping GWAS hit (P < 5 × 10 −8 ). To account for decreased mapping power but larger genome-wide abundance of mQTLs compared to eQTLs 44 , we estimated priors empirically and used a robust multi-method approach (Extended Data Fig. 8).
Across all GWASs, tissues and QTL types, we identified a total of N = 12,922 significant (regional colocalization probability > 0.3 and PP4 > 0.3) QTL-GWAS colocalizations (named simply 'colocalizations'; Supplementary Table 5). We observed that mQTL colocalizations were more abundant than eQTL colocalizations for almost all (91%) of GWASs (Fig. 4). Alike observations from other studies 20,21,57 , the observed overlap between eQTL-and mQTL-GWAS colocalizations is moderate, with 27% (749/2,734) of GWAS hits colocalizing with both QTL types in the same tissue (e/mQTL-shared), 55% of hits colocalizing with at least one mQTL but with no eQTLs (mQTL-specific) and 18% of hits colocalizing with at least one eQTL but with no mQTLs (eQTL-specific). The larger fraction of mQTL-specific colocalizations is consistent across trait groups (Extended Data Fig. 9a). These results highlight the importance of integrating different types of -omics data from different tissue sources, and considering secondary QTL signals, to maximize the possibility of identifying molecular links to inheritable traits.
For 19% (144/749) of the colocalized e/mQTL-GWAS shared loci, the mQTL-GWAS association shows greater colocalization probability than the eQTL-GWAS association in at least one tissue and/or independent QTL colocalization. We observe multiple instances where the lead e/mQTL variants are in low/moderate linkage disequilibrium and the mQTL mirrors the GWAS association substantially more optimally than eQTL does, as in the hypertension-associated MYO9B locus (Fig. 5d). We also observe cases where the GWAS locus harbors association signals compatible with the existence of multiple independent potentially causal variants, where GWAS colocalization with eQTLs and mQTLs contributes to differentiate and characterize these independent variants by their distinct colocalization patterns, as in the EFEMP2 locus linked to asthma (Fig. 5e). Together, these results provide evidence that integration of e/mQTL-GWAS colocalization signals can aid the fine-mapping of trait-associated variants and better characterize the molecular mechanisms underlying complex traits, as shown 2,61,62 .

Trait-linked mQTLs exhibit molecular regulatory pleiotropy and enrichment in trait-relevant tissues
Identification of QTL associations with traits in relevant tissues can provide insight into their underlying genetic and molecular mechanisms 6 . By analyzing the observed proportions of mQTL-GWAS colocalizations per tissue, we identified 18/65 traits with a disproportionate (test of equal proportions FDR < 0.05) amount of colocalizations in at least one tissue (Fig. 6a). Overall, the tissue with the largest proportion of colocalizations per trait matched the prior given current biological knowledge. For instance, blood clot and cell count traits were enriched in colocalizations derived from whole blood mQTLs, and breast cancer was nominally (Fisher's exact test P = 0.02) enriched in breast-derived mQTLs. For traits where the observed tissue link is less obvious, the observed enrichment can be spurious or it could point to an uncharacterized role of a specific tissue in the trait's biology. Together, these results suggest that mQTLs are potentially informative of complex traits' relevant tissues and can thereby aid the characterization of trait etiology, as observed for eQTLs [5][6][7] . However, this scenario should be   Article https://doi.org/10.1038/s41588-022-01248-z reevaluated with a mQTL catalog based on larger sample sizes and covering a wider range of potentially causal tissues. It is expected that many QTLs that impact traits do so by exhibiting molecular regulatory pleiotropy (that is, altering multiple molecules and/or molecular phenotypes). It has been shown that eQTLs where the lead eQTL variant (eVariant) regulates multiple eQTL genes (eGenes) are more likely to yield a trait association than eQTLs that regulate a single eGene 4 . However, the effect that regulatory pleiotropy plays in mQTL-trait associations has not been extensively characterized. Here, we observe that mQTL-GWAS colocalizations are enriched in mVariants regulating multiple, as opposed to single, mCpGs (OR = 2.65, Fisher's exact test P < 2.2 × 10 −16 ). Among trait-linked mCpGs that colocalize with at least one eGene, we observed enrichment (OR = 1.40, P = 6.3 × 10 −8 ) for trait colocalizations involving multiple mCpGs and eGenes (tier 4 in Extended Data Fig. 7a), as in HOXD locus (Fig. 3a). These results suggest that mQTLs that exhibit regulatory pleiotropy have increased probability of impacting a complex trait.

Trait-linked mQTLs are preferentially located in regulatory regions and exhibit decreased methylation
To better understand the DNAm features that contribute to mQTL impact on traits, we characterized DNAm levels of trait-linked mCpGs and their overlap with open chromatin regions. Compared to mCpGs without identifiable trait links, we observe an enrichment (Fisher's exact test P < 2.2 × 10 −16 , OR = 1.50) of trait-linked mCpGs in open chromatin regions; with which 53% (1,791/3,381) of trait-linked mCpGs overlap, with e/mQTL-shared GWAS-colocalized loci exhibiting a stronger enrichment (OR = 1.96) compared to mQTL-specific loci (OR = 1.36). Trait-linked mCpGs tend to have lower DNAm levels compared to trait-agnostic mCpGs (Wilcoxon rank-sum test P < 0.05), whereas trait-linked eGenes tend to be highly expressed (Extended Data  Fig. 9b). These results suggest that genetically regulated DNAm loci that play a role in trait etiology tend to correspond to lowly-methylated CpG sites in active chromatin regions.
Typically, variants contributing to the genetic basis of a trait are thought to act by affecting gene regulation. Concordantly, we observe an enrichment (Fisher's exact test P < 0.05) of trait-linked mCpGs in gene regulatory elements (OR = 1.63, P < 2.2 × 10 −16 ), compared to mCpGs without an identified trait link; 71% (2,390/3,381) of trait-linked mCpGs fall into this category. However, mCpGs corresponding to mQTL-specific colocalizations show a different profile compared to e/mQTL-shared colocalizations (Fig. 6b). Although e/mQTL-

Integration of trait-linked mQTLs with functional annotation enables identification of candidate causal genes
We detect a considerable amount of GWAS hits colocalizing with at least one mQTL but with no eQTLs (Fig. 4). To identify candidate causal genes for such loci, we integrated mQTL-specific trait-linked mCpGs with functional maps, that is curated promoter-and enhancer-gene target predictions 63 Table 7). Among highly supported loci (≥3 mCpG-gene links), we identify well-known gene-trait relationships, including APOB with cholesterol and ABO with blood traits. We also identified gene-trait links that are less characterized, such as RUNX1 with asthma and TMEM72 with red blood cell counts, with biological evidence compatible with predicted trait links (Supplementary Note). This strategy enabled the identification of candidate causal gene links for 1,129 GWAS loci (Fig. 6c). Overall, combining mQTL with functional To further identify trait-gene links, we jointly modeled GWAS, mQTL and eQTL signal across 61 multi-context eQTL maps derived from tissue and cell types across contexts not profiled in GTEx eQTL maps, including induced pluripotent stem cells and stimulated immune cells (Methods). We identified e/mQTL-GWAS colocalized clusters (PR > 0.8, PP > 0.5) for 56% (824/1,461) of loci previously identified as mQTL-specific when only considering GTEx matching-tissue eQTLs ( Fig. 6c and Supplementary Table 7). Clusters that lacked e/mQTLs colocalizations in the same tissue type tend to be present in relatively few eQTL contexts (Fig. 6d), potentially reflecting more spatiotemporally-specific trait-impacting regulatory events. Among identified context-specific colocalized eGenes (Supplementary Note), CLPTM1L (in adipose tissue) and TERT (in induced pluripotent stem cells) eQTLs colocalize with different GWAS signals for breast cancer subtypes (Extended Data Fig. 10), compatible with a putative causal TERT expression link hypothesized for this locus and trait (Fig. 5a). Our results illustrate the value of the mQTL-trait maps generated herein, integrated with additional functional genomic and multi-context eQTLs maps, to facilitate the identification of trait-linked candidate genes.

Discussion
Most prior studies that have characterized DNAm in relation to gene expression patterns have focused on blood cells 25,40 , with some exceptions 20 . In contrast, our study has characterized the genome-wide DNAm profile of diverse, healthy human tissue types in a systematic manner. In line with prior studies 20,65 , we observe that the aggregated contribution of mQTLs to human traits, in terms of number of identifiable associations, is larger than eQTLs, despite mQTLs being derived from substantially smaller sample sets. Consequently, we demonstrate that mQTLs can reveal a substantial number of molecular links to traits otherwise missed by eQTL-GWAS colocalization approaches. For these cases, mQTLs provide evidence of regulatory mechanisms underlying GWAS findings in absence of eQTL-based links to specific genes, and can pinpoint putative candidate genes. Our results indicate that in addition to expression and protein QTLs 66 , systematic integration of trait-linked mQTL associations with functional genomic maps can help guide the design of variant-to-function studies for complex traits.  (Fig. 4). Enrichment values correspond to maximum-likelihood estimates of the OR, derived from N = 270,783 trait-association tested mCpGs, and whiskers represent the 95% confidence interval of the enrichment value. enh., enhancer. c, Venn diagram represents overlap between GWAS hits with absence or presence of identified gene links by integrative approach. Analyzed GWAS hits were derived from the previous pairwise colocalization approach, considering mQTL-specific loci, that is mQTL-GWAS colocalizations in absence of eQTL-GWAS colocalization in the same GTEx tissue source (Fig. 4). d, Number of eQTL contexts/catalogs per eQTL-mQTL-GWAS-colocalized cluster (y axis, square-root transformed) derived from the multivariate colocalization approach considering 61 eQTL contexts/catalogs beyond GTEx (Methods). Observations are stratified by QTL-GWAS colocalization group derived from the previous pairwise colocalization approach, that is mQTL-GWAS colocalizations in presence (e/mQTL shared) or absence (mQTL specific) of eQTL-GWAS colocalization in the same GTEx tissue source (Fig. 4).
Article https://doi.org/10.1038/s41588-022-01248-z The lack of observed eQTLs for mQTL-specific GWAS colocalizations raises questions about how genetically regulated DNAm impacts trait-related biology. It is possible that DNAm-phenotype links involve additional molecular phenotypes other than gene expression 67 . However, our results suggest that eQTLs explain a substantial fraction of the missing links. That is, genes with low expression levels in bulk tissue or weak eQTL signal, as well as context-specific eGenes, are only revealed as eGenes for DNAm-trait linked loci when considering eQTL maps across multiple, distinct contexts (other than bulk tissue samples from adult nondiseased donors). Hence, we conclude that generally, trait-associated variants are more likely to result in detectable changes to DNAm than gene expression. In that sense, we highlight the usefulness of mQTL maps to proxy 'elusive' disease-associated eQTLs 68 , like in 69 . Taken together our results emphasize, supported by the trait-linked regulatory pleiotropy patterns observed, the importance of integrating multiple -omics, and the particular relevance of multi-tissue mQTL maps, to pinpoint molecular links and candidate genes to traits.
Importantly, maps of eQTMs, mQTLs and mQTL-GWAS links provided here are not comprehensive due to limited power, biased (gene-centric) set of CpG sites surveyed and lack of cellular resolution in these complex tissues (Experimental design limitations, Supplementary Note). The active versus passive role of DNAm variant-trait links 32 , and vertical versus horizontal pleiotropy scenarios 70 , are not resolved. In future studies, whole-genome bisulfite sequencing 44 and single-cell DNAm profiling 71 can contribute to address these limitations, and the active role of DNAm on corresponding gene targets can be interrogated with recently developed CRISPR-derived experimental approaches 72 . Larger GTEx DNAm cohort sample sizes would enable tissue-specific trans-mQTL mapping and further contribute to linking mCpGs to phenotypic variation, as observed in blood trans-mQTL studies 34,73 .
Altogether, the dataset generated herein constitutes the largest cohort with multi-tissue DNAm data generated to date. It allows for integrative -omics analyses by enhancing existing transcriptome-, epigenome-and proteome-based GTEx datasets 4,39,74,75 and provides the research community with a valuable resource to investigate the inherited susceptibility to human disease and complex traits from both cross-tissue and multi-omics perspectives. Our results contribute to a better understanding of the human DNA methylome and its relationship with both the transcriptome and complex traits.

Online content
Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41588-022-01248-z. Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.

Sample collection and ethical approval
The GTEx research protocol was reviewed by Chesapeake Bay Review, Roswell Park Comprehensive Cancer Center's Office of Research Subject Protection and the institutional review boards at the University of Pennsylvania. Informed consent, standards and protocols for optimizing postmortem tissue collection and donor recruitment are detailed elsewhere 76 . Families of participant donors were not compensated. Analyses of GTEx DNA samples at the University of Chicago were not considered human subjects research by the institutional review board at the University of Chicago, because only deidentified data on deceased individuals were involved.

Statistics and reproducibility
No statistical method was used to predetermine sample size, as for each DNAm profiled tissue, sample size was chosen based on availability of biosamples. Samples were randomly allocated in plates; tissue types were not batched by plate. The experiments were not randomized. The investigators were not blinded to allocation during the experiments and outcome assessment, because the data are not from controlled randomized studies. We used all data passing standard quality controls (QCs), resulting in 987 samples.

Generation of epigenome-wide DNAm data
Epigenome-wide DNAm was derived from 1,000 tissue samples from 9 unique tissue types obtained from 424 GTEx subjects, mostly (~83%) from European ancestry given available genotype data 4 (Tables S4 and  S5), where variants were defined at minor allele frequency (MAF) > 5% in 1000 Genomes Project. A total of 754,054 CpGs passed QC, could be mapped and lifted over from GRCh37 to GRCh38 human genome build and were retained for further analysis.
Raw DNAm values were background-adjusted using the single sample normal-exponential out-of-band (ssnoob) method 80,81 with dye bias correction implemented in minfi (v.1.36.0). DNAm β values were normalized using the beta mixture quantile (BMIQ) method, implemented in ChAMP v.2.8.6, adjusting for type I/II probe bias 82 ; output β values for each CpG were used in downstream analysis. After normalization, we removed one additional sample with an array-derived genotype profile not matching the whole-genome sequencing (WGS)-derived one. Principal-component analysis was conducted on DNAm β values within each tissue type, and three samples were removed for being outliers with respect to the top five principal components of the corresponding tissue. A total of 987 samples passed final QC and were retained for further analysis.

Estimation of tissue similarity based on DNAm and gene expression
Mean DNAm and gene expression levels for each tissue type were used to calculate tissue distances (1 − Spearman rank-correlation coefficient (Spearman's rho)) in a pairwise fashion. For DNAm, Spearman's rho was calculated on the mean DNAm for each CpG (across all samples for a given tissue) in M-value units, and each CpG was weighted by its cross-tissue variance. The entire set of post-QC profiled CpGs (N = 754,054) was considered for the analysis. For gene expression, Spearman's rho was calculated on the mean expression for each gene (across all samples for a given tissue) in log 2 (x + 1) transformed transcript per million units, and each gene was weighted by its cross-tissue variance. Only genes expressed in at least one of the nine tissues (N = 37,686) were considered for the analysis. We performed hierarchical clustering of the tissues using pvclust 2.0.0 (ref. 83 ) with complete linkage and nboot = 1,000.

eQTM mapping
We define an eQTM as the association of CpG DNAm with proximal gene expression, considering a ±1.5 Mb window centered on the CpG locus, thereby enabling the identification of short-and long-range CpG-gene associations. We analyzed samples with available DNAm, expression and genotype data (Supplementary Table 1), comprising a total of 490 samples, from 25 (testis) to 131 (lung) per tissue, excluding kidney cortex due to insufficient sample size (N = 5). We obtained DNAm residuals by regressing DNAm-derived probabilistic estimation of expression residuals (PEER) factors 84 used in cis-mQTL mapping (see below) from inverse-normalized-transformed DNAm levels, and gene expression residuals by regressing expression-derived PEER factors used for eQTL mapping in 4 from inverse-normalized-transformed gene expression levels. For each CpG site in each tissue, we calculated Spearman correlation of DNAm with gene expression residuals of proximal (±1.5 Mb window) genes. We applied Bonferroni multiple testing correction to nominal P values, accounting for multiple genes tested per CpG. Then, we adjusted for multiple CpGs tested by applying q-value multiple testing correction 85 to the set of top Bonferroni-adjusted P values per CpG. We defined significant CpGs involved in eQTMs (eCpGs) at FDR < 0.05, and for those, significant eQTMs were defined at Bonferroni-adjusted P value < 0.05. To overcome QTM-mapping limited power due to per-tissue available sample sizes, we used an approach to perform a cross-tissue QTL analysis by leveraging QTL signal across tissues 42 , implemented in the R package mashr (v.0.2.6). We assessed eQTM replication for all tissues in the FUSION Skeletal Muscle Study cohort 20 . Replication was assessed by means of π1, which enables the estimation of the true positive rate of findings derived from a discovery dataset in a replication dataset 86 . See additional details of eQTM identification and characterization in Supplementary Note.

QTL mapping for methylation and gene expression QTLs (mQTLs, eQTLs)
We define mQTLs as proximal variants (that is, in cis) to a CpG with a significant genotype effect on its DNAm estimates, considering a ± 500 kb window from the CpG locus. To assess mQTLs, we considered quality controlled (QC-ed) inverse-normalized DNAm data, generated and presented here as part of the eGTEx project, and QC-ed genotype data derived from GTEx v8 4 filtered at variant MAF > 0.01 per tissue. For each variant-CpG pair, we used an adaptation of v.2.184 FastQTL 87 , available at https://github.com/broadinstitute/ gtex-pipeline/tree/master/qtl, to fit a linear regression model (core_ PEERs) separately in each tissue, and tested for significance of genotype on methylation estimates while adjusting for additional known and unknown factors: where Y is the inverse-normal-transformed DNAm levels, β 0 is the intercept and β is the corresponding effect size. β G is the effect size of genotype on DNAm.
C represents a subset of covariates that were used in cis-eQTL mapping 4 . These covariates include five genotype principal components, two covariates derived from the generation of genotype data by WGS and biological sex status. The WGS covariates are described elsewhere 4 and represent the WGS sequencing platform (HiSeq 2000 or HiSeq X) and WGS library construction protocol (PCR based or PCR-free).
PEER represents PEER factors 84 derived from DNAm. See additional details of PEER generation in Supplementary Note.
We corrected for multiple testing of variants per CpG 87 and multiple CpGs tested 86 , defining significant mQTL CpGs (mCpGs) at FDR < 0.05. A similar approach was used to identify eQTLs. Conditional QTL analysis was used to identify multiple independent mQTLs and eQTLs for mCpGs and eGenes, respectively, as well as corresponding lead variants for each independent QTL locus. This approach accounts for allelic heterogeneity, where distinct genetic variants at a locus simultaneously and independently affect methylation at a given CpG site. For that, we applied a stepwise regression procedure as described elsewhere 43 (details can be found in Supplementary Note). Similarly as in eQTM mapping, we used mashr to identify mQTLs and eQTLs detectable after leveraging QTL signal across tissues (Supplementary Note). Across the article, we refer to CpGs and genes with at least one significant mQTL or eQTL as mCpGs and eGenes, respectively, and to mQTL and eQTL significant variants as mVariants and eVariants, respectively. We assessed mQTL replication and/or validation of mQTL tissue-specific patterns for all tissues in the FUSION Skeletal Muscle Study 20 (https://www.ebi.ac.uk/birney-srv/FUSION/), the ROSMAP brain 23 (http://mostafavilab.stat.ubc.ca/xQTLServe/) and the GoDMC studies 33 (http://mqtldb.godmc.org.uk/) (Supplementary Note, Supplementary Table 3).
To evaluate the adequateness of the mQTL mapping model of choice, two alternative mQTL mapping models were benchmarked: For each model, we assessed the number of significant (FDR < 0.05) mCpGs, the proportion of true positives, and the replication rate in the FUSION cohort. We observe that the chosen model (core_PEERs) outperforms alternative model choices, as it detects on average across tissues 15-230% more mCpGs than the other models tested, with similar replication rates (Supplementary Table 3). See additional details of mQTL and eQTL identification, characterization and validation in Supplementary Note.

QTL enrichment in genomic annotations
Functional enrichment analyses were performed using torus 88 (v.1.0.0.dev), similarly to GTEx Consortium 4 . In brief, the command 'torus -d ${qtl_statistics} -annot ${annotation_file} -est-fastqtl' was used, where ${qtl_statistics} correspond to QTL-mapping data for eQTLs (full eQTL-tested gene set per tissue) or mQTLs (subset of 21,000 mQTL-tested CpGs), and ${annotation_file} corresponds to QTL-tested annotated variants. Variants were annotated with Ensembl's Variant Effect Predictor using datasets from several sources: gene regulatory element and open chromatin annotations were obtained from ENCODE5, gene body annotations were obtained from Ensembl, chromatin state predictions from ROADMAP, TF binding and CGI annotations were obtained from UCSC browser (see Obtention and processing of functional genomic data). For each tissue, torus-derived enrichment Article https://doi.org/10.1038/s41588-022-01248-z estimates correspond to point estimates (derived by maximum likelihood estimation) of the log of odds ratio. Enrichment across tissues was evaluated by modeling single-tissue enrichment estimates (log of odds ratio) with a random-effects model (rma function, metafor R package v.2.0.0). Single-tissue and cross-tissue enrichment estimates are provided in Supplementary Table 4.

Colocalization of mQTLs with eQTLs
We investigated the associations between mQTLs and eQTLs by means of QTL effect size colocalization with coloc 89 (v4.0.0) using default priors. For each significantly (PP4 > 0.5) colocalized mQTL-eQTL pair in each tissue, the top-colocalized e/mVariant was defined as the one with the largest PP4 value. See additional details of mQTL-eQTL colocalization in Supplementary Note.

Characterization of colocalized mQTL-eQTL loci
Identification of mQTL-eQTL regulatory pleiotropy. Considering QTLs, we define regulatory pleiotropy as the event of a variant or a set of variants in a QTL region impacting multiple eGenes and/or mCpGs. To characterize the pleiotropic nature of mQTLs, we quantified the number of mCpGs and eGenes involved in loci harboring eQTL-mQTL (e/mQTL) colocalizations, and classified mCpGs involved in at least one significant e/mQTL colocalization (PP4 > 0.50) in four tiers, depending on their mCpG and eGene connectivity level (Extended Data Fig. 6). Of note, no inference of the nature of the pleiotropic effect, whether vertical or horizontal, is made in this classification. Tiers are defined as follows: (a) tier 1: mCpGs that colocalize with a single eGene and vice versa (1:1 connectivity); b) tier 2: mCpGs that colocalize with multiple eGenes, and each one of those eGenes uniquely colocalizes with that single mCpG (1:m connectivity); (c) tier 3: mCpGs that colocalize with a single eGene, which colocalizes with multiple mCpGs (n:1 connectivity); and (d) tier 4: mCpGs that colocalize with a multiple eGenes, where at least one of the eGenes colocalize with multiple mCpGs (n:m connectivity).

Colocalization of QTL with GWAS signal
Colocalization of ovary cancer GWAS with QTL signal of HOXD pleiotropic locus. We hypothesized that the mCpGs and eGenes in the HOXD region may be linked to ovarian cancer risk and tested it by means of QTL-GWAS colocalization. Ovary cancer GWAS summary statistics 90 were obtained from the Ovarian Cancer Association Consortium (OCAC) website http://ocac.ccge.medschl.cam.ac.uk/ and were filtered (not imputed) as in Barbeira et al. 5 . Considering OCAC GWAS along with QTL statistics from the set of pleiotropic mCpGs and eGenes identified in the HOXD locus (Fig. 3a), we performed colocalization analysis with coloc 89 (v4.0.0) using default priors. Considering QTLs statistics, colocalization was performed based on effect size and associated standard error values; P values and corresponding variant MAFs were used for OCAC GWAS data. The probability of one causal variant associated with both traits (PP4) was used to identify significant (PP4 > 0.50) colocalizations.
Determination of GWAS significant loci. To investigate possible associations between genetically regulated molecular and complex traits, including disease and 'healthy' phenotypes, we used GWAS summary statistics of 87 GWASs; data production and QC are described in detail in Barbeira et al. 5 . We identify significant GWAS hit loci (that is, genomic windows containing GWAS signal) similarly as described in Barbeira et al. 5 . In brief, the GWAS summary statistics were split into 1,702 approximately linkage disequilibrium-independent regions 91 (Supplementary Table 5). Each region was categorized as a significant GWAS hit locus (GWAS hit) provided it encompassed a non-imputed GWAS significant (P < 5 × 10 −8 ) variant.
Determination of QTL-GWAS significant loci. For each GWAS and each QTL tissue pair, colocalization was performed with coloc (v.4.0.0) and fastenloc (v.1.0). The latter is an improved version of enloc 92 and was recently described in multiple studies 93,94 . coloc approach. For each GWAS, at each GWAS locus, we identified overlapping (>1 bp) mCpG loci from each of the 9 analyzed tissues, considering per-tissue significant (FDR < 0.05) mCpGs resulting from the single-tissue QTL-mapping approach. For each overlapping mCpG-GWAS region pair, we applied coloc 89 to mQTL along with GWAS summary statistics. For mCpGs with secondary QTLs, that is multiple independent mQTLs, conditional -to independent lead mVariants -QTL signals were also tested for colocalization. By this, we prevent putative colocalizations to be missed or miscalculated, as coloc assumes a single variant to be causal of QTL-GWAS effects 94 . Prior probabilities of a variant yielding (a) a mQTL association (p2), (b) a GWAS association (p1) and (c) a mQTL and a GWAS association (p12) were estimated from fastenloc enrichment values in an analogous manner as done with enloc in 5 and provided in Supplementary Table 6. Only the regions with at least 50 variants in common between the GWAS and mCpG loci were tested for colocalization. Both for QTLs and GWAS statistics, colocalization was performed on effect size (effect size) and associated standard error (effect size s.e.) values. Used GWAS statistics were imputed from available z-score (z), allele frequency (f) and sample size values (N) by effect size ≈ z/(f(1 − f)N)^(1/2) and effect size s.e. ≈ effect size/z. fastenloc approach. First, GWAS posterior inclusion probability values were obtained with torus from available z-scores. Then, for each tissue, the DNAm levels, genotypes and mQTL-mapping covariates for CpGs and corresponding cis windows considered in the mQTL analysis were processed with dap-g 95 (v.1.0.0). Next, we used fastenloc to obtain regional colocalization probabilities for all tuples of interest (GWAS hit, trait, tissue, mCpG) by subsetting corresponding GWAS hits, and significant (single-tissue mQTL set: FDR < 0.05) mCpGs from the genome-wide fastenloc output. An equivalent approach was used for eQTL-GWAS colocalization; the dap-g precomputed eQTL annotations from GTEx (v8) data were obtained from the fastenloc github repository: https://github.com/xqwen/fastenloc. Evaluation of mQTL-GWAS colocalization approach is described in Supplementary Note.

Characterization of signatures of trait-linked mCpGs
Tissue-specific enrichment in mQTL-GWAS colocalizations. To identify traits with a tissue-specific colocalization enrichment profile, we performed a multisample test for equality of proportions without continuity correction. For each trait with >5 mQTL-GWAS colocalizations, we compared the observed proportions of mQTL-GWAS colocalizations per tissue to the overall proportion of colocalizations for all tissues. We identified traits with a disproportionate amount of colocalizations in at least one tissue at Bonferroni-adjusted P < 0.05. For these traits, scaled-by-trait colocalization proportions per tissue are shown in Fig. 6a, where tissues and traits were clustered via complete-linkage hierarchical clustering based on euclidean distance. A two-sided Fisher's exact test was conducted to determine significant tissue-trait enrichments at Bonferroni-adjusted P < 0.01, which are indicated by crossed cells in Fig. 6a. This approach has several limitations. We did not observe tissues of relevance for several traits, possibly due to the examination of mQTLs derived from only 9 tissues. For certain tissue-trait pairs, the limitedness of the mQTL catalog examined herein, and the low number of colocalizations identified, may have also impacted the completeness and accuracy of the identified tissue-specific enrichment profiles. This analysis would benefit from a more powered and exhaustive mQTL catalog.
Characterization of molecular regulatory pleiotropy of trait-linked mCpGs. We evaluated the enrichment of mVariants corresponding to mQTL-GWAS colocalizations in mVariants regulating multiple versus single mCpGs. For each mCpG tested for mQTL-GWAS colocalization, Article https://doi.org/10.1038/s41588-022-01248-z we kept the corresponding colocalization result with the highest coloc PP4 value and the corresponding mVariant. We classified mVariants as pleiotropic if they were estimated to be associated with multiple mCpGs, or as non-pleiotropic, if they were associated with a single mCpG. Association was defined considering dap-g fine-mapping estimates used in the fastenloc mQTL-GWAS colocalization approach, and mQTL credible sets were defined at 90% confidence. We classified mVariants as being also eVariants if they were estimated to be associated with at least one eGene, estimating eQTL credible sets in an analogous manner to mQTL sets. Significant enrichment of trait-linked versus trait-unlinked mVariants in multiple-mCpG and eVariant sets was defined at two-sided Fisher's exact test P < 0.05. Additionally, we subsetted mCpGs tested for mQTL-GWAS colocalization that were involved in e/mQTL colocalizations, and classified them into mCpG-eGene pleiotropy tiers (Extended Data Fig. 6). Enrichment of trait-linked versus trait-unlinked eGene and mCpGs in each tier was estimated at two-sided Fisher's exact test P < 0.05. Given its complex linkage disequilibrium structure and genotype imputation inaccuracies, the major histocompatibility complex (MHC) locus is challenging to analyze. To ensure that pleiotropy results are not biased by MHC, we evaluated pleiotropy excluding CpG cis loci overlapping the MHC region, as well as corresponding eQTL-mQTL colocalized eGenes. The pleiotropy distributions were not significantly different (test of equal proportions P < 0.05). We evaluated the putative enrichment of increased disease-risk alleles in gains or losses of DNAm, both across and within GWAS traits. From the set of 87 GWAS tested for mQTL colocalization, we subsetted colocalization results corresponding to 31 GWAS disease traits with for which at least one significant mQTL colocalization was identified. For each GWAS trait and mCpG tested for mQTL-GWAS colocalization, we kept the corresponding colocalization result with the highest coloc PP4 value. For each GWAS trait, we classified tested mCpGs into methylation gain or loss categories, based on the sign of the allelic mQTL effect associated with increased disease risk (that is, positive GWAS effect). Enrichment of each disease in methylation gains and losses was determined at Bonferroni-corrected two-sided Fisher's exact test P < 0.05. No disease was identified as enriched for DNAm gains or losses. Despite individual traits not being significantly enriched for DNAm change biases, it is possible that a global bias could be detected when considering all traits at once. To test this hypothesis, each disease was labeled as 'DNAm-gain' if the corresponding number of colocalized mCpGs associated with a DNAm gain outnumbered the DNAm-loss ones and labeled as 'DNAm-loss' otherwise. An exact binomial test was performed to assess whether the proportion of DNAm-gain or DNAm-loss labeled diseases deviated significantly (P < 0.05) from 50%, which was not the case (P = 0.86).

Characterization of
Characterization of trait-linked mCpGs overlap with gene regulatory regions. Gene regulatory element annotations were derived from ENCODE5 cCREs catalog and gene body element annotations were obtained from GENCODE (see Obtention and processing of functional genomic data). To annotate mCpGs for cCRE and gene body elements, we extended the span of their genomic location by ±100 bps and checked for overlap (≥1 bp) with element regions. Trait-linked mCpGs were classified as e/mQTL shared or mQTL specific (Supplementary Note: Evaluation of mQTL-GWAS colocalization). Enrichment significance of e/mQTL-shared or mQTL-specific trait-linked mCpGs in each gene regulatory region or gene body element category was estimated at two-sided Fisher's exact test P < 0.05.
Characterization of trait-linked mCpGs overlap with functional maps. To evaluate the overlap of mQTL-specific trait-linked mCpGs with genomic regions linked to specific genes, we considered functional maps based on curated regulatory-region to gene target predictions 63,64 and eQTM predictions generated herein. We extended the genomic location span of mCpGs by ±100 bp, and checked for overlap (≥1 bp) with eCpGs and regulatory regions (gene links). Enrichment of mQTL-specific trait-linked mCpG in gene links was evaluated with two-sided Fisher's exact test, considering gene links observed for non-colocalized mCpGs. For each trait/gene-linked mCpG, we annotated corresponding regulatory-region gene targets and eCpG-correlated genes; a particular gene was annotated to a mCpG if overlap was found in one or more functional maps. Annotations were collapsed at a GWAS hit level by adding up the number of predicted mCpG-gene links corresponding to the locus, hence providing an estimate of mCpGs that support a gene candidate per trait-associated locus (Supplementary Table 6).

Colocalization of trait-linked mCpGs with multi-context eQTL maps.
To expand on the identification of eGenes underlying mQTL-GWAS specific GWAS hits (Fig. 4), we considered non-GTEx multi-context eQTL resources and performed GWAS-mQTL-eQTL multivariate colocalization with HyPrColoc 96 (v.1.0.0), which implements an efficient deterministic Bayesian algorithm to detect colocalization across multiple traits (for example, multi-omic traits from different contexts). The eQTL mappings were obtained from eQTL Catalogue 37 (https:// www.ebi.ac.uk/eqtl/) release 4, eQTLGen 97 (https://eqtlgen.org/ cis-eqtls.html) and i2QTL 98 (https://doi.org/10.5281/zenodo.4005576) Consortium. The eQTL map set comprised 61 gene-level eQTL mappings identified in cohorts other than GTEx, most of derived from RNA sequencing, encompassing a wide range of tissues and cells, including simulated blood cells, highly powered pluripotent cells and whole blood datasets (Supplementary Table 8). For each mQTL-GWAS specific GWAS hit (Fig. 4), we jointly modeled (1) the mQTL signal corresponding to the top colocalizing mCpG, (2) the external eQTL datasets and eQTL signal for genes with nominal (P < 1 × 10 −3 ) eQTL signal for a variant matching the top colocalizing mQTL-GWAS variant and (2) the GWAS signal of the colocalized trait. Only cases with >50 overlapping variants were considered. A single external eQTL dataset was modeled at a time. We refer to PR as the regional association probability that all the traits share an association with one or more variants within the region and PA as the alignment probability that the shared associations between all traits are owing to a single shared putative causal variant. The posterior probability of full colocalization (PPFC = PR*PA) represents the posterior probability that all traits share a causal variant within that region. Clusters were identified when PR > 0.8 and PPFC > 0.5 (PA > 0.625) were satisfied. Prior.1 represents the prior probability that a variant is associated with a single trait, and (1 − prior.2) represents the prior probability that a variant is associated with an additional trait given that it is associated with one trait. We set prior.1 = 1 × 10 −4 , and prior.2 = 0.98. We performed sensitivity analysis by setting prior.2 to 0.95, 0.98 and 0.99, and we also calculated the corresponding empirical false positive rate by permuting eQTL with respect to mQTL and GWAS estimates (Supplementary Note); significant colocalization instances are provided (Supplementary Table 8).

Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Code availability
Code for QTL and eQTM mapping, functional enrichment, and colocalization, as well as code to to generate manuscript figures, is available at the github repository (https://github.com/meritxellop/ eGTEx_mQTLs_eQTLs_GWAS) and archived at zenodo (https://doi. org/10.5281/zenodo.7106660) 99 .  . (b-c) Cross-tissue sharing of mCpGs and eGenes. Crosstissue mean percent of mCpGs and eGenes per tissue-sharing category is shown in parentheses. Of note, testis is an outlier for eQTL tissue specificity, as 23.5% of eGenes were not detected in any other tissue. Avg.: average (mean). (c) Cross-tissue sharing of mQTL tissue-leveraged effect magnitudes (y axis) per gene regulatory region (x axis, 36 data points per box plot). P-values of paired two-sided Wilcoxon signed rank tests are shown for corresponding pairwise comparisons; p-value of Kruskal-Wallis rank sum test is shown for the three-way comparison. Enh.: Enhancer. (d, e) Validation of tissue-specific mQTLs in muscle, blood and brain mQTL external cohorts, see (b) for tissue color legend. In (d), for each cohort, the average of absolute mQTL effect sizes and corresponding standard error is displayed for each set of tissue-specific mQTLs identified in each GTEx tissue (x axis). In (e), Spearman correlation between external and GTEx mQTL effect sizes, and associated standard error, is shown. Fisher's z transformation was applied to Spearman correlation coefficients; standard errors were calculated based on the transformed coefficients. In (d,e), the number of QTL associations (N) tested for each pairwise comparison is as follows: N FUSION,GTEx = 195|4142|336,4630,287,4643|1428|360|369, N GoDMC,GTEx = 22|1353|47|1478|70|1019|292|83|102 and N ROSMAP,GTEx = 57|2428|156|2587|163|2673|791|214|109, where GTEx corresponds to Breast|Colon|Kidney|Lung|Muscle|Ovary|Prostate|Testis|Blood tissues, respectively.