Despite the successes of human genome-wide association studies, the causal genes underlying most metabolic traits remain unclear. We used outbred heterogeneous stock (HS) rats, coupled with expression data and mediation analysis, to identify quantitative trait loci (QTLs) and candidate gene mediators for adiposity, glucose tolerance, serum lipids, and other metabolic traits. Physiological traits were measured in 1,519 male HS rats, with liver and adipose transcriptomes measured in >410 rats. Genotypes were imputed from low-coverage whole-genome sequencing. Linear mixed models were used to detect physiological and expression QTLs (pQTLs and eQTLs, respectively), using both single nucleotide polymorphism (SNP)– and haplotype-based models for pQTL mapping. Genes with cis-eQTLs that overlapped pQTLs were assessed as causal candidates through mediation analysis. We identified 14 SNP-based pQTLs and 19 haplotype-based pQTLs, of which 10 were in common. Using mediation, we identified the following genes as candidate mediators of pQTLs: Grk5 for fat pad weight and serum triglyceride pQTLs on Chr1, Krtcap3 for fat pad weight and serum triglyceride pQTLs on Chr6, Ilrun for a fat pad weight pQTL on Chr20, and Rfx6 for a whole pancreatic insulin content pQTL on Chr20. Furthermore, we verified Grk5 and Ktrcap3 using gene knockdown/out models, thereby shedding light on novel regulators of obesity.

Obesity and overweight are major risk factors for poor metabolic health, including dyslipidemia, cardiovascular disease, and diabetes (1). Both genetic and environmental factors contribute to obesity and related metabolic traits, with up to 70% of the population variance attributable to genetics (2). Although human genome-wide association studies (GWAS) have associated many genomic regions with lipid levels (3), fasting glucose and insulin (4,5), diabetes (6), BMI (7), and waist-to-hip ratio (WHR) (8,9), causal genes in most of these loci remain unknown. This is due, in part, to the challenge of collecting gene expression and other -omic phenotypes from relevant human tissues.

Animal models of obesity have several advantages over human studies; environmental and dietary factors are controllable, and -omic data such as transcriptomics can be collected from relevant tissues under controlled conditions. In particular, the heterogeneous stock (HS) rats are an established model for the genetic study of obesity (10,11), and we have recently shown functional concordance in adiposity-associated adipose tissue transcripts between HS rats and humans (12). The HS are descended from eight fully sequenced inbred founder strains (13,14). Their genomes are highly heterozygous, each chromosome being a fine-grained mosaic of the founder haplotypes, suitable for high-resolution genetic mapping (15,16).

In this study, we use the HS rat model to dissect the genetic architecture of metabolic and adiposity-related phenotypes. This work significantly extends our previous study (10), in which we mapped physiological quantitative trait loci (pQTLs) for adiposity traits. First, we use many more animals, assess additional metabolic traits, and measure gene expression in two metabolically active tissues. Second, instead of genotyping the animals at a limited set of single nucleotide polymorphisms (SNPs) with an array, we genotype at a comprehensive genome-wide set of sites using low-coverage whole-genome sequencing (LC-WGS) followed by imputation. Third, we use mediation analysis to identify candidate causal genes that act in a tissue-specific manner, showing that many pQTLs are likely regulated by multiple genes. Finally, we validate two novel genes, Krtcap3 in liver and Grk5 in adipose, using in vivo knockout (KO) and cell gene silencing models, respectively.

Animals and Phenotypes

The outbred National Institutes of Health HS rat colony was initiated in 1984 from eight inbred founder strains—ACI/N, BN/SsN, BUF/N, F344/N, M520/N, MR/N, WN/N, and WKY/N—and maintained to minimize inbreeding (13). Rats used in this study were maintained at the Medical College of Wisconsin (NMcwi:HS; RGD_2314009), weaned at 21 days of age, and housed two per cage in microisolation cages in a conventional facility using autoclaved bedding (Sani-chips from P.J. Murphy). They were given ad libitum access to autoclaved Teklad 5010 diet (Harlan Laboratories) and provided reverse-osmosis water chlorinated to 2 to 3 ppm. We ran 1,519 male rats from 323 families through the experimental pipeline below using ∼12 animals/batch, with a new batch being saved weekly. We measured the following phenotypes as described below, in Keele et al. (10), and in Supplementary Table 1:

  • Intraperitoneal glucose tolerance test (IPGTT)

  • Body weight (BW)

  • Basal fasting glucose (Gluc0)

  • Fasting serum insulin (Ins0)

  • Glucose area under the curve during the IPGTT (GluAUC)

  • Insulin AUC during the IPGTT (InsAUC)

  • Quantitative insulin-sensitive check (QUICKI)

  • BW at sacrifice at 17 weeks old (sacBW)

  • Body length with tail at sacrifice at 17 weeks old (BLTail)

  • Body length without tail at sacrifice at 17 weeks old (BLNoTail)

  • Fasting serum total cholesterol (CHOL)

  • Fasting serum triglycerides (TRIG)

  • Heart weight at sacrifice (Heart)

  • Pancreas weight at sacrifice (Pancreas)

  • Retroperitoneal visceral fat pads at sacrifice (RetroFat)

  • Epididymal visceral fat pads at sacrifice (EpiFat)

  • BW gain as the difference between BW at sacrifice and BW at IPGTT (i.e., sacBW − BW) (BWGain)

  • Whole pancreas insulin content (μg/g) (WPIC)

  • Left kidney weight at sacrifice (LftKid)

  • Right kidney weight at sacrifice (RigKid)

  • Both kidneys weight at sacrifice (TotKid)

Tissues collected at sacrifice were also adjusted to BW by dividing by sacBW (Heart/BW, pancreas/BW, RetroFat/BW, EpiFat/BW, LftKid/GW, RigKid/BW, and TotKid/BW).

The IPGTT was conducted at 16 weeks of age. After an overnight fast, we measured BW, Gluc0, and Ins0 from each rat. We then injected rats with 1 mg/kg glucose, collecting blood glucose and serum for insulin at 15, 30, 60, and 90 min postinjection. We calculated GluAUC as follows:

graphic

where serum glucose was measured at t = 0, 15, 30, 60, 90 min. InsAUC was calculated using the same formula, where serum insulin replaced serum glucose in the equation.

We calculated QUICKI, a measure of insulin sensitivity (17), as [log(Ins0) + log(Gluc0)]−1.

At 17 weeks of age, a subset of 1,144 rats was euthanized after an overnight fast, after which we measured sacBW, BLTail, and BLNoTail. We used an enzymatic detection method to measure CHOL and TRIG from trunk blood. We dissected and weighed the heart, pancreas, RetroFat, and EpiFat and calculated BWGain. WPIC and IPGTT serum insulin were measured using an ELISA kit (Alpco Diagnostics, Salem, NH). Liver and adipose tissues were snap-frozen for expression analyses. In order to collect urinary protein levels for another study (18), we euthanized 375 rats at 24 weeks and weighed LftKid, RigKid, and TotKid. All protocols were approved by the institutional animal care and use committee at Medical College of Wisconsin.

Prior to genetic analysis, tissue weights were adjusted to account for BW by dividing by sacBW. Each phenotype was transformed to approximate normality, as summarized in Supplementary Table 1. CHOL and TRIG were log transformed, while all other traits except GluAUC (see explanation below) were rank-inverse normal transformed (RINT) as defined below:

graphic

where Φ−1 () is the inverse cumulative distribution function of the standard normal distribution, yi is the original phenotype, and n is the number of observations. We found that RINT for GluAUC introduced more extreme values than were originally present, which might generate false-positive associations. Therefore, we applied a more conservative RINT transformation (TRINT), in which data are mapped to a normal distribution with tails truncated at the θ percentile:

graphic

θ = 0.01 was chosen for GluAUC as it showed reasonable-looking values that appear normal without outliers. This setting limits the most extreme points to ∼2.33 SD from the mean.

RNA Sequencing

RNA from liver and adipose tissue was extracted from subsets of 430 and 415 rats, respectively, selected to maximize genetic diversity (e.g., no more than 1 or 2 rats per family), while encompassing the phenotypic spectrum of fat pad weights. Liver tissue RNA sequencing (RNA-seq) was run by the Genomics Core at the Medical College of Wisconsin, while adipose tissue RNA-seq was run by the Genomics Core at Wake Forest University School of Medicine.

For liver, poly-A libraries were prepared on the Illumina NeoPrep platform using the TruSeq Stranded mRNA Library Prep Kit for NeoPrep (catalog number NP-202–1001; Illumina). For adipose tissue, Ribo depletion libraries were prepared using the TruSeq Stranded Total RNA with Ribo-Zero Gold Preparation kit (Illumina).

Libraries were run on an Illumina HiSeq 2500 to obtain 37-bp paired-end reads for liver and 75-bp single-end reads for adipose tissue. We used STAR (v2.6.1a) (19) to align reads to the reference Rn6.0, PICARD (v2.5.0) to remove PCR duplicates, and featureCount in R package Rsubread (20) to compute gene-level expression counts, which were later normalized by sequencing depth, gene length, and RNA composition using the DESeq2 R package (v1.24.0). We excluded very lowly expressed genes with average reads per sample less than one. The normalized expression of 18,358 genes from adipose and 16,796 genes from liver were then RINT-transformed for expression QTL (eQTL) mapping and further analyses.

LC-WGS and Genotype Imputation

DNA from the HS rats was sequenced at BGI using 143-bp paired-end Illumina reads with a mean coverage depth of 0.24×, followed by imputation using STITCH (v1.3.7) (21). The STITCH algorithm uses a FAST-PHASE–type approach, optimized for low-coverage sequence data. Preexisting Illumina high-coverage (depth of 24−28×) sequence with 143-bp paired-end reads from the 8 founder strains (14) provided a haplotype reference panel. After quality control (imputation info score >0.4 and Hardy-Weinberg equilibrium P value >10−6), we retained 4,865,047 imputed SNPs. Imputed genotypes were then compared with a subset of 989 rats genotyped with a 10K Rat SNP array (10) in which 3,253 SNPs were in common. The concordance and the correlation R2 were 0.98 and 0.91, respectively, indicating the imputed genotypes were generally accurate. For genetic mapping, we used Plink (v1.9) (22) to remove SNPs in high linkage disequilibrium (LD) (r2 > 0.95), leaving 125,611 SNPs for genetic mapping.

Statistical Analyses

Founder Haplotype Dosages and Kinship Matrices

Each HS rat chromosome is a mosaic of the eight-founder haplotype reference panel, which STITCH (21) represents as founder dosages to account for uncertainty in the haplotype assignment, either because of identity by state between founders or ambiguities caused by heterozygosity or genotyping error. The STITCH haplotype dosages, which are initialized to be the known founder haplotypes, are dynamically updated during imputation and therefore need not be identical to the original founders. Therefore, we recomputed founder haplotype dosages from the imputed 126K tagging SNP genotypes plus the original founder haplotype reference panel using R/qtl2.

To account for the unequal relatedness between HS rats, we estimated a kinship matrix in two ways. The first was a SNP-based additive kinship matrix, KSNP, computed from the tagging SNPs, which was used to estimate heritability and genetic correlations for the physiological traits, and in genetic mapping of the expression traits. Heritability and genetic correlations were estimated in bivariate linear mix models using model 1 below. The second was a haplotype-based kinship matrix, KHAP, used for genetic mapping of physiological traits. This was estimated in R/qtl2 (v0.24) (23) using the aforementioned founder haplotype dosage information from either all chromosomes (used for SNP-based mapping) or all but one chromosome (used for haplotype-based mapping).

Genetic Mapping and Analyses of Physiological Traits

The normalized physiological phenotypes and 126K tagging SNPs were used for all genetic analyses. We conducted both SNP-based and haplotype-based genetic mapping to identify pQTLs.

We mapped SNP-based pQTLs using imputed SNP dosages from STITCH with the kinship matrix KHAP via the miqtl R package (v1.1.2) (10). The mixed model for SNP-based pQTL mapping at the SNP indexed by position s is:

graphic

where y is the vector of normalized phenotypes, X is a design matrix of fixed effect covariates, including coat color, euthanasia order, tissue harvest order, and technician (Supplementary Table 1), α is the vector of regression coefficients for these fixed effects, and gs is the vector of SNP dosages with regression coefficient β. The error vector e has variance-covariance matrix V=KHAPσg2 + Iσe2, where σg2 and σe2 are, respectively, the additive genetic and environmental variance components for the trait. Pearson correlations between physiological traits were estimated in R. We tested the association between phenotype and SNP using generalized least squares by premultiplying the above equation by the inverse of the matrix square root of V. We estimated genome-wide significance thresholds by parametric bootstrap, performing genome scans on 500 simulations from a fitted null model (i.e., no QTL) (24,25). Confidence intervals (CIs) on pQTL position were defined using LD drop between the peak and neighboring SNPs using a threshold R2 of 0.5.

Haplotype-based mixed-model genetic mapping was performed using R/qtl2 with the haplotype dosages, together with the kinship KHAP, with:

graphic

where Hs is the matrix of haplotype dosages at site s. The genome-wide significance threshold was estimated using permutation, and pQTL positional CIs were defined by extending the location of the peak SNP to both sides until all of the SNPs within that region had a logarithm of the odds (LOD) score higher than the peak LOD score minus 1.5 (i.e., 1.5 LOD drop). We converted LOD scores to logP units given that under the null hypothesis, 2 loge10 ∗ LOD is distributed as a χ72 random variable.

Heritability and genetic correlations were estimated using GCTA (26) in a bivariate linear mixed model similar to model 1 but with the SNP effects excluded. Heritability was estimated as below:

graphic

Estimation of Founder Haplotype Effects at pQTLs

We refined our estimates of founder haplotype effects at each detected QTL using Diploffect (v0.11) (27) to determine which founders carried nonzero haplotype effects on physiological phenotypes (pQTLs) and gene expression (eQTLs) and whether the QTL is potentially multiallelic.

Genetic Mapping of Liver and Adipose Expression Levels

We used Matrix eQTL (v2.3) (28) to map eQTLs with the same genotypes and kinship matrix KSNP restricted to those rats with expression data. We tested transformed normalized gene expression levels for association with SNP dosages in a linear mixed model, as with SNP-based pQTL mapping, except the choice of covariates (Supplementary Table 1). We focused on cis-eQTLs in order to identify candidate mediators of pQTLs. A cis-eQTL for a gene was recorded if at least two associated SNPs with logP >5 mapped within 1 Mb of that gene’s location.

For adipose genes, using principal component (PC) analysis for gene expression, we identified 28 adipose genes for which expression is highly variable between individuals. Many of these genes are annotated to be expressed in muscle rather than adipose cells, and also, one gene (MPZ) is a marker for nerve cells. To control for tissue heterogeneity, we included nerve and muscle proxy variables as covariates in the association model for adipose genes. The nerve proxy was the transformed expression of the MPZ gene, and the muscle proxy was the first eigenvector in PC analysis of the transformed expression of the other 27 genes, which showed high variance. The first 10 PCs from genotypes were also included as fixed effects, to control for environmental and genetic effects, while batch (i.e., the pools in which samples were sequenced together) was modeled as random effects. For liver genetic analysis, the covariates included batch (defined as the combination of sequence date and lane, which we modeled as a random effect), RNA integrity number, the first PC from expression data, and the first 10 PCs from genotypes (as fixed effects).

Mediation Analysis

For those genes that map as a cis-eQTL and fall within a pQTL, we conducted mediation analysis. Specifically, we assessed whether a gene with a cis-eQTL containing at least two associated SNPs that falls within the interval of a pQTL (as defined by R2 = 0.5) acts as an intermediate for the physiological trait using the methodology in Keele et al. (10) and Baron and Kenny (29). Briefly, consider gene z with a cis-eQTL mapping within the interval of a pQTL. One potential explanation is that the eQTL and pQTL are driven by the same causal variant and that the action of the pQTL on the phenotype is partially or fully mediated through the transcript level of gene z. The evidence for this is evaluated by comparisons among three models:

graphic
graphic
graphic

where y is the normalized physiological phenotype, tz is the normalized expression for gene z, and gtop is the SNP genotype dosage for the top SNP within the pQTL. Positive evidence for mediation is reported when model 3 is significantly more predictive of the phenotype than model 4. Because this test is applied to all nt genes with cis-eQTL within a given pQTL, the nt P values are subject to multiple test corrections using the Benjamini-Hochberg false discovery rate, with significant mediators defined as genes with q values <0.1. For significant mediators, full mediation is implied when model 3 fails to be significantly more predictive than model 5 (nominal P value <0.05); otherwise, mediation is considered partial.

Transduction of 3T3-L1 Cells to Knock Down Grk5 Expression

We grew 3T3-L1 cells, a mouse preadipocyte cell line, under standard culture conditions (30). The Grk5 gene was silenced by infecting 3T3-L1 preadipocytes with lentiviral particles to deliver mouse gene-specific shRNA expression vectors (sc-39043-V; Santa Cruz Biotechnology, Santa Cruz, CA) in the presence of polybrene (8 µg/mL; sc-134220; Santa Cruz Biotechnology) according to the manufacturer’s protocol. Control shRNA, lentiviral Particles-A (sc-108080; Santa Cruz Biotechnology), were used as a negative control. Cells successfully transduced and stably expressing shRNA were selected using 2 µg/mL puromycin (A1113803; Gibco, Thermo Fisher Scientific). The RNA of control and puromycin-selected Grk5 knockdown 3T3-L1 (n = 3/genotype) cells were isolated and reverse transcribed into cDNA for real-time PCR quantification of Grk5 (TaqMan primer number Mm00517039_m1) and Gapdh (an endogenous control; TaqMan primer number Mm99999915_g1) to determine the extent of Grk5 knockdown.

Adipocyte Differentiation and TRIG Measurement

The 3T3-L1 preadipocytes were differentiated into adipocytes for 7 days as described previously (31,32). Cellular TRIG was measured from nonsaponified samples after solubilization in Triton X-100 using a commercially available kit (Wako L-Type Triglyceride M test) and normalized to protein content (Pierce BCA Protein Assay Kit) (32).

Creation of a Krtcap3 KO Rat Model

Krtcap3 KO rats were created on the WKY background using Crispr-Cas9 as previously described (33). Wild-type (WT) and KO rats were bred from heterozygous breeder pairs. Rats were maintained on normal chow (Lab Diet, Prolab RMH3000, cat. no. 5P00) and weighed at 6 weeks of age.

Data and Resource Availability

Phenotype data are provided in Supplementary Table 2 and have been deposited in the Rat Genome Database (RGD; https://rgd.mcw.edu, access number RGD:153344611). Genotype data have been deposited in RGD (RGD:153344611). RNA-seq data have been deposited in the Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/) under accession numbers GSE200453 (liver) and GSE196572 (adipose).

Metabolic Traits Are Highly Correlated and Show Strong Heritability

SNP-based heritabilities, genetic correlations, and phenotype correlations of the normalized physiological traits are plotted in Fig. 1 and tabulated in Supplementary Table 3. Statistics for total tissue weights (prior to dividing by BW) are shown in Fig. 1, while those for tissue weight fractions are in Supplementary Fig. 1. Heritabilities ranged between 0.18 and 0.60 and are typical for rat studies of these phenotypes (10). Heritability was highest for BW and fat pad weights (RetroFat and EpiFat) and lowest for BWGain and BLNoTail (0.54 to 0.6 vs. 0.18 to 0.22). Phenotypic (Rp) and genotypic (Rg) correlations between traits tracked each other closely but on average, Rg is 19% larger than Rp. Fat pad weight traits, body length (BLTail and BLNoTail), and BW were significantly positively correlated (0.4 < Rg < 0.95). These measures were more strongly correlated with insulin traits (0.25 < Rg < 0.73) than with glucose traits (−0.03 < Rg < 0.28). Fat pad weights were also significantly positively correlated with other metabolic traits: Rg ∼0.45 with TRIG and Rg ∼0.35 with CHOL.

pQTLs

All pQTL exceeding the 10% genome-wide significance thresholds and the strain distribution patterns (SDP) at the top SNPs of each pQTL are shown in Table 1 and Fig. 2. We mapped 14 SNP-based pQTLs and 19 haplotype-based pQTLs, of which 10 overlap. pQTL mapping interval widths range between 0.59 and 8.12 Mb for SNP-based pQTLs and between 0.07 and 7.61 Mb for haplotype-based pQTLs. Several pQTLs detected by both SNP- and haplotype-based mapping and, based on haplotype plots (Supplementary Fig. 2), are likely pleiotropic: pQTLs on Chr1:180Mb associate with TRIG, EpiFat/BW, and RetroFat/BW; pQTLs on Chr6:27Mb associate with CHOL, TRIG, and RetroFat/BW; pQTLs on Chr8:50Mb associate with CHOL and TRIG; pQTLs on Chr10:101Mb associate with WPIC and EpiFat/BW; and pQTLs on Chr19:52Mb associate with BW, BLNoTail, and EpiFat. We also identified a pQTL on Chr16 that associates with both BW and RetroFat/BW using only haplotype-based mapping. Other pQTLs mapped to overlapping regions but are unlikely to be pleiotropic based on having differing haplotype effects (Supplementary Fig. 2): Chr10:101Mb for BWGain overlaps the Chr10:101Mb QTL for WPIC and EpiFat/BW; and Chr20:6Mb maps TotKid/BW and EpiFat. Fat pad weights that were not divided by sacBW mapped to generally the same regions as those divided by sacBW, as shown in Supplementary Table 4. Several glucose and insulin measures lacked genome-wide significant pQTLs (Supplementary Fig. 3) despite strong heritability. We identified a likely misassembly in the Rn6.0 reference for Heart weight on Chr5, where two pQTL separated by 50 Mb were in strong LD (shown as diamonds in Fig. 2). On reference Rn3.4, these loci are linked, and the pQTLs represent a single region.

Cis-Expression QTLs

We analyzed 16,796 and 18,358 transcripts with detectable expression levels in liver and adipose, respectively. Using SNP-based mixed models, we mapped 2,267 liver and 2,226 adipose cis-eQTLs with Benjamini-Hochberg false discovery rates of 0.00024 for liver and 0.00028 for adipose (Supplementary Tables 1, 5, and 6).

Mediation Analysis to Identify Candidate Genes

We performed mediation analysis to identify candidate causal genes with cis-eQTLs in either tissue that overlapped a SNP-based pQTL. The number of genes with overlapping cis-eQTLs in each pQTL ranged from 1 to 15 in liver and from 1 to 14 in adipose (average number of 5 genes), cumulatively representing 238 candidate genes with cis-eQTLs within pQTLs (Supplementary Table 7). We tested if the expression of each gene could fully or partially explain the effect of the pQTL. We identified 4 full mediators (3 in liver and 1 in adipose) and 41 partial mediators (17 in liver and 24 in adipose), listed in Table 2. Many of the partial mediators were correlated with the phenotypic trait (Supplementary Table 8).

Full Mediator: Grk5

There were 12 genes within the pQTL on Chr1:280–282 Mb for RetroFat/BW, EpiFat/BW, and TRIG, but only liver Cacul1 and adipose Grk5 cis-eQTLs overlapped this region (Supplementary Table 5). Mediation analysis identified only adipose Grk5 as a full mediator for RetroFat/BW and EpiFat/BW and a partial mediator for TRIG (Fig. 3A–C). The estimated haplotype effects at the pQTLs showed HS founders WKY and BUF increase both fat pad weight and TRIG and Grk5 expression, consistent with the positive correlation between Grk5 and the phenotypes (0.11 < Rp < 0.16) (Fig. 3D). Thus, BUF and WKY haplotypes likely increase the expression of adipose Grk5, which then increases fat pad weight and TRIG levels.

Compared with control shRNA-treated cells, day 0 undifferentiated Grk5 shRNA-treated 3T3-L1 preadipocytes had reduced Grk5 mRNA levels by 50% (P = 0.029). Consistent with the findings from mediation analysis, there was significantly decreased adipocyte TRIG content in day 7 differentiated Grk5 knockdown 3T3-L1 adipocytes relative to control adipocytes (P = 0.003) (Fig. 3E).

Full Mediator: Krtcap3

There were 123 genes within the pQTL on Chr6 for RetroFat/BW, CHOL, and TRIG, of which 15 liver genes and 14 adipose genes had cis-eQTLs (Supplementary Table 7). Of these, Krtcap3 was the only full mediator for RetroFat/BW and a partial mediator for TRIG (Fig. 4A and B). The common top SNP for these traits (chr6_29615083) was private to founder WKY, which is associated with decreased fat pad weight and increased Krtcap3 expression. This was consistent with the negative correlation between Krtcap3 liver expression and associated phenotypes (Rp ∼−0.25) (Fig. 4C), suggesting that the WKY haplotype increases the expression of liver Krtcap3, leading to reduced accumulation of fat pad and TRIG. We used CRISPR-Cas9 to KO Krtcap3 on the WKY background strain and found that both male and female KO rats exhibit increased BW at 6 weeks of age relative to wild-type (WT) rats (T22 = 2.348, P = 0.028; and T22 = 2.274, P = 0.033, respectively) (Fig. 4D). Subsequent work has shown that Krtcap3 KO rats exhibit increased BW over time, with female rats showing increased fat pad weights but improved insulin sensitivity (33), providing support for this gene’s role in adiposity. Several additional genes were identified as partial mediators, with different genes acting in adipose versus liver (Table 2), indicating multiple causal genes may underlie this locus.

Full Mediator: Rfx6

There were 25 genes within the 1.4 Mb pQTL on Chr20:31.43–32.83 Mb for WPIC, of which 4 liver-expressed genes and 3 adipose-expressed genes mapped as cis-eQTLs. Of these, liver Rfx6 was a full mediator for WPIC. The founder haplotypes ACI and WKY were associated with decreased WPIC and Rfx6 expression, consistent with the positive correlation between Rfx6 and WPIC (Rp ∼0.15) (Fig. 5). These data suggest that the ACI and WKY haplotypes decreased Rfx6 expression, leading to decreased WPIC.

Full Mediator: Ilrun

There were 27 genes within the 0.79 Mb pQTL on Chr20:7.28–8.08 Mb for EpiFat/BW, of which 7 liver and 10 adipose genes mapped as cis-eQTLs. Although expressed in both tissues, only liver Ilrun was a full mediator for EpiFat/BW. The founder haplotype WKY was associated with both decreased EpiFat/BW and Ilrun liver expression. Their positive correlation (Rp ∼0.17) (Fig. 6) suggests WKY reduces Ilrun expression, which decreases EpiFat/BW.

Although many loci have been identified in human GWAS for BMI and WHR (8,9), the causal genes underlying many of these loci remain unknown. In this study, we addressed this challenge in outbred HS rats by collecting transcriptome data in two metabolically active tissues from the same animals for which we had in-depth metabolic phenotypes, enabling us to assess causality. We mapped 23 pQTLs for 13 metabolic traits, including 6 pleiotropic loci. We used mediation analysis to identify candidate causal genes, showing that multiple tissue-specific mediators underlie several loci. We identified four full mediators and verified two novel genes (Krtcp3 and Grk5) experimentally. Our work demonstrates the effectiveness of imputation from LC-WGS for genotyping outbred populations and the power of mediation analysis for identifying candidate causal genes that regulate metabolic phenotypes.

In addition to identifying causal mediators, this study extends our previous work (10,11) by investigating a broader spectrum of metabolic traits, including fasting lipids, glucose tolerance, insulin sensitivity, and WPIC. We also used both SNP-based and haplotype-based mapping, the latter of which can reveal which founder alleles have contrasting phenotypic effects and whether they have a more complex structure, potentially suggesting the presence of multiple linked causal variants (34). In addition, we genotyped using LC-WGS followed by imputation, thereby increasing the number of markers from previous work (10). We recommend the LC-WGS/imputation strategy as an economic way to genotype populations; it also provides data for future structural variation analysis (35).

By increasing the numbers of SNPs and animals, we replicated known adiposity loci (10) at higher significance and mapped additional loci. We also replicated a subset of BW and fat pad weight pQTLs from a mega-analysis in HS rats (11). In addition to 13 loci for BW and fat pad traits, we identified 3 loci for CHOL, 4 for TRIG, 2 loci for WPIC, and 4 loci for tissue weights. Surprisingly, we mapped no genome-wide significant loci for the glucose and insulin traits, despite strong heritability, suggesting these traits are highly polygenic and larger sample sizes are needed. We identified full mediators underlying four loci, with partial mediators underlying several additional loci. Liver and adipose mediators often differ, demonstrating that causal genes are likely tissue-specific. We critically discuss the evidence for causality for each full mediator below.

Adipose Grk5 is a full mediator for EpiFat/BW and RetroFat/BW on rat Chr1 and a partial mediator for TRIG (11). Previous work has shown that Grk5 null mice show protection against diet-induced obesity, reduced expression of adipogenesis genes, glucose tolerance, and insulin resistance compared with WT littermates (36,37). Human GWAS have associated Grk5 with type 2 diabetes in Chinese cohorts (38,39) and lipid levels in Europeans (7). Our current work demonstrates experimentally that Grk5 knockdown in 3T3-L1 cells leads to decreased adipose TRIG accumulation in adipocytes, supporting a role for this gene in adiposity. We have since demonstrated that Grk5 is critical for adipogenesis in 3T3-L1 preadipocytes, potentially by controlling ERK activity (M. Seramur, B. McDonald, X. Wang, J. Chou, A. Craddock, L.C. Solberg Woods, C.C. Key, unpublished observations). Within the same pQTL, there is a nonsynonymous coding variant in the gene Prlhr private to the BUF and WKY founders and previously suggested as a candidate (10,11). Although Prlhr is associated with feeding behavior (40,41), Prlhr is mainly expressed in brain and adrenal tissue, so our expression data could not shed light on its role as a mediator, and it is possible that both adipose Grk5 and brain Prlhr are causal.

Liver Ktrcap3 is a full mediator for RetroFat/BW at the Chr6 pQTL, replicating previous work (10). Krtcap3 does not fall within human GWAS loci for BMI or WHR, but has been identified as a pleiotropic gene for obesity, dyslipidemia, and type 2 diabetes using a multivariate analysis (42). The current work shows that Krtcap3 KO rats have increased body weight at 6 weeks of age relative to WT rats. We have since shown that both male and female Krtcap3 KO rats have increased body weight over time relative to WT rats, with females showing increased adiposity and insulin sensitivity (33). Interestingly, Krtcap3 is unlikely the only gene acting in this locus. Adcy3 also falls in this region and was previously nominated as a functional candidate in the rat (10,11) via a coding variant as opposed to expression levels. Adcy3 is associated with monogenic and polygenic obesity in humans (4346) and has been shown to act via coding variation that impacts protein folding or function (47,48) and altered expression (43). Based on our functional validation of Krtcap3 (33) coupled with strong support for Adcy3 in human studies, it is likely both genes are causal for adiposity. In addition, our work identifies several partial mediators in both liver and adipose tissue, emphasizing the complexity of the Chr6 adiposity locus.

Liver Ilrun expression is a full mediator for EpiFat/BW on rat Chr20. Ilrun was previously named LOC294154 in rats and C6orf106 in humans. C6orf106 associates with anthropometric traits, including height (49), weight (50), BMI (51), and body fat ratio (52) in humans and might drive persistent thinness versus severe obesity (53). The mouse ortholog is implicated in lipid metabolism and hepatic lipoprotein production (54), making this a highly attractive candidate at this locus.

Liver Rfx6 expression is a full mediator for WPIC at a separate Chr20 locus. Rfx6 plays a role in controlling the formation of pancreatic islets and production of insulin in mice and humans (55). Inactivation of Rfx6 in adult β-cells impairs insulin secretion via altered glucose sensing (56). In addition, Rfx6 associates with monogenic diabetes of the young (57). Although Rfx6 is not expressed in livers of mice or humans, it is expressed in liver of rats (58,59). We suspect that the liver signal that we found in this study is not tissue-specific and may be representative of the role of this gene in the pancreas. Together, this evidence suggests that genetic variation in Rfx6 contributes to metabolic dysfunction and/or diabetes susceptibility.

Using multiple metabolic phenotypes and transcriptome data from relevant tissues measured in the same animals, we identify four candidate causal genes and experimentally validate two of them. This work demonstrates the power of mediation analysis for delineating genetic effects underlying multiple related phenotypes. Furthermore, it emphasizes tissue specificity of genetic effects and the potentially complex underpinnings of several metabolic loci due to multiple candidate gene drivers.

This article contains supplementary material online at https://doi.org/10.2337/figshare.21295908.

R.M. and L.C.S.W. are co-last authors.

Acknowledgments. The authors thank Apurva Chitre (Department of Psychiatry, University of California, San Diego) for assistance in creating Porcupine plot, Michael Scott (School of Biological Sciences, University of East Anglia) for helping with STITCH, and Leilei Cui (Genetics Institute, University College London) for helping with eQTL mapping.

Funding. This study was funded by National Institutes of Health grants R01 DK 106386 (to L.C.S.W.), R01 DK120667 (to L.C.S.W.), and R01 DK 088975 (to L.C.S.W.) and National Institute of General Medical Sciences grant R35 GM127000 (to W.V.).

Duality of Interest. No potential conflicts of interest relevant to this article were reported.

Author Contributions. W.V., R.M., and L.C.S.W. conceived experiments. K.H. and L.C.S.W. conducted animal experiments. K.H. and O.S. extracted DNA or RNA. M.T., A.C., and G.H. were responsible for RNA-seq. T.H.-L., W.L.C., and G.R.K. performed statistical analysis. T.H.-L., W.V., R.M., and L.C.S.W. were responsible for interpretation and presentation of analyzed results. S.K.D., B.M., N.K.S., and C.-C.C.K. conducted 3T3-L1 knockdown experiments. L.C.S.W. oversaw all experimental work. W.V. and R.M. oversaw statistical analysis. T.H.-L., W.V., R.M., and L.C.S.W. wrote the manuscript. All authors approved the final version of the manuscript. L.C.S.W. is the guarantor of this work and, as such, had full access to all of the data in the study and takes responsibility for the integrity of the data and the accuracy of the data analysis.

1.
GBD 2015 Obesity Collaborators
.
Health effects of overweight and obesity in 195 countries over 25 years
.
N Engl J Med
2017
;
377
:
13
27
2.
Maes
HH
,
Neale
MC
,
Eaves
LJ
.
Genetic and environmental factors in relative body weight and human adiposity
.
Behav Genet
1997
;
27
:
325
351
3.
Klarin
D
,
Damrauer
SM
,
Cho
K
, et al.;
Global Lipids Genetics Consortium
;
Myocardial Infarction Genetics (MIGen) Consortium
;
Geisinger-Regeneron DiscovEHR Collaboration
;
VA Million Veteran Program
.
Genetics of blood lipids among ∼300,000 multi-ethnic participants of the Million Veteran Program
.
Nat Genet
2018
;
50
:
1514
1523
4.
Chen
J
,
Spracklen
CN
,
Marenne
G
, et al.;
Lifelines Cohort Study
;
Meta-Analysis of Glucose and Insulin-related Traits Consortium (MAGIC)
.
The trans-ancestral genomic architecture of glycemic traits
.
Nat Genet
2021
;
53
:
840
860
5.
Horikoshi
M
,
Mägi
R
,
van de Bunt
M
, et al.;
ENGAGE Consortium
.
Discovery and fine-mapping of glycaemic and obesity-related trait loci using high-density imputation
.
PLoS Genet
2015
;
11
:
e1005230
6.
Vujkovic
M
,
Keaton
JM
,
Lynch
JA
, et al.;
HPAP Consortium
;
Regeneron Genetics Center
;
VA Million Veteran Program
.
Discovery of 318 new risk loci for type 2 diabetes and related vascular outcomes among 1.4 million participants in a multi-ancestry meta-analysis
.
Nat Genet
2020
;
52
:
680
691
7.
Yengo
L
,
Sidorenko
J
,
Kemper
KE
, et al.;
GIANT Consortium
.
Meta-analysis of genome-wide association studies for height and body mass index in ∼700000 individuals of European ancestry
.
Hum Mol Genet
2018
;
27
:
3641
3649
8.
Liu
C-T
,
Monda
KL
,
Taylor
KC
, et al
.
Genome-wide association of body fat distribution in African ancestry populations suggests new loci
.
PLoS Genet
2013
;
9
:
e1003681
9.
Pulit
SL
,
Stoneman
C
,
Morris
AP
, et al.;
GIANT Consortium
.
Meta-analysis of genome-wide association studies for body fat distribution in 694 649 individuals of European ancestry
.
Hum Mol Genet
2019
;
28
:
166
174
10.
Keele
GR
,
Prokop
JW
,
He
H
, et al
.
Genetic fine-mapping and identification of candidate genes and variants for adiposity traits in outbred rats
.
Obesity (Silver Spring)
2018
;
26
:
213
222
11.
Chitre
AS
,
Polesskaya
O
,
Holl
K
, et al
.
Genome wide association study in 3,173 outbred rats identifies multiple loci for body weight, adiposity, and fasting glucose
.
Obesity (Silver Spring)
2020
;
28
:
1964
1973
12.
Crouse
WL
,
Das
SK
,
Le
T
, et al
.
Transcriptome-wide analyses of adipose tissue in outbred rats reveal genetic regulatory mechanisms relevant for human obesity
.
Physiol Genomics
2022
;
54
:
206
219
13.
Hansen
C
,
Spuhler
K
.
Development of the National Institutes of Health genetically heterogeneous rat stock
.
Alcohol Clin Exp Res
1984
;
8
:
477
479
14.
Ramdas
S
,
Ozel
AB
,
Treutelaar
MK
, et al
.
Extended regions of suspected mis-assembly in the rat reference genome
.
Sci Data
2019
;
6
:
39
15.
Woods
LCS
,
Mott
R
.
Heterogeneous stock populations for analysis of complex traits
.
Methods Mol Biol
2017
;
1488
:
31
44
16.
Solberg Woods
LC
,
Palmer
AA
.
Using heterogeneous stocks for fine-mapping genetically complex traits
.
Methods Mol Biol
2019
;
2018
:
233
247
17.
Solberg Woods
LC
,
Holl
KL
,
Oreper
D
,
Xie
Y
,
Tsaih
SW
,
Valdar
W
.
Fine-mapping diabetes-related traits, including insulin resistance, in heterogeneous stock rats
.
Physiol Genomics
2012
;
44
:
1013
1026
18.
Keele
GR
,
Prokop
JW
,
He
H
, et al
.
Sept8/SEPTIN8 involvement in cellular structure and kidney damage is identified by genetic mapping and a novel human tubule hypoxic model
.
Sci Rep
2021
;
11
:
2071
19.
Dobin
A
,
Davis
CA
,
Schlesinger
F
, et al
.
STAR: ultrafast universal RNA-seq aligner
.
Bioinformatics
2013
;
29
:
15
21
20.
Love
MI
,
Huber
W
,
Anders
S
.
Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2
.
Genome Biol
2014
;
15
:
550
21.
Davies
RW
,
Flint
J
,
Myers
S
,
Mott
R
.
Rapid genotype imputation from sequence without reference panels
.
Nat Genet
2016
;
48
:
965
969
22.
Purcell
S
,
Neale
B
,
Todd-Brown
K
, et al
.
PLINK: a tool set for whole-genome association and population-based linkage analyses
.
Am J Hum Genet
2007
;
81
:
559
575
23.
Broman
KW
,
Gatti
DM
,
Simecek
P
, et al
.
R/qtl2: software for mapping quantitative trait loci with high-dimensional data and multiparent populations
.
Genetics
2019
;
211
:
495
502
24.
Valdar
W
,
Holmes
CC
,
Mott
R
,
Flint
J
.
Mapping in structured populations by resample model averaging
.
Genetics
2009
;
182
:
1263
1277
25.
Solberg Woods
LC
,
Holl
K
,
Tschannen
M
,
Valdar
W
.
Fine-mapping a locus for glucose tolerance using heterogeneous stock rats
.
Physiol Genomics
2010
;
41
:
102
108
26.
Yang
J
,
Lee
SH
,
Goddard
ME
,
Visscher
PM
.
GCTA: a tool for genome-wide complex trait analysis
.
Am J Hum Genet
2011
;
88
:
76
82
27.
Zhang
Z
,
Wang
W
,
Valdar
W
.
Bayesian modeling of haplotype effects in multiparent populations
.
Genetics
2014
;
198
:
139
156
28.
Shabalin
AA
.
Matrix eQTL: ultra fast eQTL analysis via large matrix operations
.
Bioinformatics
2012
;
28
:
1353
1358
29.
Baron
RM
,
Kenny
DA
.
The moderator-mediator variable distinction in social psychological research: conceptual, strategic, and statistical considerations
.
J Pers Soc Psychol
1986
;
51
:
1173
1182
30.
Zebisch
K
,
Voigt
V
,
Wabitsch
M
,
Brandsch
M
.
Protocol for effective differentiation of 3T3-L1 cells to adipocytes
.
Anal Biochem
2012
;
425
:
88
90
31.
Cuffe
H
,
Liu
M
,
Key
C-CC
, et al
.
Targeted deletion of adipocyte Abca1 (ATP-binding cassette transporter A1) impairs diet-induced obesity
.
Arterioscler Thromb Vasc Biol
2018
;
38
:
733
743
32.
Key
C-CC
,
Liu
M
,
Kurtz
CL
, et al
.
Hepatocyte ABCA1 deletion impairs liver insulin signaling and lipogenesis
.
Cell Rep
2017
;
19
:
2116
2129
33.
Szalanczy
AM
,
Goff
E
,
Seshie
O
, et al
.
Keratinocyte-associated protein 3 is a novel gene for adiposity with differential effects in males and females
.
Front Genet
2022
;
13
:
942574
34.
Baud
A
,
Hermsen
R
,
Guryev
V
, et al.;
Rat Genome Sequencing and Mapping Consortium
.
Combined sequence-based and genetic mapping analysis of complex traits in outbred rats
.
Nat Genet
2013
;
45
:
767
775
35.
Imprialou
M
,
Kahles
A
,
Steffen
JG
, et al
.
Genomic rearrangements in Arabidopsis considered as quantitative traits
.
Genetics
2017
;
205
:
1425
1441
36.
Wang
F
,
Wang
L
,
Shen
M
,
Ma
L
.
GRK5 deficiency decreases diet-induced obesity and adipogenesis
.
Biochem Biophys Res Commun
2012
;
421
:
312
317
37.
Wang
L
,
Shen
M
,
Wang
F
,
Ma
L
.
GRK5 ablation contributes to insulin resistance
.
Biochem Biophys Res Commun
2012
;
429
:
99
104
38.
Li
H
,
Gan
W
,
Lu
L
, et al.;
DIAGRAM Consortium
;
AGEN-T2D Consortium
.
A genome-wide association study identifies GRK5 and RASGRP1 as type 2 diabetes loci in Chinese Hans
.
Diabetes
2013
;
62
:
291
298
39.
Xia
Z
,
Yang
T
,
Wang
Z
,
Dong
J
,
Liang
C
.
GRK5 intronic (CA)n polymorphisms associated with type 2 diabetes in Chinese Hainan Island
.
PLoS One
2014
;
9
:
e90597
40.
Matta
CA
.
Maternal treatment with teratogen causes congenital malformations in mouse embryos
.
Folia Morphol (Praha)
1990
;
38
:
12
18
41.
Yang
J
,
Wu
X
,
Wu
X
, et al
.
The multiple roles of XBP1 in regulation of glucose and lipid metabolism
.
Curr Protein Pept Sci
2017
;
18
:
630
635
42.
Chen
Y-C
,
Xu
C
,
Zhang
J-G
, et al
.
Multivariate analysis of genomics data to identify potential pleiotropic genes for type 2 diabetes, obesity and dyslipidemia using Meta-CCA and gene-based approach
.
PLoS One
2018
;
13
:
e0201173
43.
Saeed
S
,
Bonnefond
A
,
Tamanini
F
, et al
.
Loss-of-function mutations in ADCY3 cause monogenic severe obesity
.
Nat Genet
2018
;
50
:
175
179
44.
Grarup
N
,
Moltke
I
,
Andersen
MK
, et al
.
Loss-of-function variants in ADCY3 increase risk of obesity and type 2 diabetes
.
Nat Genet
2018
;
50
:
172
174
45.
Nordman
S
,
Abulaiti
A
,
Hilding
A
, et al
.
Genetic variation of the adenylyl cyclase 3 (AC3) locus and its influence on type 2 diabetes and obesity susceptibility in Swedish men
.
Int J Obes
2008
;
32
:
407
412
46.
Speliotes
EK
,
Willer
CJ
,
Berndt
SI
, et al.;
MAGIC
;
Procardis Consortium
.
Association analyses of 249,796 individuals reveal 18 new loci associated with body mass index
.
Nat Genet
2010
;
42
:
937
948
47.
Stergiakouli
E
,
Gaillard
R
,
Tavaré
JM
, et al
.
Genome-wide association study of height-adjusted BMI in childhood identifies functional variant in ADCY3
.
Obesity (Silver Spring)
2014
;
22
:
2252
2259
48.
Toumba
M
,
Fanis
P
,
Vlachakis
D
, et al
.
Molecular modelling of novel ADCY3 variant predicts a molecular target for tackling obesity
.
Int J Mol Med
2022
;
49
:
10
49.
Weedon
MN
,
Lango
H
,
Lindgren
CM
, et al.;
Diabetes Genetics Initiative
;
Wellcome Trust Case Control Consortium
;
Cambridge GEM Consortium
.
Genome-wide association analysis identifies 20 loci that influence adult height
.
Nat Genet
2008
;
40
:
575
583
50.
Cho
H-W
,
Jin
H-S
,
Eom
Y-B
.
A genome-wide association study of novel genetic variants associated with anthropometric traits in Koreans
.
Front Genet
2021
;
12
:
669215
51.
Ahmad
S
,
Poveda
A
,
Shungin
D
, et al
.
Established BMI-associated genetic variants and their prospective associations with BMI and other cardiometabolic traits: the GLACIER Study
.
Int J Obes
2016
;
40
:
1346
1352
52.
Rask-Andersen
M
,
Karlsson
T
,
Ek
WE
,
Johansson
Å
.
Genome-wide association study of body fat distribution identifies adiposity loci and sex-specific genetic effects
.
Nat Commun
2019
;
10
:
339
53.
Riveros-McKay
F
,
Mistry
V
,
Bounds
R
, et al.;
Understanding Society Scientific Group
.
Genetic architecture of human thinness compared to severe obesity
.
PLoS Genet
2019
;
15
:
e1007603
54.
Bi
X
,
Kuwano
T
,
Lee
PC
, et al
.
ILRUN, a human plasma lipid GWAS locus, regulates lipoprotein metabolism in mice
.
Circ Res
2020
;
127
:
1347
1361
55.
Smith
SB
,
Qu
H-Q
,
Taleb
N
, et al
.
Rfx6 directs islet formation and insulin production in mice and humans
.
Nature
2010
;
463
:
775
780
56.
Piccand
J
,
Strasser
P
,
Hodson
DJ
, et al
.
Rfx6 maintains the functional identity of adult pancreatic β cells
.
Cell Rep
2014
;
9
:
2219
2232
57.
Patel
KA
,
Kettunen
J
,
Laakso
M
, et al
.
Heterozygous RFX6 protein truncating variants are associated with MODY with reduced penetrance
.
Nat Commun
2017
;
8
:
888
58.
National Library of Medicine, National Center for Biotechnology Information
.
Rfx6 regulatory factor X, 6 [Rattus norvegicus (Norway rat)], 23 June 2022
.
59.
Yu
Y
,
Fuscoe
JC
,
Zhao
C
, et al
.
A rat RNA-Seq transcriptomic BodyMap across 11 organs and 4 developmental stages
.
Nat Commun
2014
;
5
:
3230
Readers may use this article as long as the work is properly cited, the use is educational and not for profit, and the work is not altered. More information is available at https://www.diabetesjournals.org/journals/pages/license.