Joint mouse-human phenome-wide association to test gene function and disease risk.

Overview of phenome data for the BXD cohort.(a) Five pairs of isogenic BXD cohort strains—BXD43 to BXD102. There are now approximately 100 readily available BXD strains and another 50 that are almost fully inbred. Almost all current phenome data is restricted to the parents, F1 hybrids (B6D2F1 and D2B6F1) and BXD1 through BXD102. (b) Phenome data categorized by type, including classic phenotypes (top and see comprehensive definitions in Supplementary Data 1), metabolic and proteomic trait data (middle), and independent mRNA expression assays (bottom, n=86 unique eQTL data sets, see Supplementary Note 1 Summary Expression Phenome). (c) Body weight data for BXD strains on high-fat (grey) and low-fat (black) diets taken from Wu et al.10. (d) Expression of Bckdhb mRNA and its protein in six tissues for the five BXD strains. Protein data from Wu et al.10

Overview of experimental phenome-wide association.(a) Workflow for the PheWAS analysis. D2 genome was sequenced using Illumina and SOLiD platforms and aligned to the reference genome. Approximately five million sequence and structural variants were identified. A set of high-impact variants including 11,979 missense, 58 nonsense, 70 splice site, 45 frame shift and 276 CNVs spanning at least one gene was used for classic and molecular PheWAS analyses. (b) Whole-genome read coverage. The distribution was generated using combined mapped reads from Illumina and SOLiD, and mapped reads from Illumina and SOLiD platforms individually. (c) Circos plot of sequence and structure variants segregating in the BXD cohort. The outmost circle represents SNP density per 100 kb window (black at the lowest density and orange at the highest density). The second circle represents indel density per 100 kb. The innermost circle represents CNVs. Blue and green ticks indicate D2 losses and gains, respectively. (d) Venn diagram of missense variants predicted to be deleterious by SIFT and PolyPhen2. (e) False positive rate of strong variants. A small subset of each type of variants was selected for validation by Sanger resequencing.

Association analysis for a missense variant in Fh1.(a) Structure of the Fh1 gene showing lyase 1 and fumarase c domains, the former of which contains a missense mutation. (b) Combined eQTL mapping of Fh1 mRNA across 15 tissues. The eigenvalues associated with the first principal component map to Fh1 with a likelihood ratio statistic of >20. The solid red line represents genome-wide significance. The red triangle indicates the genomic position of Fh1. (c) Manhattan plot of an expression phenome scan of molecular traits linked to the Fh1 locus in midbrain. The y-axis shows the −log10 q values of association, and the x-axis shows positions of 55,681 assays generated using Agilent SurePrint array. (d) QTL heat map of mRNAs involved in the unfolded protein response. The x-axis lists mouse chromosome numbers. Each horizontal line represents the QTL map for a single transcript in midbrain. Transcripts are grouped into three major categories—genes involved in the canonical UPRmt, cytoplasmic heat shock response (HSR) and unfolded protein response in the endoplasmic reticulum (UPRer). A subset of UPRmt genes at the top are strongly modulated by the Fh1 locus on Chr 1 (the intense colours to the upper left). In contrast, none of the UPRer subsets are modulated by Fh1. (e) Principal component analysis plot for six UPRmt transcripts (left). The first two components explain ∼67% of the variance in expression in hypothalamus. There is significant correlation between UPRmt expression and Fh1 (P=5 × 10−8; Pearson product-moment correlation coefficient) (right). (f) Validation that fumarate hydratase selectively controls the UPRmt in C. elegans. The left-most pair of images demonstrates effects of the fum-1 RNAi knockdown on hsp-6::gfp signal—a marker of UPRmt induction. The middle and right panels demonstrate that the fum-1 knockdown does not induce either the UPRer (hsp-4::gfp), or the cytoplasmic heat shock response (hsp-16.2::gfp). (g) Manhattan plot of a phenome scan. The y-axis shows the −log10 q values of ∼4,230 phenotypes, and the x-axis shows the 15 phenotypic categories.

Association analysis for missense variants in Col6a5.(a) Twenty missense variants in Col6a5 distributed across 10 von Willebrand factor A-type (vWFA) domains. (b) Differential mRNA expression of Col6a5 in tibias (n=4) measured by rtPCR. The D haplotype (blue, right) has far higher expression than the B haplotype (green) relative to Gapdh. (c) Phenome scan of Col6a5 (rs13480398) across mRNA assays for femur. (d) Phenome scan of Col6a5 (rs13480398) across classic phenotypes. (e) Marked difference in bone density between B6 and D2 parents. Femurs from 12-week-old mice were scanned using high-resolution micro-CT (μCT40, SCANCO Medical, Bassersdorf, Switzerland). More highly mineralized areas are indicated in red. (f) Difference in material bone density (P=0.02; two-tailed Student's t-test, n=3). (g) Human phenome scan for association of Col6a5 (rs113396273) across BioVU.

Association analysis for nonsense variant in Ahr.(a) Structure of the Ahr gene showing three domains with a nonsense mutation and five missense mutations. The nonsense mutation (*805R) leads to loss of the stop codon, and the addition of 43 C-terminal amino acids. Dotted rectangle to the right is the extended coding region in the D haplotype. (b) Expression phenome scan of Ahr (rs3711448) across liver mRNA levels in the BXDs. (c) Phenome scan of Ahr across classic phenotypes in the BXD strains. Both AHR protein level and cleft palate induced by TCDF injection are strongly linked to Ahr. (d) Manhattan plot showing the association in human between SNP (rs2066853) in AHR and classic phenotypes. The cleft palate phenotype is also associate with AHR in human clinical cohorts.

Association analysis for CNV covering Alad and Hdhd3.(a) The CNV region for Alad and Hdhd3 derived using read-depth information from genome sequencing. Red dots represent at least a two-fold increase in coverage compared with the reference genome. The x-axis shows the reference genomic position of the CNV. Two gene models (that is, Hdhd3 and Alad) are shown in the CNV plot. (b,c) Rank ordered mean expression levels of Hdhd3 and Alad across 67 BXD strains, their parental strains, and F1 crosses. Expression values are normalized on a log2 scale (mean±s.e.m.). Strains with D alleles (red) have higher levels of Alad and Hdhd3 compared with B alleles (green). F1 hybrids (blue) are intermediate. The comparison between B and D alleles for Alad and Hdhd3 are shown in an inset boxplot. (d) The phenome scan of the BXD cohort highlights several interesting potential phenotypes including pain response (thermal nociception), brain deoxycorticosterone levels, and antigenic activity in the spleen. Two triangles represent pigmentation traits that we know are associated with a variant in the linkage disequilibrium block. (e) Manhattan plot obtained after phenome scan of the BioVU EHR data showing the association in humans between a SNP (rs1800435) in ALAD and chronic pain syndrome and several other classic phenotypes.

Compensation for Actc1 loss by upregulation of Acta1 and Acta2.(a) Structure of the cardiac actin gene, Actc1, and the tandem duplication (dotted box) in strains with the D haplotype. (b) Actc1 expression in heart is controlled by a strong cis-eQTL that corresponds precisely to the location of the Actc1 gene on Chr 2 at 114 Mb (solid red triangle). X-axis represents the megabase coordinate on Chr 2 while the y-axis represents the LOD linkage score. Thin horizontal lines provide genome-wide significance thresholds (upper red line at P<0.05 and lower grey line at P<0.63). (c) Mean expression of Actc1, Acta1 and Acta2 highlights the partial compensation for B haplotypes (green, n=24) and D haplotypes (blue, n=12). Expression values are normalized on a log2 scale (mean±s.e.m.). Heart expression data is from GeneNetwork data set GN485 (EPFL/LISP BXD CD-HFD Heart Affy Mouse Gene 2.0 ST (Jan14) RMA). (d,e) Both skeletal muscle and smooth muscle actins (Acta1 and Acta2) also map precisely to Actc1.