Genetic Mapping of Multiple Traits Identifies Novel Genes for Adiposity, Lipids, and Insulin Secretory Capacity in Outbred Rats

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.

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:

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)]⁻¹.

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:

where Φ⁻¹ () is the inverse cumulative distribution function of the standard normal distribution, y_i 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:

θ = 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 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.

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⁻⁶), 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 R² 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) (r² > 0.95), leaving 125,611 SNPs for genetic mapping.

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, K_SNP, 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, K_HAP, 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).

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 K_HAP 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:

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 g_s 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 R² of 0.5.

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

where H_s 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 log_e10 ∗ 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:

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.

We used Matrix eQTL (v2.3) (28) to map eQTLs with the same genotypes and kinship matrix K_SNP 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).

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 R² = 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:

where y is the normalized physiological phenotype, t_z is the normalized expression for gene z, and g^top 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 n_t genes with cis-eQTL within a given pQTL, the n_t 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.

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.

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).

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.

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).

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 (R_p) and genotypic (R_g) correlations between traits tracked each other closely but on average, R_g is 19% larger than R_p. Fat pad weight traits, body length (BLTail and BLNoTail), and BW were significantly positively correlated (0.4 < R_g < 0.95). These measures were more strongly correlated with insulin traits (0.25 < R_g < 0.73) than with glucose traits (−0.03 < R_g < 0.28). Fat pad weights were also significantly positively correlated with other metabolic traits: R_g ∼0.45 with TRIG and R_g ∼0.35 with CHOL.

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.

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).

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).

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 < R_p < 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).

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 (R_p ∼−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 (T₂₂ = 2.348, P = 0.028; and T₂₂ = 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.

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 (R_p ∼0.15) (Fig. 5). These data suggest that the ACI and WKY haplotypes decreased Rfx6 expression, leading to decreased WPIC.

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 (R_p ∼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 (43–46) 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.

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.