Leveraging inter-individual transcriptional correlation structure to infer discrete signaling mechanisms across metabolic tissues

Abstract

Inter-organ communication is a vital process to maintain physiologic homeostasis, and its dysregulation contributes to many human diseases. Given that circulating bioactive factors are stable in serum, occur naturally, and are easily assayed from blood, they present obvious focal molecules for therapeutic intervention and biomarker development. Recently, studies have shown that secreted proteins mediating inter-tissue signaling could be identified by 'brute force' surveys of all genes within RNA-sequencing measures across tissues within a population. Expanding on this intuition, we reasoned that parallel strategies could be used to understand how individual genes mediate signaling across metabolic tissues through correlative analyses of gene variation between individuals. Thus, comparison of quantitative levels of gene expression relationships between organs in a population could aid in understanding cross-organ signaling. Here, we surveyed gene-gene correlation structure across 18 metabolic tissues in 310 human individuals and 7 tissues in 103 diverse strains of mice fed a normal chow or high-fat/high-sucrose (HFHS) diet. Variation of genes such as FGF21, ADIPOQ, GCG, and IL6 showed enrichments which recapitulate experimental observations. Further, similar analyses were applied to explore both within-tissue signaling mechanisms (liver PCSK9) and genes encoding enzymes producing metabolites (adipose PNPLA2), where inter-individual correlation structure aligned with known roles for these critical metabolic pathways. Examination of sex hormone receptor correlations in mice highlighted the difference of tissue-specific variation in relationships with metabolic traits. We refer to this resource as gene-derived correlations across tissues (GD-CAT) where all tools and data are built into a web portal enabling users to perform these analyses without a single line of code (gdcat.org). This resource enables querying of any gene in any tissue to find correlated patterns of genes, cell types, pathways, and network architectures across metabolic organs.

Keywords: endocrine; genetics; genomics; human; mouse; organ cross-talk; systems genetics.

Figure 1.. Web tool overview and inter-individual correlation structure of established endocrine proteins.

(A) Web server structure for user-defined interactions, as well as server and shiny app implementation scheme for gene-derived correlations across tissues (GD-CAT). (B) All genes across the 18 metabolic tissues in 310 individuals were correlated with expression of ADIPOQ in subcutaneous adipose tissue, where a q-value cutoff of q<0.1 showed the strongest enrichments with subcutaneous and muscle gene expression (pie chart, left). Gene set enrichment analysis (GSEA) was performed using the bicor coefficient of all genes to ADIPOQ using gene ontology biological process annotations and network construction of top pathways using clusterprofiler, where pathways related to fatty acid oxidation were observed in adipose (left) and chemotaxis/ECM remodeling in skeletal muscle (right). (B–D) The same q-value binning, top within-tissue and top peripheral enrichments were applied to intestinal GCG (C), liver FGF21 (D), and muscle IL6 (E). For these analyses all 310 individuals (across both sexes) were used and q-value adjustments calculated using a Benjamini-Hochberg FDR adjustment.

Figure 1—figure supplement 1.. Performance across four methods of cell-type deconvolution where relative proportions of cells (y-axis) are shown for all cell types annotated in single-cell reference (x-axis) in liver.

Figure 1—figure supplement 2.. Performance across four methods of cell-type deconvolution where relative proportions of cells (y-axis) are shown for all cell types annotated in single-cell reference (x-axis) in heart.

Figure 1—figure supplement 3.. Performance across four methods of cell-type deconvolution where relative proportions of cells (y-axis) are shown for all cell types annotated in single-cell reference (x-axis) in skeletal muscle.

Figure 1—figure supplement 4.. Pancreatic INS expression correlations across tissues in gene-by-tissue expression project (GTEx) were binned according to q<0.1 (top) and corresponding pancreatic gene set enrichment analysis (GSEA) network graph is shown (bottom).

Figure 2.. Tissue-specific contributions to pan-organ gene-gene correlation structure.

(A) Heatmap showing all the number of gene-gene correlations across tissues which achieve significance relative to total number of genes in each pathway at biweigth midcorrelation Student’s p-value <1e-3 (left column), p-value <1e-6 (left middle column) of BH-corrected q-value <0.1 (right middle column) or BH-corrected q-value <0.01 (right column). Within-tissue correlations are omitted from this analysis. (B–D) Genes corresponding to each KEGG pathway shown were correlated both within and across all other organs where the number of genes which meet each Student’s p-value threshold are shown (y-axis). Tissues (x-axis) are rank-ordered by the number of genes which correlate for hsa04062 − chemokine signaling pathway at p-value <0.01 and shown for other KEGG terms, hsa04640 − hematopoietic cell lineage (C) and hsa00190 − oxidative phosphorylation (D) and additionally p-value <1e-4 (right side).

Figure 3.. Inter-individual transcript correlation structure and network architecture of liver PCSK9 and adipose PNPLA2.

(A) Distribution of pan-tissue genes correlated with liver PCSK9 expression (q<0.1), where 93% of genes were within liver (purple). (B) Gene ontology (GO) (BP) overrepresentation test for the top 500 hepatic genes correlated with PCSK9 expression in liver. (C) Undirected network constructed from liver genes (aqua) correlated with PCSK9, where those annotated for ‘cholesterol biosynthetic process’ are colored in red. (D–E) Overrepresentation tests corresponding to the top-correlated genes with adipose (subcutaneous) PNPLA2 expression residing in adipose (D) or peripherally in skeletal muscle (E). (F) Undirected network constructed from the strongest correlated subcutaneous adipose tissue (light aqua) and muscle genes (light brown) with PNPLA2 (black), where genes corresponding to GO terms annotated as ‘fatty acid beta oxidation’ or ‘muscle contraction’ are colored purple or red, respectively. For these analyses all 310 individuals (across both sexes) were used and q-value adjustments calculated using a Benjamini-Hochberg FDR adjustment. Network graphs generated based in biweight midcorrelation coefficients, where edges are colored blue for positive correlations or red for negative correlations. Network edges represent positive (blue) and negative (red) correlations and the thicknesses are determined by coefficients. They are set for a range of bicor=0.6 (minimum to include) to bicor=0.99.

Figure 4.. Hybrid mouse diversity panel (HMDP) tissue- and diet-specific correlations of sex hormone receptors.

The top 10 phenotypic traits are shown for correlations with expression of androgen receptor (A), estrogen receptor 1 (B), or estrogen receptor 2 (C), colored by direction in the hybrid mouse diversity panel. Positive correlations are shown in light blue and negative correlations as sunset orange, where phenotypes (y-axis) are ordered by significance (x-axis, -log10(p-value) of correlation). Correlations are segregated by whether sex hormone receptors are expressed by gonadal adipose tissue (left two columns) in ~100 HMDP strains fed a high-fat/high-sucrose (HFHS) diet (left), normal chow diet (left middle), or liver-expressed receptors fed an HFHS diet (right middle) or normal chow diet (right). Dashed lines show a Student’s correlation p-value (from bicor) of p=0.05 (purple) or p=0.001 (beige).