Parent-of-origin effects propagate through networks to shape metabolic traits

Abstract
Editor's evaluation
Introduction
Results
Discussion
Materials and methods
Data availability
References
Article and author information
Metrics

Abstract

Parent-of-origin effects are unexpectedly common in complex traits, including metabolic and neurological traits. Parent-of-origin effects can be modified by the environment, but the architecture of these gene-by-environmental effects on phenotypes remains to be unraveled. Previously, quantitative trait loci (QTL) showing context-specific parent-of-origin effects on metabolic traits were mapped in the F₁₆ generation of an advanced intercross between LG/J and SM/J inbred mice. However, these QTL were not enriched for known imprinted genes, suggesting another mechanism is needed to explain these parent-of-origin effects phenomena. We propose that non-imprinted genes can generate complex parent-of-origin effects on metabolic traits through interactions with imprinted genes. Here, we employ data from mouse populations at different levels of intercrossing (F₀, F₁, F₂, F₁₆) of the LG/J and SM/J inbred mouse lines to test this hypothesis. Using multiple populations and incorporating genetic, genomic, and physiological data, we leverage orthogonal evidence to identify networks of genes through which parent-of-origin effects propagate. We identify a network comprised of three imprinted and six non-imprinted genes that show parent-of-origin effects. This epistatic network forms a nutritional responsive pathway and the genes comprising it jointly serve cellular functions associated with growth. We focus on two genes, Nnat and F2r, whose interaction associates with serum glucose levels across generations in high-fat-fed females. Single-cell RNAseq reveals that Nnat expression increases and F2r expression decreases in pre-adipocytes along an adipogenic trajectory, a result that is consistent with our observations in bulk white adipose tissue.

Editor's evaluation

We all learned simple Mendelian Punnett Squares in Junior High or earlier when studying simple Mendelian traits. But we also all know that the world is so much richer and more complex than this. The current article explores some of that complexity, opening rich insights into intergenerational effects, offering the opportunity for mathematical thinking and further hypothesis testing, and opening up exciting new hypotheses to test. As Professor Stephen Stearns wrote, "Many of us do not do science only, or even primarily, to achieve practical results. We do it because we are fascinated with neat ideas. Evolutionary medicine is full of them, including parent-of-origin pattern." Let us enjoy the wonder.

https://doi.org/10.7554/eLife.72989.sa0

Introduction

Parent-of-origin effects, where the phenotypic effect of an allele depends on whether the allele is inherited maternally or paternally, are epigenetic phenomena associated with a wide range of complex traits and diseases (Lawson et al., 2013). Thus, the functional impact of a specific genetic variant can depend on its parental origin. The best characterized parent-of-origin effect is genomic imprinting, an epigenetic process in which either the maternally or paternally inherited allele is silenced, typically through DNA methylation. In humans, there are 107 verified imprinted genes and in mice there are 124, of which ~ 70% overlap (Jirtle, 2012). Despite the rarity of imprinted genes, parent-of-origin effects on complex traits and diseases are relatively common, suggesting that canonical imprinting mechanisms are not sufficient to account for these phenomena (Mozaffari et al., 2019; Zeng et al., 2019). With so few imprinted genes, what mechanisms underlie these parent-of-origin effects? We hypothesize that a small number of imprinted genes can generate a large number of parent-of-origin effects through interactions with non-imprinted genes.

In this study, we use four populations at different levels of intercrossing of the LG/J and SM/J inbred mouse lines to test the hypothesis that non-imprinted genes can contribute to parent-of-origin effects on metabolic phenotypes through epistatic interactions with imprinted genes. Multiple populations (F₀, F₁, F₂, F₁₆) allow us to refine our search space and provide orthogonal evidence supporting putative networks of interacting genes. Metabolic traits were previously mapped in a F₁₆ generation of an advanced intercross between LG/J and SM/J (Cheverud et al., 2011; Lawson et al., 2010; Lawson et al., 2011a; Lawson et al., 2011b). We generated visceral white adipose tissue gene expression profiles from 20 week-old F₁ animals in order to match the age of the F₁₆ LG/J x SM/J advanced intercross population. F₁ reciprocal cross (LxS and SxL) mice were subjected to the same high and low-fat diets and phenotyping protocols as the previously-studied F₁₆ mice to keep environmental contexts consistent. We identified genes showing parent-of-origin-dependent allele-specific expression (ASE), characterized interactions among these genes and biallelic genes that are differentially expressed by reciprocal cross (DE), and correlated interacting ASE and DE gene pairs with metabolic phenotypes in the F₁ population. Pairs that significantly associated with phenotypic variation were tested for epistasis on correlated traits in the F₁₆ population.

We identified an epistatic network that forms a nutritional environment responsive pathway mediated through calcium signaling. This network contributes to metabolic variation by balancing proliferation, differentiation, and apoptosis in adipocytes. The genes comprising this network jointly serve functions associated with growth in multiple tissues, which is consistent with the evolutionary hypothesis that sexual conflict underlies some parent-of-origin effects (Mochizuki et al., 1996). We focus on two key interacting genes: Nnat (neuronatin), a canonically imprinted gene, and F2r (coagulation factor II receptor), a biallelic gene showing significant DE by cross in F₁ high-fat-fed female animals. Co-expression of these two genes associates with variation in basal glucose levels, and this association persists across generations. Further, single-cell RNAseq reveals that Nnat expression increases and F2r expression decreases in pre-adipocytes along an adipogenic trajectory, a pattern consistent both with their expression in bulk white adipose tissue and with their respective roles in adipogenesis. Our results demonstrate that incorporating orthogonal lines of evidence including genotype, allele specific expression, total gene expression, single-cell expression, and phenotype from different populations varying in their degree of intercrossing is a powerful way to identify putative mechanisms and test hypotheses underlying parent-of-origin effects on phenotype.

Results

Non-imprinted genes interact with imprinted genes and effect metabolic phenotypes

We test the hypothesis that non-imprinted genes can mediate complex parent-of-origin effects on phenotypes through genetic interactions with imprinted genes using a F₁ reciprocal cross model of the LG/J and SM/J inbred mice (LxS and SxL). In this model the effects of parental origin on an allele can be tested directly and isolated from sequence dependent cis-regulatory differences. We validated our findings in LG/J and SM/J parentals (F₀) as well as in F₂ and F₁₆ intercrosses of LGxSM (Figure 1). The parental F₀ animals serve to anchor variation in allele-specific expression that is a function of allelic identity (L or S). Incorporating the F₂ and F₁₆ populations into our validations ensures that the interactions we observe are not solely a function of linkage in the F₁ animals. We generated mRNA expression profiles in white adipose tissue from 20-week-old F₁ reciprocal cross animals. These animals were subjected to the same high and low-fat diets and phenotyping protocols as the previously studied F₁₆ animals (Cheverud et al., 2011; Lawson et al., 2011a, Lawson et al., 2010; Carson et al., 2020; Miranda et al., 2020). We identified two classes of genes: (1) imprinted genes and (2) non-imprinted genes with parent-of-origin effects on total expression.

Figure 1

Download asset Open asset

Proposed model for propagation of parent-of-origin effects through gene-gene interactions.

Parent-of-origin effects should be partitioned into *cis* mechanisms and *trans* mechanisms (A) An example of a *cis* parent-of-origin effect is a system with three regulatory elements: promoter, insulator, and enhancer. Activation of transcription requires the enhancer to act upon the promoter. Enhancer activity is blocked by the insulator when it has been bound by CTCF. CTCF cannot be bound when methylated. In this system, the insulator is selectively methylated when inherited maternally, so methylation of the maternally inherited insulator blocks CTCF binding, allowing the enhancer to activate transcription. Because the paternally inherited insulator is not methylated, it is bound by CTCF which blocks enhancer activity, silencing transcription. This canonical genomic imprinting mechanism interacts with genetic variation in the three regulatory features. For example, if one allele produces stronger enhancer activity (Alt) than the other, individuals inheriting the Alt allele maternally would have elevated expression compared to those that inherit the same allele paternally. These *cis* genetic effects do not occur in isolation. Due to the highly interconnected nature of biological systems, there are downstream effects. We refer to these as *trans* parent-of-origin effects. (B) An example of a *trans* parent-of-origin effect is a system with two genes each having its own promoter. The first gene is canonically imprinted, and the activity of the gene promoter is blocked by DNA methylation. The imprinted gene’s promoter is methylated when inherited maternally. Consequently, the paternally inherited allele is almost exclusively expressed. As before, when genetic variation in a regulatory feature interacts with these epigenetic mechanisms, we see parent-of-origin effects on expression of the imprinted gene. In this example, the imprinted gene regulates expression of a non-imprinted gene. Despite the non-imprinted gene being agnostic to parental origin, its expression nonetheless depends on the parental origin of alleles at the imprinted locus. (C) Summary of our experimental design. Expression patterns of genes showing allele-specific expression (ASE) such as imprinted genes are shaped by parental genotypes and environment (e.g. nutrition). Downstream gene expression is a function of their genotype and the expression of upstream ASE genes. Altered parent-of-origin dependent total gene expression of ASE genes leads to differential expression of downstream genes varying only in allelic parent-of-origin (DE). Phenotype is most directly affected by expression of DE genes. Variation in DE gene expression leads to corresponding variation in phenotype. Mouse populations used to probe parts of this model are labeled F₀ (inbred lines), F₁ (reciprocal cross of inbred lines), F₂ (intercross of F₁ mice), and F₁₆ (advanced intercross of inbred lines).

To test our model, we identified genes showing parent-of-origin dependent allele specific expression (ASE). We identified 23 genes showing significant ASE (Figure 2A; Supplementary file 1). Of these 23 genes, 17 are canonically imprinted genes, two are not reported as imprinted genes but are located in known imprinted domains, and four are novel. Next, we identified genes showing differential total expression between individuals varying only in allelic parent-of-origin (DE between reciprocal crosses, SxL vs LxS). We identified 33 genes that are significantly DE in at least one sex or dietary context (Figure 2A; Supplementary file 2). A larger set of genes show signatures of parent-of-origin effects at the total gene expression level, but do not meet the statistical rigor demanded by the massive multiple tests burden incurred by a genome-wide scan accounting for sex, diet, and parent-of-origin (see Materials and methods).

Figure 2 with 13 supplements see all

Download asset Open asset

Genes showing parent-of-origin effects at the allele specific and/or total expression levels covary with each other and with metabolic traits.

(A) Mean parent-of-origin effect score across contexts. Effect size of ASE is calculated as the mean allelic bias (L / L + S) of SxL animals minus LxS animals. Effect size of DE is measured by log₂(Fold Change) between LxS and SxL crosses. The single context with largest magnitude fold change is plotted for each gene. Dashed lines represent minimum acceptable effect size cut-offs within a context. Genes showing significant ASE and sufficiently large parent-of-origin effect score are shown in blue. Genes showing significant DE and sufficiently large fold change in some sex or dietary context are shown in lime. Genes showing both ASE and DE are shown in teal. Genes not meeting cut-offs are shown in gray. The two genes showing significant ASE but falling short of parent-of-origin effect score requirements are a case of context dependent bipolar parent-of-origin effect scores (i.e. paternally expressed in one context and maternally expressed in its opposite). (B) Parent-of-origin effect network constructed from ASE and DE gene pairs. (C) Significantly overrepresented ontologies after multiple tests correction in the parent-of-origin effect network. Terms are color coded by ontology domain. GO biological process (yellow), GO cellular component (orange), and Mammalian phenotype (purple). Circle size denotes the number of genes with each term. (D) Correlation of parent-of-origin effect network genes with metabolic traits. Only genes and phenotypes with at least one significant correlation after multiple test corrections are shown. The heatmap is broken up into subnetworks with the ASE gene as the first separated row followed by associated DE genes in subsequent rows. Columns correspond to metabolic traits. Coloration of each cell denotes the Pearson’s correlation coefficient value.

To identify interactions between gene sets, we constructed a network comprised of genes that could initiate a parent-of-origin effect on phenotype (ASE) and genes that may mediate the effect onto phenotype (DE). Interacting gene pairs were predicted by modeling the expression of biallelic genes that are significantly DE by reciprocal cross as a function of the expression of genes showing significant parent-of-origin-dependent ASE, their allelic bias (L_bias), diet, sex, and the diet-by-sex interaction. Genes showing parent-of-origin effects form a highly interconnected network comprised of 52 genes forming 217 gene pairs (Figure 2B)(Supplementary file 3). Most of these interactions are trans-chromosomal. We identified two genes that could serve as initiation points of propagating parent-of-origin effects through this network. These two genes, Nnat (neuronatin) and Cdkn1c (cyclin dependent kinase inhibitor 1 C), are both canonically imprinted and differentially expressed by reciprocal cross (Supplementary file 1).

Functional over representation analysis was performed and seven terms were significantly overrepresented at an FDR ≤ 0.05 (Figure 2C; Zhang et al., 2005). Enriched terms suggest this network plays a role in signaling and genetic imprinting (Supplementary file 4). In order to identify which phenotypes might be affected by genes in this network, gene expression was correlated with metabolic phenotypes collected for the F₁ animals (Figure 2D). Seventy-four ASE/DE/phenotype sets were identified as candidates for subsequent testing (Supplementary file 5).

Epistasis in an F₁₆ advanced intercross identifies a diet-responsive network affecting adipogenesis

To validate the interactions we identified in F₁ animals, we tested for imprinting-by-imprinting epistasis in an F₁₆ population. Imprinting-by-imprinting epistasis occurs when the parent-of-origin effect at a locus is dependent on the parent-of-origin of alleles at another locus. This allowed us to determine if the effect of parent-of-origin at DE genotype on phenotype is dependent upon the parent-of-origin at ASE genotype. This orthogonal approach allows us to connect genotype at these loci to phenotype as predicted in the F₁ candidates. Nine epistatic interactions replicated in the F₁₆ population (n = 1002 animals, FDR ≤ 0.1; Figure 3A; Supplementary file 6). These interactions were comprised of three ASE genes showing parent-of-origin (Cdknlc, Nnat, Plcd1), six genes that are DE by cross (Car3, F2r, Hexb, Hmger, Srgn, Tnfrsf11a) and four phenotypes (basal glucose level, AUC calculated from a glucose tolerance test, serum cholesterol, necropsy weight). Together, these nine genes form a putative diet-responsive network affecting adipogenesis (Figure 3B).

Figure 3 with 3 supplements see all

Download asset Open asset

Interacting genes form a diet-responsive network affecting adipogenesis.

(A) There are nine significant imprinting-by-imprinting epistatic ASE/DE/phenotype sets in the F₁₆ advanced intercross population (n = 1002). Interactions are shown as lines connecting ASE (yellow) and DE genes (purple). Chromosome number is shown around the plot. (B) The epistatic parent-of-origin effect network is comprised of key steps in a putative pathway regulating differentiation and survival of adipocytes. This pathway was constructed by incorporating previously published cellular functions. The pathway members are color coded in blue for ASE genes (*Plcd1*, *Nnat*, and *Cdkn1c*) and green for DE genes (*F2r*, *Hexb*, *Hmgcr*, *Car3*, *Tnfrsf11a*, and *Srgn*). The network breaks down into potentiation, transduction, and response. *Nnat* and *Hexb* potentiate signaling by managing availability and accumulation of calcium necessary for signal transduction. Once a signal is received, *F2r* and *Plcd1* transduce it by activating second messengers to initiate a response. This response initiates an adipogenesis cellular program that affects expression of *Cdkn1c*, *Hmgcr*, *Car3*, *Tnfrsf11a*, and *Srgn*.

The network can be broken down into signal potentiation, transduction, and response. Nnat (neuronatin) and Hexb (beta-hexosaminidase subunit beta) fall into the potentiation group. These genes play a role in managing the availability and accumulation of calcium necessary for signal transduction. Nnat is a paternally expressed canonically imprinted gene which encodes a proteolipid protein that localizes to the ER (Li et al., 2010). Nnat is diet-responsive and its overexpression in 3T3L1 pre-adipocytes promotes adipogenesis through increased free cytosolic calcium (Young et al., 2005). In pre-neural stem cells, Nnat binds sarco/endoplasmic reticulum Ca²⁺-ATPase (SERCA) to block Ca²⁺ uptake into the ER thereby increasing cytosolic Ca²⁺ levels (Lin et al., 2010). In addition to Nnat, Hexb regulates the uptake and accumulation of Ca²⁺ in the ER via SERCA (Pelled et al., 2003). Upon the arrival of a signal, F2r (coagulation factor II receptor) and Plcd1 (1-phosphatidylinositol 4,5-bisphosphate phosphodiesterase delta-1) in the transduction group initiate the adipogenesis cellular program. F2r is a G-protein-bound receptor that promotes phosphoinositide hydrolysis (Soh et al., 2010). Variation in the human F2R gene is associated with obesity (Kichaev et al., 2019). G-protein-coupled receptors transmit external signals into the cell where they are then propagated by secondary messenger systems, one of which is mediated by Plcd1 (Nakamura et al., 2005; McDonald and Mamrack, 1995). The downstream effect of Plcd1-mediated signaling is the efflux of calcium into the cytosol from the ER, thereby increasing cytosolic Ca²⁺ levels (Thatcher, 2010; Berridge, 2016). Increased cytosolic Ca²⁺ in pre-adipocytes promotes phosphorylation of cAMP-response element-binding protein (CREB), which promotes activity of CCAAT/enhancer-binding protein (C/EBP) transcription factors, activating adipogenesis, altering the expression of Cdkn1c (cyclin dependent kinase inhibitor 1 C), Hmgcr (3-hydroxy-3-methylglutaryl-CoA reductase), Car3 (carbonic anhydrase 3), Tnfrsf11a (TNF receptor superfamily member 11 a), and Srgn (serglycin).

Cdkn1c is a canonically imprinted maternally expressed gene that inhibits cell proliferation (Kang et al., 2008). Increased expression of Cdkn1c is protective against diet-induced obesity in mice (Van de Pette et al., 2018), and in humans increased caloric intake results in decreased CDKN1C expression (Franck et al., 2011). Hmgcr is the rate-limiting enzyme in cholesterol biosynthesis (Burg and Espenshade, 2011; Jo and Debose-Boyd, 2010) and converts HMG-CoA into mevalonate, which is essential for adipocyte survival (Yeh et al., 2018). Srgn is an adipocytokine thought to be part of a feedback loop with Tnfα (tumor necrosis factor alpha), mediating paracrine cross-talk between macrophages and adipocytes (Lemire et al., 2007; Imoto-Tsubakimoto et al., 2013; Schick et al., 2001; Zernichow et al., 2006). Srgn is known to play a role in osteoblast-mediated bone mineralization (Bigdeli et al., 2010), which along with osteoclast-driven bone deconstruction drives bone remodeling (Aubin, 1992). Osteoblasts share a lineage with adipocytes, and the quantity of osteoblasts is inversely proportional to that of marrow adipose tissue (Rodríguez et al., 2008; Prockop, 1997; Ali et al., 2005; Akune et al., 2004; Cho et al., 2011; Rosen and Bouxsein, 2006; Turner et al., 2018). Tnfrsf11a is a cell surface protein that regulates differentiation of osteoclasts (Nakagawa et al., 1998). Osteoprotegerin (OPG) is a decoy receptor for TNFRSF11A thereby inhibiting osteoclastogenesis and bone resorption (Matsuo et al., 2020). OPG is expressed during differentiation of 3T3L1 adipocytes (An et al., 2007). Expression of OPG is induced by Tnfα in 3T3L1 adipocytes and is associated with obesity in humans (Holecki et al., 2007; Erol et al., 2016; Zaky et al., 2022).

The exact function of OPG/Tnfrsf11a outside of osteoclastogenesis is unknown, but the function of osteoclasts is to break down bone tissue during bone resorption. Bone resorption regulates the level of blood calcium. The bioavailability of calcium in the blood potentially alters ER calcium stores, creating cross-talk between bone cells and white adipose tissue calcium signaling. Osteoclasts break down bone by acidifying mineralized bone, orchestrated by osteoblasts that have become embedded in the matrix they produce (osteocytes). Oxidative stress on osteocytes from the bone acidification process is prevented by Car3. Car3 is an enzyme that catalyzes the conversion of carbonic acid to CO₂ and water. Its expression in white adipose is negatively correlated with, and responsive to, long-term obesity in mice and humans (Stanton et al., 1991; Font-Clos et al., 2017). Car3 does not protect against diet-induced obesity and is not necessary for fatty acid synthesis (Renner et al., 2017). As such its exact function in adipocytes is unknown.

Nnat and F2r covary in white adipose tissue and their interaction associates with variation in basal glucose levels across generations

To better understand how these interactions affect phenotype, we focused on the negative correlation of the imprinted gene, Nnat, and the biallelic gene, F2r, in the above network in high fat-fed females, the cohort with the most significant covariation in the F₁ animals (FDR = 0.05). Nnat and F2r show significant imprinting-by-imprinting epistasis for basal glucose levels in the F₁₆ population (FDR = 6.00e^–16; Figure 4A and B). To validate gene expression patterns, we combined F₁ biological replicates and F₀ high fat-fed female animals (F₁ n = 13 and F₀ n = 12) and again observe that F2r and Nnat are each significantly differentially expressed between reciprocal heterozygotes, that is by cross (Figure 4C and D). Further, the co-expression of Nnat and F2r also persists in the F₀/F₁ population (Figure 4E).

Figure 4 with 1 supplement see all

Download asset Open asset

*Nnat* and *F2r* covary across generations.

(A) Breeding scheme for the F₁₆ Advanced Intercross between the LG/J and SM/J inbred strains. (B) Significant imprinting-by-imprinting epistasis associated with variation in basal glucose (n = 1002). The parent-of-origin effects of *F2r* on basal glucose depend on the parent-of-origin effects at *Nnat*. (C) Expression of *Nnat* across genotypes in a combined F₀/F₁ population of high fat-fed females (n = 25). (D) Expression of *F2r* across genotypes in a combined F₀/F₁ population of high-fat-fed females (n = 25). (E) Significant correlation between *Nnat* and *F2r* expression in the F₀/F₁ mice (F₁ n = 13; F₀ n = 12). (F) and (G) Correlations between basal glucose and *Nnat* and *F2r* in the F₀/F₁ mice (F₁ n = 13; F₀ n = 12). (H) Significant correlation between *Nnat* and *F2r* expression in the high fat-fed female F₂ mice (n = 14). (I) and (J) Correlations between basal glucose and *Nnat* and *F2r* are not individually significant in the F₂ mice. Alleles are ordered maternal | paternal within the genotype classes.

A limitation of identifying covariation patterns in F₁ and F₀ populations is that all loci are linked. This makes it difficult to determine which ASE genes truly co-express with DE genes. While incorporation of orthogonal F₁₆ genotypes and phenotypes helps reduce false discoveries, a population with randomized genetic background for which we have expression data is needed to replicate these results. To that end, F₂ animals were generated and Nnat and F2r gene expression levels were measured via qPCR (n = 14). We found that F2r and Nnat are significantly co-expressed in high-fat-fed female F₂ animals (Figure 4H).

F2r expression significantly positively correlates with basal glucose levels in the RNA-sequenced high-fat-fed female F₁ animals (r = 0.514, FDR = 0.01; Supplementary file 5). F2r expression is also significantly positively correlated with basal glucose in the combined F₀/F₁ population (Figure 4G). A negative trend between Nnat expression and basal glucose level is observed but not statistically significant in the combined F₀/F₁ animals (Figure 4F). Correlation of F2r’s and Nnat’s individual expression with basal glucose in F₂ mice follows the same pattern as in the F₀/F₁’s. Bootstrapping to calculate confidence intervals shows that the correlation differences between F₀/F₁ and F₂ are not significant (Figure 4I and J; Figure 4—figure supplement 1). However, the product of Nnat and F2r expression (Nnat x F2r) is significantly predictive of basal glucose (p = 0.045, R² = 0.29). This indicates that expression of Nnat and F2r, as a function of their genotypes and allelic parent-of-origin, are not individually sufficient to explain variation in basal glucose levels. But together they are able to explain a significant amount of phenotypic variation. This is precisely what our epistatic model would predict.

Finally, studying the F₂ animals allows us to determine if maternal mitochondrial ancestry contributes significantly to Nnat or F2r expression or to variation in basal glucose. We find mitochondrial genome identity does not significantly covary with F2r expression (p = 0.198), Nnat expression (p = 0.365), or basal glucose (p = 0.388).

Single-cell RNAseq reveals that Nnat expression increases and F2r expression decreases in pre-adipocytes along an adipogenic trajectory

To determine what cell types express Nnat and F2r and whether the directionality of the Nnat imprinted → F2r target correlation persists along the adipogenic trajectory, we turned to single-cell RNAseq. We used publicly available scRNAseq data collected from stromal vascular cells isolated from C57BL/6 J epididymal adipose tissue (Burl et al., 2018). Cell type identity was assigned using previously reported markers for this data set (Adipoq = differentiating mesenchymal stem cells; Pdgfra = mesenchymal stem cells; Csf1r = macrophage; Cdh5 = vascular endothelial cells; Acta2 = vascular smooth muscle cells; Cd2 = B cells) (Supplementary file 7; Figure 5—figure supplement 1). The adipogenic trajectory refers to cells transitioning from pre-adipocytes (mesenchymal stem cells) to cells differentiating into adipocytes. Clusters along this trajectory were identified by the opposing expression patterns of Pdgfra and Adipoq (Figure 5A-D and I). We found that Nnat expression increases along the trajectory while F2r expression decreases (Figure 5E–F and H). Further there is a negative association between Nnat and F2r expression within adipocytes along the trajectory (Figure 5G). This pattern is consistent with the negative correlation we observe between Nnat and F2r in the bulk white adipose tissue. Because available scRNAseq data do not match the exact sex/ diet/ genetic background contexts of the LGxSM mice, there will be unaccounted for differences between the data sets. However, the observed consistent pattern indicates that the pathway structure persists across sex/ diet/ genetic backgrounds.

Figure 5 with 5 supplements see all

Download asset Open asset

*Nnat* expression increases and *F2r expression decreases* in pre-adipocytes along an adipogenic trajectory.

(A) *Adipoq* is a marker of adipocytes whose expression (purple) increases along the trajectory. (B) *Pdgfra* is a marker of mesenchymal stem cells whose expression (pink) decreases along the trajectory. (C) Cells in clusters expressing one or both *Adipoq* and *Pdgfra* fall along an adipogenic trajectory. (D) Intensity of expression of *Adipoq* and *Pdgfra* indicated by coloration. (E) *Nnat* expression (blue) increases along the trajectory. (F) *F2r* expression (teal) decreases along the trajectory. (G) Negative association between *Nnat* and *F2r* expression within adipocytes along the trajectory. (H) Intensity of expression of *Nnat* and *F2r* indicated by coloration. (I) The adipogenic trajectory is broken into subclusters of cells with no *Adipoq* expression (cluster 0) to high *Adipoq* expression (cluster 2).

In addition to interrogating Nnat and F2r in single cells along an adipogenic trajectory, we found that eight of the nine genes comprising the epistatic parent-of-origin effect network described above are differentially expressed along the trajectory, and they associate with cell types that are consistent with their respective roles in adipose tissue (Figure 5—figure supplement 2).

Discussion

Epistatic interactions between imprinted and non-imprinted genes can influence complex traits when the genotypic effects of one gene depends on the parent-of-origin of alleles at another (Lawson et al., 2013; Wolf and Cheverud, 2009). Here, we examined epistatic interactions associated with parent-of-origin effects on dietary-obesity traits in white adipose using a simple yet powerful F₁ reciprocal cross mouse model. Although these parent-of-origin dependent allele-specific expression biases are consistent with imprinting mechanisms, we cannot rule out that maternal and/or paternal effects also contribute to the phenomena we observe (Hager et al., 2008).

Interactions between imprinted and non-imprinted genes have previously been shown to contribute to variation in metabolic phenotypes. For example, the maternally expressed transcription factor KLF14 (kruppel-like factor 14) regulates biallelic gene expression related to adiposity (Parker-Katiraee et al., 2007; Small et al., 2011). Mapping studies have identified two SNPs (rs4731702, rs972283) upstream of KLF14 associated with type II diabetes and cholesterol levels (Voight et al., 2010; Teslovich et al., 2010). Both variants have maternally restricted cis-regulatory associations with KLF14 expression in adipose tissue (Kong et al., 2009). eQTL analysis found that rs4731702 is also enriched for trans-associations with KLF family transcription factor binding sites in subcutaneous white adipose tissue, suggesting that KLF14 may be a master transcriptional regulator in adipose tissue (Small et al., 2011). Whether additional pairs of imprinted and biallelic genes are similarly co-expressed and associate with phenotypic variation remains an open question that has not been thoroughly investigated in large landmark functional genomics studies including ENCODE, GTEx, and GWAS, leaving a significant gap in our knowledge. Interactions between imprinted and biallelic genes could explain some of the observed parent-of-origin effect patterns associated with regions lacking obvious candidate genes, as described in a recent survey of 97 complex traits measured in outbred mice (Mott et al., 2014).

Our model asserts that parent-of-origin effects start at ASE genes and are transduced through DE genes onto phenotype. This is illustrated in the interaction between Nnat and F2r. If a cis-regulatory effect interacts with epigenetic modifications (i.e. imprinting) at Nnat, then Nnat expression of genotypic classes are affected by paternal allele identity (Lawson et al., 2013). Between the LG/J and SM/J alleles at Nnat, the LG/J allele is more highly expressed. If our model is correct, the downstream DE gene should show a corresponding pattern (Figure 1B). In the case of Nnat and F2r, which have strong negative correlated expression, when the LG/J allele is inherited paternally at Nnat, the higher expression of Nnat should correspond with lower expression of F2r. This is what we observe (Figure 4). If this relationship is true, we should see persistent co-expression of Nnat and F2r across genetic backgrounds (F₀, F₁, F₂), which we do (Figure 4). This supports a biologically meaningful relationship between Nnat and F2r. Our model further predicts that the DE genes should more closely affect phenotype (Pierce et al., 2014; Shan et al., 2019; Lutz and Hokanson, 2015). In the case of Nnat and F2r, we expect F2r to more strongly associate with basal glucose levels than Nnat, which we observe (Figure 4).

There is a clear relationship between Nnat and F2r in adipogenesis, but the specifics of how this relationship extends to glucose homeostasis are unclear. One possibility is that by altering SERCA function, Nnat affects not only the formation of new adipocytes, but also the beiging of adipocytes. The SERCA channel is uncoupled in beige adipocytes as part of a UCP1-independent form of non-shivering thermogenesis (Ikeda and Yamada, 2020). Non-shivering thermogenesis consumes a significant amount of energy, thereby altering glucose homeostasis (Carson et al., 2020). This hypothesis links these genes to physiological processes that are consistent with evolutionary hypotheses about the prevalence of parent-of-origin effects. Nnat and F2r are members of a putative network we identified that is coordinated by interactions between ASE and DE genes. From the literature, we found that the genes in this epistatic network function in key steps in a pathway regulating differentiation and survival of adipocytes in response to nutritional environment (Figure 3B). Specifically, there is evidence that it plays a critical role in the induction of adipogenesis. This alone demonstrates how parent-of-origin effects can move through networks along molecular pathways. Beyond proof-of-principle this network provides a clue to the puzzle of the prevalence of parent-of-origin effects.

The constituents of this single network appear to play vastly different physiological roles depending on the tissue. In white adipose the network appears to play some role in balancing proliferation, differentiation, and apoptosis as we describe above. In pancreatic ß-cells, members of this network affect insulin secretion (Millership et al., 2018). In bone, members of this network affect the balance of cartilage/bone growth and reabsorption. These three physiological processes may at first seem unrelated, but they share one key commonality – they are jointly critical to growth. This is consistent with the sexual conflict hypothesis attributed to parent-of-origin effects (Patten et al., 2014; Babak et al., 2015). The of size of progeny in placental mammals can have opposing fitness consequences for mothers/ maternal relatives and fathers/ paternal relatives. The fitness of fathers and paternal relatives, particularly in the case of multi-paternity litters, is improved with larger progeny (Mochizuki et al., 1996; Babak et al., 2015; Fowden and Moore, 2012; Patten et al., 2014; Wilkins and Haig, 2003; Haig, 1997). This comes at a fitness disadvantage to the mother. The fitness of mothers is improved by progeny of a manageable size, allowing her to produce multiple litters.

According to this model, imprinting evolved in part to allow one parental lineage to hijack parts of a nutritional environment response pathway driving growth in a direction favorable to maximize the fitness of that lineage. Key processes in such a pathway driving growth would include the secretion of growth factors, construction of cartilage and bone, and the accumulation of energy stores. We present a network that appears to play a role in all three processes. If the sexual conflict hypothesis is true, then the most parsimonious place for imprinting to evolve would be in key regulatory points that affect as many aspects of growth as possible. This is consistent with the network we identified, a single pathway affecting many aspects of growth. This hints at the possibility that parent-of-origin effects are common because of the multi-purpose nature of the pathways in which genomic imprinting manifests and parent-of-origin effects propagate.

By leveraging the reciprocal F₁ hybrids, we are able to integrate parent-of-origin-dependent allele specific expression and parent-of-origin-dependent differential expression with F₁₆ phenotypes. By doing so, we identify plausible candidates for functional validation and describe discrete molecular networks that may contribute to the observed parent-of-origin effects on metabolic phenotypes. The genes and interactions we present here represent a set of actionable interacting candidates that can be probed to further identify the machinery driving these phenomena and make predictions informed by genomic sequence. The frameworks we have developed account for the genetic, epigenetic, and environmental components underlying these parent-of-origin effects, thereby improving our ability to predict complex phenotypes from genomic sequence. We focused on metabolic phenotypes in this study, but the patterns we identified may translate to other complex traits where parent-of-origin effects have been implicated.

Materials and methods

Key resources table

Reagent type (species) or resource	Designation	Source or reference	Identifiers	Additional information
Other	High fat diet	Teklad	TD88137	42% kcal from fat
Other	Low-fat diet	Research Diets	D12284	15% kcal from fat
Commercial assay, kit	RNeasy Lipid Tissue Kit	QIAgen	74,804
Commercial assay, kit	RiboZero kit	Illumina	20040529
Commercial assay, kit	DNA 1000LabChip	Agilent	5067–1504
Commercial assay, kit	High-Capacity cDNA Reverse Transcription Kit	Thermo Fisher	4368814
Sequence-based reagent	Nnat forward primer	This paper		Detailed information is found in the methods section
Sequence-based reagent	Nnat reverse primer	This paper		Detailed information is found in the methods section
Sequence-based reagent	F2r forward primer	This paper		Detailed information is found in the methods section
Sequence-based reagent	F2r reverse primer	This paper		Detailed information is found in the methods section
Sequence-based reagent	L32 forward primer	This paper		Detailed information is found in the methods section
Sequence-based reagent	L32 reverse primer	This paper		Detailed information is found in the methods section
Software, algorithm	R	R	3.6.1
Software, algorithm	STAR	STAR	DOI: 10.1093/bioinformatics/bts635
Software, algorithm	FASTQC	FASTQC	other	https://www.bioinformatics.babraham.ac.uk/projects/fastqc/
Software, algorithm	EdgeR	CRAN	DOI: 10.1093/bioinformatics/btp616
Software, algorithm	WEB-based Gene Set Analysis Toolkit	WEB-based Gene Set Analysis Toolkit	DOI: 10.1093/nar/gkz401
Software, algorithm	Seurat	Seurat	DOI: 10.1038/nbt.3192
Strain, strain background (Mus musculius)	SM/J	The Jackson Laboratory	000687
Strain, strain background (Mus musculus)	LG/J	The Jackson Laboratory	000675

Mouse husbandry and phenotyping

Request a detailed protocol

LG/J and SM/J founders (F₀) were obtained from The Jackson Laboratory (Bar Harbor, ME). F₁ reciprocal cross animals were generated by mating LG/J mothers with SM/J fathers (LxS) and the inverse (SxL). F₂ intercrossed animals were generated by mating LxS mothers with SxL fathers and the inverse. After weaning at 21 days, animals were separated into sex-specific cages of 3–5 animals per cage and randomly placed on high-fat (42% kcal from fat; Teklad TD88137) or low-fat (15% kcal from fat; Research Diets D12284) isocaloric diets. Feeding was ad libitum. There were 96 experimental F₁ animals in total, with 48 animals for each cross (LxS and SxL). Within each cross, there were 24 high-fat-fed animals (12 males; 12 females) and 24 low-fat-fed animals (12 males; 12 females). The F₂ animals were generated for a different study, following the same weaning protocol and diets, and we used data from the high fat-fed females (n = 14) for validation in the the current study (Carson et al., 2020). Additionally, we used data generated from founder F₀ (LG/J (n = 6) and SM/J (6)) high fat-fed female animals, also generated for a different study and subjected to the same weaning protocol and diets (Carson et al., 2020). The barrier animal facilities at WUSM follow a 12/12 hr light/dark schedule, all water is autoclaved and changed weekly, and all cages are changed weekly.

All animals were weighed weekly from three weeks of age until sacrifice. At 19 weeks of age, body composition was determined by MRI and a glucose tolerance test was performed after a 4 hr fast. At 20 weeks of age, animals were given an overdose of sodium pentobarbital after a 4 hr fast and blood was collected via cardiac puncture. Euthanasia was achieved by cardiac perfusion with phosphate-buffered saline. After cardiac perfusion, the reproductive fat pad was harvested, flash frozen in liquid nitrogen, and stored at –80 °C.

Study design

Request a detailed protocol

The weaning, phenotyping protocols, and diets were chosen to reproduce the protocols and diets used in studies of the F16 Advanced Intercross of the LG/J x SM/J inbred mouse lines that were used in previously published mapping studies (Cheverud et al., 2011; Lawson et al., 2011a; Lawson et al., 2011b, Lawson et al., 2010). The experimental unit for the current study is the individual mouse. For the RNA sequencing, a single animal was randomly chosen from each cage using a random number generator in R. All other animals served as biological replicates. Mice from multiple cages representing different crosses, generations, diets, and sexes, were necropsied at the same time to avoid batch effects. Library prep and RNA sequencing was performed blinded by the WUSM Genome Technology and Access Center.

Genomes and annotations

Request a detailed protocol

LG/J and SM/J indels and SNVs were leveraged to construct strain-specific genomes using the GRC38.72-mm10 reference as a template (Nikolskiy et al., 2015). This was done by replacing reference bases with alternative (LG/J | SM/J) bases using custom python scripts. Ensembl R72 annotations were adjusted for indel-induced indexing differences for both genomes.

RNA sequencing

Request a detailed protocol

Total RNA was isolated from adipose tissue using the RNeasy Lipid Tissue Kit (QIAgen) (n = 32, 4 animals per sex/diet/cross cohort). RNA concentration was measured via NanoDrop and RNA quality/integrity was assessed with a BioAnalyzer (Agilent). RNA-Seq libraries were constructed using the RiboZero kit (Illumina) from total RNA samples with RIN scores > 8.0. Libraries were checked for quality and concentration using the DNA 1000LabChip assay (Agilent) and quantitative PCR, according to manufacturer’s protocol. Libraries were sequenced at 2 × 100 paired end reads on an Illumina HiSeq 4,000. After sequencing, reads were de-multiplexed and assigned to individual samples. RNAseq data are available through the NCBI-SRA, accession: PRJNA753198.

Library complexity

Request a detailed protocol

Complexity was measured by fitting a beta-binomial distribution to the distribution of L_bias values using the VGAM package (Yee, 2010). The shape parameters (α, β) of beta-binomial distributions were estimated and used to calculate dispersion (ρ). Dispersion values less than 0.05 indicate our libraries are sufficiently complex (Figure 2—figure supplement 1).

ρ_{s} = \frac{1}{1 + α_{s} + β_{s}}

One library was found to have insufficient complexity and was removed from the analyses.

Allele-specific expression

Request a detailed protocol

FASTQ files were filtered to remove low quality reads and aligned against both LG/J and SM/J pseudo-genomes simultaneously using STAR with multimapping disallowed (Dobin et al., 2013). Read counts were normalized via upper quartile normalization and a minimum normalized read depth of 20 was required. Alignment summaries are provided in Supplementary file 8 and Figure 2—figure supplement 2.

For each gene in each individual, allelic bias (L_bias) was calculated as the proportion of total reads for a given gene with the LG/J haplotype. Parent-of-origin-dependent allele-specific expression was detected by ANOVA using one of a number of models in which L_bias is responsive to cross and the interaction of cross with some combination of sex and diet:

M o d e l {\begin{cases} i f e a c h C r o s s c o n t e x t h a s \geq 2 s a m p l e s, L b i a s \sim C r o s s \\ i f e a c h C r o s s : S e x c o n t e x t h a s \geq 2 s a m p l e s, L b i a s \sim C r o s s + C r o s s : S e x \\ i f e a c h C r o s s : D i e t c o n t e x t h a s \geq 2 s a m p l e s, L b i a s \sim C r o s s + C r o s s : D i e t \\ i f e a c h c o n t e x t h a s \geq 2 s a m p l e s, L b i a s \sim C r o s s + C r o s s : S e x + C r o s s : D i e t + C r o s s : S e x : D i e t \end{cases}

Accurately estimating the significance of these effects and correcting for multiple tests is challenging for two reasons: (1) the complexity of the many environmental contexts and (2) the correlation of allelic bias within and between imprinted domains breaks assumptions of independence. A permutation approach is an effective way to overcome these challenges. The context data was randomly shuffled for each gene independently and analyses were rerun to generate a stable null distribution of p-values (Figure 2—figure supplement 3). False discovery rates were estimated for a given significance threshold as the proportion of significant tests under the permutated null model relative to significant tests under the real data model. A value of 1 meaning that 100% of tests at a given significance threshold are likely false positives. An FDR ≤ 0.1 was considered significant (Supplementary file 1, Figure 2—figure supplement 4).

To determine the parental direction and size of expression biases, a parent-of-origin effect score was calculated as the difference in mean L_bias between reciprocal crosses (LxS or SxL). Parent-of-origin effect scores range from completely maternally expressed (–1), to biallelic (0), to completely paternally expressed ( + 1). Parent-of-origin effect score thresholds were calculated from a critical value of α = 0.01, determined from a null distribution created by permutation Genes with significant allele-specific expression and parent-of-origin scores beyond the critical value were considered to have significant parent-of-origin-dependent allele-specific expression (Figure 2—figure supplement 5).

Differential expression

Request a detailed protocol

Differential expression by reciprocal cross was determined by first aligning reads against the LG/J and SM/J genomes simultaneously with multimapping permitted. Reads were normalized by Trimmed mean of M-values (TMM) normalization, which estimates scale factors among samples to allow for differences in RNA composition (Robinson and Oshlack, 2010). A minimum normalized read count of 10 was required. Generalized linear models accounting for diet, sex, and diet-by-sex were fit in EdgeR (Robinson et al., 2010). Differential expression was detected by a likelihood ratio test. Significance was determined for five models for each gene:

1 . E x p r e s s i o n ~ C r o s s

2 . E x p r e s s i o n ~ C r o s s : S e x

3 . E x p r e s s i o n ~ C r o s s : D i e t

4 . E x p r e s s i o n ~ C r o s s : S e x : D i e t

5 . E x p r e s s i o n ~ C r o s s + C r o s s : S e x + C r o s s : D i e t + C r o s s : S e x : D i e t

Multiple test corrections were performed by implementing the ‘qvalue’ R package to estimate false discovery rates (Figure 2—figure supplement 6). Genes with a FDR of ≤0.1 and a $|f o l d c h a n g e| \geq 1.5$ were considered significantly differentially expressed by reciprocal cross (Figure 2—figure supplement 7 and Supplementary file 2).

Gene-gene interactions

Request a detailed protocol

Networks were constructed in each tissue by pairing genes showing parent-of-origin-dependent allele-specific expression with biallelic genes that are differentially expressed by cross. Pairs were predicted by modeling the expression of biallelic genes as a function of parent-of-origin-dependent allele-specific expression, L_bias, sex, diet, and sex-by-diet. The strength of a prediction was measured through model fit, which was estimated as a mean test error with 10-fold cross-validation employed to prevent overfitting. Significance was estimated by likelihood ratio test using a chi-square distribution. Given the complexity of contexts, false discovery rates were determined by permuting the context and expression data to generate a stable null-distribution of p-value (Figure 2—figure supplement 8) Null distribution stability was evaluated by calculating the critical value for alpha = 0.05 at each genome wide iteration. The standard deviation of critical values was calculated after each iteration for the last 5 iterations. Genome-wide shuffling was done 500 times, with 759 independent randomized tests per iteration, meaning the stable null model is composed of 379,500 randomized observations. Using the null model, the ‘qvalue’ package estimated a $\hat{π_{0}}$ . This estimate was then used to estimate false discovery rates in the real data. MTE score thresholds were calculated from a critical value of α = 0.01, determined from a null distribution created by permutation (Figure 2—figure supplement 9). Connections with an FDR ≤ 0.1 (Supplementary file 9) and MTE below the critical value were considered significant (Figure 2—figure supplement 10).

Functional enrichment analysis

Request a detailed protocol

Functional enrichment of interacting genes showing parent-of-origin-dependent allele-specific expression with biallelic genes that are differentially expressed by cross was tested by over-representation analysis in the WEB-based Gene Set Analysis Toolkit v2019 (Zhang et al., 2005). We performed analyses of gene ontologies (biological process, cellular component, and molecular function), pathway (KEGG), and phenotype (Mammalian Phenotype Ontology). The list of all unique interacting genes was analyzed against the background of all unique genes expressed in white adipose. A Benjamini-Hochberg FDR-corrected p-value ≤ 0.01 was considered significant (Supplementary file 4).

Phenotype correlation

Request a detailed protocol

In order to identify which phenotypes might be affected by genes in the parent-of-origin effects network, gene expression was correlated with metabolic phenotypes collected for F₁ animals with the contexts combined. Phenotypes were log transformed when necessary, as determined by Shapiro Wilkes test to approximate normality (Figure 2—figure supplement 11). Additionally, the effects of sex and diet were residualized out leaving only the effect of cross. This was done to mirror later residualizing of phenotypes in the F16 population when testing for epistasis. Multiple test corrections were performed by implementing the ‘qvalue’ R package to estimate false discovery rates (Figure 2—figure supplement 12). The minimum Pearson’s correlation coefficient threshold was set to |0.5|. Connections with an FDR ≤ 0.05 (Supplementary file 5) and MTE below the critical value were considered significant (Figure 2—figure supplement 13).

Epistasis testing

Request a detailed protocol

The F₁₆ LxS advanced intercross population, phenotypes, genotypes, genotypic scores, and QTL mapping methods are described elsewhere (Cheverud et al., 2011; Lawson et al., 2011a; Lawson et al., 2011b, Lawson et al., 2010). We tested for epistasis in interacting pairs between genes showing parent-of-origin-dependent allele-specific expression and biallelic genes that are differentially expressed by cross. We selected F₁₆ genotyped markers that fall within 1.5 Mb up- and downstream from the geometric center of each gene, defined as the genomic position halfway between the transcription start and stop position of that gene (Supplementary file 10). For every F₁₆ animal, an ‘imprinting score’ was assigned to each marker based on that animal’s genotypic values (LL = 0, LS = 1, SL = –1, SS = 0; maternal allele is depicted first). Non-normally distributed phenotypes (as evaluated by a Shapiro-Wilk test) were log₁₀-transformed to approximate normality (Figure 3—figure supplement 1). Because of the number of epistasis tests performed and the number of contexts represented in the data, we removed the effects of sex, diet and their interaction from each F₁₆ phenotype with a covariate screen. We tested for epistasis on the residualized data using the following generalized linear model:

R_{p h e n o} ~ {B D E}_{I M P} + {A S E}_{I M P} + {B D E}_{I M P} : {A S E}_{I M P}

where R_pheno is the residual phenotype, BDE_IMP is the imprinted genotypic score for the biallelic gene that is differentially expressed by cross, ASE_IMP is the imprinted genotypic score for the gene showing parent-of-origin-dependent allele-specific expression bias, and BDE_IMP:ASE_IMP is the interaction between the two genes’ imprinted genotypic score. We employed a permutation approach to accurately estimate significance given the linkage of proximal markers. Imprinted genotypic values were randomly shuffled to generate a stable null model of p-values (Figure 3—figure supplement 2). False discovery rates were estimated for a given significance threshold as the proportion of significant tests under the permutated null model relative to significant tests under the real data model (Figure 3—figure supplement 3). An FDR ≤ 0.1 was considered significant. Epistasis was considered significant if the BDE_IMP: ASE_IMP interaction term met the significance threshold (Supplementary file 6).

Validation of Nnat and F2r expression patterns

Request a detailed protocol

Expression patterns of Nnat and F2r in white adipose were validated by qRT-PCR in high-fat-fed female LG/J and SM/J mice and in biological replicates of high-fat-fed female F₁ reciprocal cross animals (n = 6 LG/J homozygotes, n = 10 LxS and 10 SxL reciprocal heterozygotes, n = 6 SM/J homozygotes). Total RNA was extracted from adipose samples using the Qiagen RNeasy Lipid Kit. High-Capacity cDNA Reverse Transcription Kit (Thermo Fisher) was used for reverse transcription. Quantitative RT-PCR was performed with an Applied Biosystems (USA) QuantStudio 6 Flex instrument using SYBR Green reagent. Results were normalized to L32 expression using the ΔΔCt method. Nnat forward primer – CTACCCCAAGAGCTCCCTTT and reverse primer – CAGCTTCTGCAGGGAGTACC. F2r forward primer – TGAACCCCCGCTCATTCTTTC and reverse primer – CCAGCAGGACGCTTTCATTTTT. L32 forward primer – TCCACAATGTCAAGGAGCTG and reverse primer – GGGATTGGTGACTCTGATGG. Data points were considered outliers if they led to violation of normality assumptions or were considered outliers by box and whisker plots. ANOVA was used to estimate significance of differential expression by cross (1), paternal allele identity (2), mitochondrial ancestry (3).

1. E x p r e s s i o n \sim C r o s s \in {\begin{cases} L L, 0 \\ L S, - 1 \\ S L, 1 \\ S S, 0 \end{cases}

2. E x p r e s s i o n \sim P a t e r n a l A l l e l e \in {\begin{cases} L L, 0 \\ L S, 1 \\ S L, 0 \\ S S, 1 \end{cases}

3. E x p r e s s i o n \sim M i t o c h o n d r i a l a n c e s t r y \in {\begin{cases} L x S x S x L, 0 \\ S x L x L x S, 1 \end{cases}

Expression patterns were also validated by qRT-PCR in high fat-fed female F₂ animals (n = 14). Co-expression was determined by fitting a general linear model and estimating significance using the Wald test approximation of the LR test. Correlation with basal glucose was determined by fitting a general linear model and estimating significance using the Wald test approximation of the LR test. Pearson’s correlation coefficients were calculated for each gene with basal glucose. To test whether patterns in these correlations was significantly different between F₀/F₁ and F₂ populations, bootstrapping was used to calculate 90% confidence intervals for the Pearson’s correlation coefficients. 5,000 iterations were run with 10 individuals randomly selected with replacement. scRNA analysis of Nnat and F2r scRNAseq data was downloaded from SRA: SRP145475 (Burl et al., 2018). Data were processed and aligned to the C57BL/6 J reference (mm10) using Cell Ranger (Zheng et al., 2017). Analysis and cell quality control was performed using the Seurat (3.2.2)(Stuart et al., 2019) package in R (3.6.1)(R Development Core Team, 2013). Cell quality was controlled using three metrics (Luecken and Theis, 2019): (1) number of features, (2) number of counts, (3) covariation of features and counts. High quality cells were required to have between 500 and 3000 features and read counts between 1000 and 30,000. As sequencing is a process of random sampling, the number of features and the number of counts should covary. This relationship was fit to a generalized additive model. Deviation from this relationship (residuals) were computed for each cell. High-quality cells were required to have a residual within 3 standard deviations of the mean residual of all cells (Figure 5—figure supplement 3).

Seurat normalization with a scale factor of 10,000 was performed. Dimensionality reduction (UMAP) was performed (dims = 1:10, resolution = 0.15). Resolution was chosen using the clustree (0.4.3) package (Zappia and Oshlack, 2018). A range of resolutions from 0.06 to 0.18 were tested, and the highest resolution with stable clustering was chosen (Figure 5—figure supplement 4). Cell type markers were identified by differential expression analysis using the ‘MAST’ hurdle-model test (Finak et al., 2015). Genes overexpressed in a given cell type relative to all other cell types were considered cell type ‘markers’. Cell type identity was assigned using previously reported markers for this data set (Figure 5—figure supplement 1).

Cells along the adipogenic trajectory were subset and subjected to dimensionality reduction (UMAP, dims = 1:10, resolution = 0.17). A range of resolutions from 0.01 to 0.25 were tested. Using Adipoq as a marker of differentiation, we sought to identify the set of clusters that would best encapsulate the stages of differentiation. To this end for every level of resolution we calculated the mean count variance ( $\bar{C_{σ}}$ ). This is done by calculating the standard deviation ( $σ$ ) of Adipoq expression (E) within each cluster (G), referred to as the count variance ( $C_{σ}$ ). Cells with no expression of Adipoq were excluded. The mean of count variances for all clusters is calculated. This process is similar to k-means clustering, where the goal is to find that parameters which minimize the within group variation.

C o u n t V a r i a n c e = \frac{\sum_{G = 1}^{n} σ (E_{G})}{n}

We also calculated the percent expressing variance ( $\bar{P_{σ}}$ ). This was taken as the mean of the standard deviation in the percent of cells expressing Adipoq.

P e r c e n t E x p r e s s i n g V a r i a n c e = \frac{\sum_{G = 1}^{n} σ (% E > 0_{G})}{n}

The resolution 0.17 was chosen as the lowest resolution where variation is minimized and no longer significantly changes (Figure 5—figure supplement 5). Using Adipoq as a marker of adipogenesis, clusters 1 and 2 were identified as pre- and post-differentiated cells, respectively. Differential expression was analyzed using the ‘MAST’ test. Expression was compared between clusters 1 and 2 only. Multiple tests correction was performed using the Bonferroni method. We required changes in expression to show either a sufficiently large fold change ( $|\log_{2} F o l d C h a n g e| \geq 0.3$ ) OR a sufficiently large change in the percent of cells expressing the gene in question ( $p c t . ∆ \geq 0.4$ ). The change in percent of cells expressing a gene was calculated as the difference in percent of cells expressing the gene between the clusters and scaled by dividing by the larger percentage.

p c t . ∆ = \frac{p c t . 2 - p c t . 1}{m a x (p c t . 1, p c t . 2)}

Source code is available at https://github.com/LawsonLab-WUSM/POE_Epistasis, (copy archived at swh:1:rev:b39046ce35f53e0c3f15bcdefa122c274aee48b7, Lawson, 2019).

Data availability

Sequencing data are available through the NCBI-SRA under accession code PRJNA753198.

The following data sets were generated

1. Lawson HA
(2021) NCBI BioProject
ID PRJNA753198. LG/J and SM/J hybrid RNAseq.

https://www.ncbi.nlm.nih.gov/bioproject/PRJNA753198/

The following previously published data sets were used

1. Burl RB
(2018) European Nucleotide Archive
ID PRJNA470640. Single-cell RNA-sequencing of white adipose tissue stromal cells during CL-induced adipogenesis.

https://www.ebi.ac.uk/ena/browser/view/PRJNA470640?show=reads

References

1. Akune T
2. Ohba S
3. Kamekura S
4. Yamaguchi M
5. Chung U-I
6. Kubota N
7. Terauchi Y
8. Harada Y
9. Azuma Y
10. Nakamura K
11. Kadowaki T
12. Kawaguchi H
(2004) PPARgamma insufficiency enhances osteogenesis through osteoblast formation from bone marrow progenitors
The Journal of Clinical Investigation 113:846–855.

https://doi.org/10.1172/JCI19900
- PubMed
- Google Scholar
(2005) Rosiglitazone causes bone loss in mice by suppressing osteoblast differentiation and bone formation
Endocrinology 146:1226–1235.

https://doi.org/10.1210/en.2004-0735
- PubMed
- Google Scholar
1. An JJ
2. Han DH
3. Kim DM
4. Kim SH
5. Rhee Y
6. Lee EJ
7. Lim SK
(2007) Expression and regulation of osteoprotegerin in adipose tissue
Yonsei Medical Journal 48:765–772.

https://doi.org/10.3349/ymj.2007.48.5.765
- PubMed
- Google Scholar
1. Aubin JE
(1992) Osteoclast adhesion and resorption: the role of podosomes
Journal of Bone and Mineral Research 7:365–368.

https://doi.org/10.1002/jbmr.5650070402
- PubMed
- Google Scholar
1. Babak T
2. DeVeale B
3. Tsang EK
4. Zhou Y
5. Li X
6. Smith KS
7. Kukurba KR
8. Zhang R
9. Li JB
10. van der Kooy D
11. Montgomery SB
12. Fraser HB
(2015) Genetic conflict reflected in tissue-specific maps of genomic imprinting in human and mouse
Nature Genetics 47:544–549.

https://doi.org/10.1038/ng.3274
- PubMed
- Google Scholar
1. Berridge MJ
(2016) The Inositol Trisphosphate/Calcium Signaling Pathway in Health and Disease
Physiological Reviews 96:1261–1296.

https://doi.org/10.1152/physrev.00006.2016
- PubMed
- Google Scholar
(2010) Superior osteogenic capacity of human embryonic stem cells adapted to matrix-free growth compared to human mesenchymal stem cells
Tissue Engineering. Part A 16:3427–3440.

https://doi.org/10.1089/ten.tea.2010.0112
- PubMed
- Google Scholar
1. Burg JS
2. Espenshade PJ
(2011) Regulation of HMG-CoA reductase in mammals and yeast
Progress in Lipid Research 50:403–410.

https://doi.org/10.1016/j.plipres.2011.07.002
- PubMed
- Google Scholar
(2018) Deconstructing Adipogenesis Induced by β3-Adrenergic Receptor Activation with Single-Cell Expression Profiling
Cell Metabolism 28:300–309.

https://doi.org/10.1016/j.cmet.2018.05.025
- PubMed
- Google Scholar
(2020) Brown Adipose Expansion and Remission of Glycemic Dysfunction in Obese SM/J Mice
Cell Reports 33:108237.

https://doi.org/10.1016/j.celrep.2020.108237
- PubMed
- Google Scholar
1. Cheverud JM
2. Lawson HA
3. Fawcett GL
4. Wang B
5. Pletscher LS
6. R Fox A
7. Maxwell TJ
8. Ehrich TH
9. Kenney-Hunt JP
10. Wolf JB
11. Semenkovich CF
(2011) Diet-dependent genetic and genomic imprinting effects on obesity in mice
Obesity (Silver Spring, Md.) 19:160–170.

https://doi.org/10.1038/oby.2010.141
- PubMed
- Google Scholar
1. Cho SW
2. Yang J-Y
3. Her SJ
4. Choi HJ
5. Jung JY
6. Sun HJ
7. An JH
8. Cho HY
9. Kim SW
10. Park KS
11. Kim SY
12. Baek W-Y
13. Kim J-E
14. Yim M
15. Shin CS
(2011) Osteoblast-targeted overexpression of PPARγ inhibited bone mass gain in male mice and accelerated ovariectomy-induced bone loss in female mice
Journal of Bone and Mineral Research 26:1939–1952.

https://doi.org/10.1002/jbmr.366
- PubMed
- Google Scholar
1. Dobin A
2. Davis CA
3. Schlesinger F
4. Drenkow J
5. Zaleski C
6. Jha S
7. Batut P
8. Chaisson M
9. Gingeras TR
(2013) STAR: ultrafast universal RNA-seq aligner
Bioinformatics (Oxford, England) 29:15–21.

https://doi.org/10.1093/bioinformatics/bts635
- PubMed
- Google Scholar
(2016) Association of Osteoprotegerin with Obesity, Insulin Resistance and Non-Alcoholic Fatty Liver Disease in Children
Iranian Red Crescent Medical Journal 18:41873.

https://doi.org/10.5812/ircmj.41873
- PubMed
- Google Scholar
1. Finak G
2. McDavid A
3. Yajima M
4. Deng J
5. Gersuk V
6. Shalek AK
7. Slichter CK
8. Miller HW
9. McElrath MJ
10. Prlic M
11. Linsley PS
12. Gottardo R
(2015) MAST: a flexible statistical framework for assessing transcriptional changes and characterizing heterogeneity in single-cell RNA sequencing data
Genome Biology 16:1–13.

https://doi.org/10.1186/s13059-015-0844-5
- PubMed
- Google Scholar
(2017) Integrative analysis of pathway deregulation in obesity
NPJ Systems Biology and Applications 3:18.

https://doi.org/10.1038/s41540-017-0018-z
- PubMed
- Google Scholar
1. Fowden AL
2. Moore T
(2012) Maternal-fetal resource allocation: co-operation and conflict
Placenta 33 Suppl 2:e11–e15.

https://doi.org/10.1016/j.placenta.2012.05.002
- PubMed
- Google Scholar
1. Franck N
2. Gummesson A
3. Jernås M
4. Glad C
5. Svensson PA
6. Guillot G
7. Rudemo M
8. Nyström FH
9. Carlsson LMS
10. Olsson B
(2011) Identification of adipocyte genes regulated by caloric intake
The Journal of Clinical Endocrinology and Metabolism 96:E413–E418.

https://doi.org/10.1210/jc.2009-2534
- PubMed
- Google Scholar
(2008) Maternal effects as the cause of parent-of-origin effects that mimic genomic imprinting
Genetics 178:1755–1762.

https://doi.org/10.1534/genetics.107.080697
- PubMed
- Google Scholar
1. Haig D
(1997) Parental antagonism, relatedness asymmetries, and genomic imprinting
Proceedings. Biological Sciences 264:1657–1662.

https://doi.org/10.1098/rspb.1997.0230
- PubMed
- Google Scholar
(2007) The influence of weight loss on serum osteoprotegerin concentration in obese perimenopausal women
Obesity (Silver Spring, Md.) 15:1925–1929.

https://doi.org/10.1038/oby.2007.229
- PubMed
- Google Scholar
1. Ikeda K
2. Yamada T
(2020) UCP1 Dependent and Independent Thermogenesis in Brown and Beige Adipocytes
Frontiers in Endocrinology 11:498.

https://doi.org/10.3389/fendo.2020.00498
- PubMed
- Google Scholar
(2013) Serglycin is a novel adipocytokine highly expressed in epicardial adipose tissue
Biochemical and Biophysical Research Communications 432:105–110.

https://doi.org/10.1016/j.bbrc.2013.01.078
- PubMed
- Google Scholar
Software
1. Jirtle RL
(2012) No Title
Geneimprint.

http://www.geneimprint.com/site/genes-by-species
1. Jo Y
2. Debose-Boyd RA
(2010) Control of cholesterol synthesis through regulated ER-associated degradation of HMG CoA reductase
Critical Reviews in Biochemistry and Molecular Biology 45:185–198.

https://doi.org/10.3109/10409238.2010.485605
- PubMed
- Google Scholar
1. Kang JW
2. Choi Y
3. Park JH
4. Kim JS
5. Park KD
6. Baek DH
7. Seong SK
8. Choi KS
9. Lim SY
10. Kim HS
(2008) The effects of cyclin-dependent kinase inhibitors on adipogenic differentiation of human mesenchymal stem cells
Biochemical and Biophysical Research Communications 366:624–630.

https://doi.org/10.1016/j.bbrc.2007.11.080
- PubMed
- Google Scholar
1. Kichaev G
2. Bhatia G
3. Loh PR
4. Gazal S
5. Burch K
6. Freund MK
7. Schoech A
8. Pasaniuc B
9. Price AL
(2019) Leveraging Polygenic Functional Enrichment to Improve GWAS Power
American Journal of Human Genetics 104:65–75.

https://doi.org/10.1016/j.ajhg.2018.11.008
- PubMed
- Google Scholar
(2009) Parental origin of sequence variants associated with complex diseases
Nature 462:868–874.

https://doi.org/10.1038/nature08625
- PubMed
- Google Scholar
1. Lawson HA
2. Zelle KM
3. Fawcett GL
4. Wang B
5. Pletscher LS
6. Maxwell TJ
7. Ehrich TH
8. Kenney-Hunt JP
9. Wolf JB
10. Semenkovich CF
11. Cheverud JM
(2010) Genetic, epigenetic, and gene-by-diet interaction effects underlie variation in serum lipids in a LG/JxSM/J murine model
Journal of Lipid Research 51:2976–2984.

https://doi.org/10.1194/jlr.M006957
- PubMed
- Google Scholar
(2011a) Genetic effects at pleiotropic loci are context-dependent with consequences for the maintenance of genetic variation in populations
PLOS Genetics 7:e1002256.

https://doi.org/10.1371/journal.pgen.1002256
- PubMed
- Google Scholar
1. Lawson HA
2. Lee A
3. Fawcett GL
4. Wang B
5. Pletscher LS
6. Maxwell TJ
7. Ehrich TH
8. Kenney-Hunt JP
9. Wolf JB
10. Semenkovich CF
11. Cheverud JM
(2011b) The importance of context to the genetic architecture of diabetes-related traits is revealed in a genome-wide scan of a LG/J × SM/J murine model
Mammalian Genome 22:197–208.

https://doi.org/10.1007/s00335-010-9313-3
- PubMed
- Google Scholar
(2013) Genomic imprinting and parent-of-origin effects on complex traits
Nature Reviews. Genetics 14:609–617.

https://doi.org/10.1038/nrg3543
- PubMed
- Google Scholar
Software
1. Lawson HA
(2019) LawsonLab-WUSM/POE_Epistasis, version swh:1:rev:b39046ce35f53e0c3f15bcdefa122c274aee48b7
Software Heritage.

https://archive.softwareheritage.org/swh:1:dir:3e0c4952b5f93241102e0b8a1c030ac0ff41ce93;origin=https://github.com/LawsonLab-WUSM/POE_Epistasis;visit=swh:1:snp:532c07d4bba72b337aa77fcc9cd90553ecde1426;anchor=swh:1:rev:b39046ce35f53e0c3f15bcdefa122c274aee48b7
1. Lemire JM
2. Chan CK
3. Bressler S
4. Miller J
5. LeBaron RG
6. Wight TN
(2007) Interleukin-1beta selectively decreases the synthesis of versican by arterial smooth muscle cells
Journal of Cellular Biochemistry 101:753–766.

https://doi.org/10.1002/jcb.21235
- PubMed
- Google Scholar
1. Li X
2. Thomason PA
3. Withers DJ
4. Scott J
(2010) Bio-informatics analysis of a gene co-expression module in adipose tissue containing the diet-responsive gene Nnat
BMC Systems Biology 4:1–11.

https://doi.org/10.1186/1752-0509-4-175
- PubMed
- Google Scholar
1. Lin HH
2. Bell E
3. Uwanogho D
4. Perfect LW
5. Noristani H
6. Bates TJD
7. Snetkov V
8. Price J
9. Sun YM
(2010) Neuronatin promotes neural lineage in ESCs via Ca(2+) signaling
Stem Cells (Dayton, Ohio) 28:1950–1960.

https://doi.org/10.1002/stem.530
- PubMed
- Google Scholar
1. Luecken MD
2. Theis FJ
(2019) Current best practices in single-cell RNA-seq analysis: a tutorial
Molecular Systems Biology 15:e8746.

https://doi.org/10.15252/msb.20188746
- PubMed
- Google Scholar
1. Lutz SM
2. Hokanson JE
(2015) Mediation analysis in genome-wide association studies: current perspectives
Open Access Bioinformatics 7:S63643.

https://doi.org/10.2147/OAB.S63643
- Google Scholar
(2020) RANKL induces beige adipocyte differentiation in preadipocytes
American Journal of Physiology. Endocrinology and Metabolism 318:E866–E877.

https://doi.org/10.1152/ajpendo.00397.2019
- PubMed
- Google Scholar
1. McDonald LJ
2. Mamrack MD
(1995) Phosphoinositide hydrolysis by phospholipase C modulated by multivalent cations La(3+), Al(3+), neomycin, polyamines, and melittin
Journal of Lipid Mediators and Cell Signalling 11:81–91.

https://doi.org/10.1016/0929-7855(94)00029-c
- PubMed
- Google Scholar
1. Millership SJ
2. Da Silva Xavier G
3. Choudhury AI
4. Bertazzo S
5. Chabosseau P
6. Pedroni SM
7. Irvine EE
8. Montoya A
9. Faull P
10. Taylor WR
11. Kerr-Conte J
12. Pattou F
13. Ferrer J
14. Christian M
15. John RM
16. Latreille M
17. Liu M
18. Rutter GA
19. Scott J
20. Withers DJ
(2018) Neuronatin regulates pancreatic β cell insulin content and secretion
The Journal of Clinical Investigation 128:3369–3381.

https://doi.org/10.1172/JCI120115
- PubMed
- Google Scholar
(2020) Spontaneous restoration of functional β‐cell mass in obese SM/J mice
Physiological Reports 8:1–12.

https://doi.org/10.14814/phy2.14573
- Google Scholar
(1996) The evolution of genomic imprinting
Genetics 144:1283–1295.

https://doi.org/10.1093/genetics/144.3.1283
- PubMed
- Google Scholar
1. Mott R
2. Yuan W
3. Kaisaki P
4. Gan X
5. Cleak J
6. Edwards A
7. Baud A
8. Flint J
(2014) The architecture of parent-of-origin effects in mice
Cell 156:332–342.

https://doi.org/10.1016/j.cell.2013.11.043
- PubMed
- Google Scholar
1. Mozaffari SV
2. DeCara JM
3. Shah SJ
4. Sidore C
5. Fiorillo E
6. Cucca F
7. Lang RM
8. Nicolae DL
9. Ober C
(2019) Parent-of-origin effects on quantitative phenotypes in a large Hutterite pedigree
Communications Biology 2:28.

https://doi.org/10.1038/s42003-018-0267-4
- Google Scholar
1. Nakagawa N
2. Kinosaki M
3. Yamaguchi K
4. Shima N
5. Yasuda H
6. Yano K
7. Morinaga T
8. Higashio K
(1998) RANK is the essential signaling receptor for osteoclast differentiation factor in osteoclastogenesis
Biochemical and Biophysical Research Communications 253:395–400.

https://doi.org/10.1006/bbrc.1998.9788
- PubMed
- Google Scholar
1. Nakamura Y
2. Hamada Y
3. Fujiwara T
4. Enomoto H
5. Hiroe T
6. Tanaka S
7. Nose M
8. Nakahara M
9. Yoshida N
10. Takenawa T
11. Fukami K
(2005) Phospholipase C-delta1 and -delta3 are essential in the trophoblast for placental development
Molecular and Cellular Biology 25:10979–10988.

https://doi.org/10.1128/MCB.25.24.10979-10988.2005
- PubMed
- Google Scholar
1. Nikolskiy I
2. Conrad DF
3. Chun S
4. Fay JC
5. Cheverud JM
6. Lawson HA
(2015) Using whole-genome sequences of the LG/J and SM/J inbred mouse strains to prioritize quantitative trait genes and nucleotides
BMC Genomics 16:415.

https://doi.org/10.1186/s12864-015-1592-3
- PubMed
- Google Scholar
1. Parker-Katiraee L
2. Carson AR
3. Yamada T
4. Arnaud P
5. Feil R
6. Abu-Amero SN
7. Moore GE
8. Kaneda M
9. Perry GH
10. Stone AC
11. Lee C
12. Meguro-Horike M
13. Sasaki H
14. Kobayashi K
15. Nakabayashi K
16. Scherer SW
(2007) Identification of the imprinted KLF14 transcription factor undergoing human-specific accelerated evolution
PLOS Genetics 3:0030065.

https://doi.org/10.1371/journal.pgen.0030065
- PubMed
- Google Scholar
1. Patten MM
2. Ross L
3. Curley JP
4. Queller DC
5. Bonduriansky R
6. Wolf JB
(2014) The evolution of genomic imprinting: theories, predictions and empirical tests
Heredity 113:119–128.

https://doi.org/10.1038/hdy.2014.29
- PubMed
- Google Scholar
(2003) Inhibition of calcium uptake via the sarco/endoplasmic reticulum Ca2+-ATPase in a mouse model of Sandhoff disease and prevention by treatment with N-butyldeoxynojirimycin
The Journal of Biological Chemistry 278:29496–29501.

https://doi.org/10.1074/jbc.M302964200
- PubMed
- Google Scholar
1. Pierce BL
2. Tong L
3. Chen LS
4. Rahaman R
5. Argos M
6. Jasmine F
7. Roy S
8. Paul-Brutus R
9. Westra HJ
10. Franke L
11. Esko T
12. Zaman R
13. Islam T
14. Rahman M
15. Baron JA
16. Kibriya MG
17. Ahsan H
(2014) Mediation analysis demonstrates that trans-eQTLs are often explained by cis-mediation: A genome-wide analysis among 1,800 South Asians
PLOS Genetics 10:e1004818.

https://doi.org/10.1371/journal.pgen.1004818
- PubMed
- Google Scholar
1. Prockop DJ
(1997) Marrow stromal cells as stem cells for nonhematopoietic tissues
Science (New York, N.Y.) 276:71–74.

https://doi.org/10.1126/science.276.5309.71
- PubMed
- Google Scholar
Software
1. R Development Core Team
(2013) R: A Language and Environment for Statistical Computing
R Foundation for Statistical Computing, Vienna, Austria.

http://www.r-project.org
(2017) Carbonic anhydrase III (Car3) is not required for fatty acid synthesis and does not protect against high-fat diet induced obesity in mice
PLOS ONE 12:e0176502.

https://doi.org/10.1371/journal.pone.0176502
- PubMed
- Google Scholar
(2010) edgeR: a Bioconductor package for differential expression analysis of digital gene expression data
Bioinformatics (Oxford, England) 26:139–140.

https://doi.org/10.1093/bioinformatics/btp616
- PubMed
- Google Scholar
1. Robinson MD
2. Oshlack A
(2010) A scaling normalization method for differential expression analysis of RNA-seq data
Genome Biology 11:1–9.

https://doi.org/10.1186/gb-2010-11-3-r25
- PubMed
- Google Scholar
(2008) Involvement of adipogenic potential of human bone marrow mesenchymal stem cells (MSCs) in osteoporosis
Current Stem Cell Research & Therapy 3:208–218.

https://doi.org/10.2174/157488808785740325
- PubMed
- Google Scholar
1. Rosen CJ
2. Bouxsein ML
(2006) Mechanisms of disease: is osteoporosis the obesity of bone?
Nature Clinical Practice. Rheumatology 2:35–43.

https://doi.org/10.1038/ncprheum0070
- PubMed
- Google Scholar
(2001) Synthesis, secretion, and subcellular localization of serglycin proteoglycan in human endothelial cells
Blood 97:449–458.

https://doi.org/10.1182/blood.v97.2.449
- PubMed
- Google Scholar
1. Shan N
2. Wang Z
3. Hou L
(2019) Identification of trans-eQTLs using mediation analysis with multiple mediators
BMC Bioinformatics 20:14–16.

https://doi.org/10.1186/s12859-019-2651-6
- PubMed
- Google Scholar
(2011) Identification of an imprinted master trans regulator at the KLF14 locus related to multiple metabolic phenotypes
Nature Genetics 43:561–564.

https://doi.org/10.1038/ng.833
- PubMed
- Google Scholar
1. Soh UJK
2. Dores MR
3. Chen B
4. Trejo J
(2010) Signal transduction by protease-activated receptors
British Journal of Pharmacology 160:191–203.

https://doi.org/10.1111/j.1476-5381.2010.00705.x
- PubMed
- Google Scholar
(1991) Expression of CA III in rodent models of obesity
Molecular Endocrinology (Baltimore, Md.) 5:860–866.

https://doi.org/10.1210/mend-5-6-860
- PubMed
- Google Scholar
1. Stuart T
2. Butler A
3. Hoffman P
4. Hafemeister C
5. Papalexi E
6. Mauck WM
7. Hao Y
8. Stoeckius M
9. Smibert P
10. Satija R
(2019) Comprehensive Integration of Single-Cell Data
Cell 177:1888–1902.

https://doi.org/10.1016/j.cell.2019.05.031
- Google Scholar
1. Teslovich TM
2. Musunuru K
3. Smith AV
4. Edmondson AC
5. Stylianou IM
6. Koseki M
7. Pirruccello JP
8. Ripatti S
9. Chasman DI
10. Willer CJ
11. Johansen CT
12. Fouchier SW
13. Isaacs A
14. Peloso GM
15. Barbalic M
16. Ricketts SL
17. Bis JC
18. Aulchenko YS
19. Thorleifsson G
20. Feitosa MF
21. Chambers J
22. Orho-Melander M
23. Johnson T
24. Li X
25. Guo X
26. Li M
27. Shin Cho Y
28. Jin Go M
29. Jin Kim Y
30. Lee JY
31. Park T
32. Kim K
33. Sim X
34. Twee-Hee Ong R
35. Croteau-Chonka DC
36. Lange LA
37. Smith JD
38. Song K
39. Hua Zhao J
40. Yuan X
41. Luan J
42. Lamina C
43. Ziegler A
44. Zhang W
45. Zee RYL
46. Wright AF
47. Witteman JCM
48. Wilson JF
49. Willemsen G
50. Wichmann HE
51. Whitfield JB
52. Waterworth DM
53. Wareham NJ
54. Waeber G
55. Vollenweider P
56. Voight BF
57. Vitart V
58. Uitterlinden AG
59. Uda M
60. Tuomilehto J
61. Thompson JR
62. Tanaka T
63. Surakka I
64. Stringham HM
65. Spector TD
66. Soranzo N
67. Smit JH
68. Sinisalo J
69. Silander K
70. Sijbrands EJG
71. Scuteri A
72. Scott J
73. Schlessinger D
74. Sanna S
75. Salomaa V
76. Saharinen J
77. Sabatti C
78. Ruokonen A
79. Rudan I
80. Rose LM
81. Roberts R
82. Rieder M
83. Psaty BM
84. Pramstaller PP
85. Pichler I
86. Perola M
87. Penninx B
88. Pedersen NL
89. Pattaro C
90. Parker AN
91. Pare G
92. Oostra BA
93. O’Donnell CJ
94. Nieminen MS
95. Nickerson DA
96. Montgomery GW
97. Meitinger T
98. McPherson R
99. McCarthy MI
100. McArdle W
101. Masson D
102. Martin NG
103. Marroni F
104. Mangino M
105. Magnusson PKE
106. Lucas G
107. Luben R
108. Loos RJF
109. Lokki ML
110. Lettre G
111. Langenberg C
112. Launer LJ
113. Lakatta EG
114. Laaksonen R
115. Kyvik KO
116. Kronenberg F
117. König IR
118. Khaw KT
119. Kaprio J
120. Kaplan LM
121. Johansson A
122. Jarvelin MR
123. Janssens A
124. Ingelsson E
125. Igl W
126. Kees Hovingh G
127. Hottenga JJ
128. Hofman A
129. Hicks AA
130. Hengstenberg C
131. Heid IM
132. Hayward C
133. Havulinna AS
134. Hastie ND
135. Harris TB
136. Haritunians T
137. Hall AS
138. Gyllensten U
139. Guiducci C
140. Groop LC
141. Gonzalez E
142. Gieger C
143. Freimer NB
144. Ferrucci L
145. Erdmann J
146. Elliott P
147. Ejebe KG
148. Döring A
149. Dominiczak AF
150. Demissie S
151. Deloukas P
152. de Geus EJC
153. de Faire U
154. Crawford G
155. Collins FS
156. Chen YI
157. Caulfield MJ
158. Campbell H
159. Burtt NP
160. Bonnycastle LL
161. Boomsma DI
162. Boekholdt SM
163. Bergman RN
164. Barroso I
165. Bandinelli S
166. Ballantyne CM
167. Assimes TL
168. Quertermous T
169. Altshuler D
170. Seielstad M
171. Wong TY
172. Tai ES
173. Feranil AB
174. Kuzawa CW
175. Adair LS
176. Taylor HA
177. Borecki IB
178. Gabriel SB
179. Wilson JG
180. Holm H
181. Thorsteinsdottir U
182. Gudnason V
183. Krauss RM
184. Mohlke KL
185. Ordovas JM
186. Munroe PB
187. Kooner JS
188. Tall AR
189. Hegele RA
190. Kastelein JJP
191. Schadt EE
192. Rotter JI
193. Boerwinkle E
194. Strachan DP
195. Mooser V
196. Stefansson K
197. Reilly MP
198. Samani NJ
199. Schunkert H
200. Cupples LA
201. Sandhu MS
202. Ridker PM
203. Rader DJ
204. van Duijn CM
205. Peltonen L
206. Abecasis GR
207. Boehnke M
208. Kathiresan S
(2010) Biological, clinical and population relevance of 95 loci for blood lipids
Nature 466:707–713.

https://doi.org/10.1038/nature09270
- PubMed
- Google Scholar
1. Thatcher JD
(2010) The inositol trisphosphate (IP3) signal transduction pathway
Science Signaling 3:tr3.

https://doi.org/10.1126/scisignal.3119tr3
- PubMed
- Google Scholar
(2018) Metabolic Coupling Between Bone Marrow Adipose Tissue and Hematopoiesis
Current Osteoporosis Reports 16:95–104.

https://doi.org/10.1007/s11914-018-0422-3
- PubMed
- Google Scholar
(2018) Loss of Imprinting of Cdkn1c Protects against Age and Diet-Induced Obesity
International Journal of Molecular Sciences 19:E2734.

https://doi.org/10.3390/ijms19092734
- PubMed
- Google Scholar
1. Voight BF
2. Scott LJ
3. Steinthorsdottir V
4. Morris AP
5. Dina C
6. Welch RP
7. Zeggini E
8. Huth C
9. Aulchenko YS
10. Thorleifsson G
11. McCulloch LJ
12. Ferreira T
13. Grallert H
14. Amin N
15. Wu G
16. Willer CJ
17. Raychaudhuri S
18. McCarroll SA
19. Langenberg C
20. Hofmann OM
21. Dupuis J
22. Qi L
23. Segrè AV
24. van Hoek M
25. Navarro P
26. Ardlie K
27. Balkau B
28. Benediktsson R
29. Bennett AJ
30. Blagieva R
31. Boerwinkle E
32. Bonnycastle LL
33. Bengtsson Boström K
34. Bravenboer B
35. Bumpstead S
36. Burtt NP
37. Charpentier G
38. Chines PS
39. Cornelis M
40. Couper DJ
41. Crawford G
42. Doney ASF
43. Elliott KS
44. Elliott AL
45. Erdos MR
46. Fox CS
47. Franklin CS
48. Ganser M
49. Gieger C
50. Grarup N
51. Green T
52. Griffin S
53. Groves CJ
54. Guiducci C
55. Hadjadj S
56. Hassanali N
57. Herder C
58. Isomaa B
59. Jackson AU
60. Johnson PRV
61. Jørgensen T
62. Kao WHL
63. Klopp N
64. Kong A
65. Kraft P
66. Kuusisto J
67. Lauritzen T
68. Li M
69. Lieverse A
70. Lindgren CM
71. Lyssenko V
72. Marre M
73. Meitinger T
74. Midthjell K
75. Morken MA
76. Narisu N
77. Nilsson P
78. Owen KR
79. Payne F
80. Perry JRB
81. Petersen A-K
82. Platou C
83. Proença C
84. Prokopenko I
85. Rathmann W
86. Rayner NW
87. Robertson NR
88. Rocheleau G
89. Roden M
90. Sampson MJ
91. Saxena R
92. Shields BM
93. Shrader P
94. Sigurdsson G
95. Sparsø T
96. Strassburger K
97. Stringham HM
98. Sun Q
99. Swift AJ
100. Thorand B
101. Tichet J
102. Tuomi T
103. van Dam RM
104. van Haeften TW
105. van Herpt T
106. van Vliet-Ostaptchouk JV
107. Walters GB
108. Weedon MN
109. Wijmenga C
110. Witteman J
111. Bergman RN
112. Cauchi S
113. Collins FS
114. Gloyn AL
115. Gyllensten U
116. Hansen T
117. Hide WA
118. Hitman GA
119. Hofman A
120. Hunter DJ
121. Hveem K
122. Laakso M
123. Mohlke KL
124. Morris AD
125. Palmer CNA
126. Pramstaller PP
127. Rudan I
128. Sijbrands E
129. Stein LD
130. Tuomilehto J
131. Uitterlinden A
132. Walker M
133. Wareham NJ
134. Watanabe RM
135. Abecasis GR
136. Boehm BO
137. Campbell H
138. Daly MJ
139. Hattersley AT
140. Hu FB
141. Meigs JB
142. Pankow JS
143. Pedersen O
144. Wichmann H-E
145. Barroso I
146. Florez JC
147. Frayling TM
148. Groop L
149. Sladek R
150. Thorsteinsdottir U
151. Wilson JF
152. Illig T
153. Froguel P
154. van Duijn CM
155. Stefansson K
156. Altshuler D
157. Boehnke M
158. McCarthy MI
159. MAGIC investigators
160. GIANT Consortium
(2010) Twelve type 2 diabetes susceptibility loci identified through large-scale association analysis
Nature Genetics 42:579–589.

https://doi.org/10.1038/ng.609
- PubMed
- Google Scholar
1. Wilkins JF
2. Haig D
(2003) What good is genomic imprinting: the function of parent-specific gene expression
Nature Reviews. Genetics 4:359–368.

https://doi.org/10.1038/nrg1062
- PubMed
- Google Scholar
1. Wolf JB
2. Cheverud JM
(2009) A framework for detecting and characterizing genetic background-dependent imprinting effects
Mammalian Genome 20:681–698.

https://doi.org/10.1007/s00335-009-9209-2
- PubMed
- Google Scholar
1. Yee TW
(2010) The VGAM Package for Categorical Data Analysis
Journal of Statistical Software 32:i10.

https://doi.org/10.18637/jss.v032.i10
- Google Scholar
1. Yeh Y-S
2. Jheng H-F
3. Iwase M
4. Kim M
5. Mohri S
6. Kwon J
7. Kawarasaki S
8. Li Y
9. Takahashi H
10. Ara T
11. Nomura W
12. Kawada T
13. Goto T
(2018) The Mevalonate Pathway Is Indispensable for Adipocyte Survival
IScience 9:175–191.

https://doi.org/10.1016/j.isci.2018.10.019
- PubMed
- Google Scholar
1. Young HS
2. Won HK
3. Moon C
4. Yun HH
5. Eun SY
6. Joo HL
7. Joo SC
8. Song J
9. Jung MH
(2005) Ectopic expression of Neuronatin potentiates adipogenesis through enhanced phosphorylation of cAMP-response element-binding protein in 3T3-L1 cells
Biochemical and Biophysical Research Communications 337:481–489.

https://doi.org/10.1016/j.bbrc.2005.09.078
- Google Scholar
Software
(2022)
Circulating Osteoprotegerin Level in Relation to Obesity in Middle Aged Females

Sapub.
1. Zappia L
2. Oshlack A
(2018) Clustering trees: a visualization for evaluating clusterings at multiple resolutions
GigaScience 7:1–9.

https://doi.org/10.1093/gigascience/giy083
- PubMed
- Google Scholar
1. Zeng Y
2. Amador C
3. Xia C
4. Marioni R
5. Sproul D
6. Walker RM
7. Morris SW
8. Bretherick A
9. Canela-Xandri O
10. Boutin TS
11. Clark DW
12. Campbell A
13. Rawlik K
14. Hayward C
15. Nagy R
16. Tenesa A
17. Porteous DJ
18. Wilson JF
19. Deary IJ
20. Evans KL
21. McIntosh AM
22. Navarro P
23. Haley CS
(2019) Parent of origin genetic effects on methylation in humans are common and influence complex trait variation
Nature Communications 10:1383.

https://doi.org/10.1038/s41467-019-09301-y
- PubMed
- Google Scholar
1. Zernichow L
2. Åbrink M
3. Hallgren J
4. Grujic M
5. Pejler G
6. Kolset SO
(2006) Serglycin Is the Major Secreted Proteoglycan in Macrophages and Has a Role in the Regulation of Macrophage Tumor Necrosis Factor-α Secretion in Response to Lipopolysaccharide
Journal of Biological Chemistry 281:26792–26801.

https://doi.org/10.1074/jbc.M512889200
- PubMed
- Google Scholar
1. Zhang B
2. Kirov S
3. Snoddy J
(2005) WebGestalt: an integrated system for exploring gene sets in various biological contexts
Nucleic Acids Research 33:W741–W748.

https://doi.org/10.1093/nar/gki475
- PubMed
- Google Scholar
1. Zheng GXY
2. Terry JM
3. Belgrader P
4. Ryvkin P
5. Bent ZW
6. Wilson R
7. Ziraldo SB
8. Wheeler TD
9. McDermott GP
10. Zhu J
11. Gregory MT
12. Shuga J
13. Montesclaros L
14. Underwood JG
15. Masquelier DA
16. Nishimura SY
17. Schnall-Levin M
18. Wyatt PW
19. Hindson CM
20. Bharadwaj R
21. Wong A
22. Ness KD
23. Beppu LW
24. Deeg HJ
25. McFarland C
26. Loeb KR
27. Valente WJ
28. Ericson NG
29. Stevens EA
30. Radich JP
31. Mikkelsen TS
32. Hindson BJ
33. Bielas JH
(2017) Massively parallel digital transcriptional profiling of single cells
Nature Communications 8:1–12.

https://doi.org/10.1038/ncomms14049
- PubMed
- Google Scholar

Article and author information

Author details

Juan F Macias-Velasco

Department of Genetics, Washington University School of Medicine, Saint Louis, United States

Contribution
Formal analysis, Methodology, Visualization, Writing - original draft, Writing – review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-2827-4647
Celine L St Pierre

Department of Genetics, Washington University School of Medicine, Saint Louis, United States

Contribution
Visualization, Writing – review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-5465-6601
Jessica P Wayhart

Department of Genetics, Washington University School of Medicine, Saint Louis, United States

Contribution
Investigation

Competing interests
No competing interests declared
Li Yin

Department of Medicine, Washington University School of Medicine, Saint Louis, United States

Contribution
Investigation

Competing interests
No competing interests declared
Larry Spears

Department of Medicine, Washington University School of Medicine, Saint Louis, United States

Contribution
Investigation

Competing interests
No competing interests declared
Mario A Miranda

Department of Genetics, Washington University School of Medicine, Saint Louis, United States

Contribution
Validation

Competing interests
No competing interests declared
Caryn Carson

Department of Genetics, Washington University School of Medicine, Saint Louis, United States

Contribution
Investigation

Competing interests
No competing interests declared
Katsuhiko Funai

Diabetes and Metabolism Research Center, University of Utah, Salt Lake City, United States

Contribution
Investigation

Competing interests
No competing interests declared
James M Cheverud

Department of Biology, Loyola University, Chicago, United States

Contribution
Resources, Writing – review and editing

Competing interests
No competing interests declared
Clay F Semenkovich

Department of Medicine, Washington University School of Medicine, Saint Louis, United States

Contribution
Resources, Writing – review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-1163-1871
Heather A Lawson

Department of Genetics, Washington University School of Medicine, Saint Louis, United States

Contribution
Conceptualization, Funding acquisition, Methodology, Supervision, Writing - original draft, Writing – review and editing

For correspondence
lawson@wustl.edu

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-3550-5485

Funding

Diabetes Research Center

Clay F Semenkovich

Diabetes Research Center

Heather A Lawson

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

This work was supported by the Washington University Department of Genetics, the Diabetes Research Center at Washington University (grant P30DK020579), the NIH NIDDK (grant K01 DK095003), the NIEHS (grant U24ES026699), and the NIH NHGRI (grant T32-GM007067). The authors declare no conflicts of interest.

Ethics

Mouse colony was maintained at the Washington University School of Medicine and all experiments were approved by the Institutional Animal Care and Use Committee in accordance with the National Institutes of Health guidelines for the care and use of laboratory animals. Protocol #20-0384.

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.