Protein Phosphatase 1 Regulatory Subunit 3B Genotype at rs4240624 Has a Major Effect on Gallbladder Bile Composition.

The protein phosphatase 1 regulatory subunit 3B (PPP1R3B) gene is a target of farnesoid X receptor (FXR), which is a major regulator of bile acid metabolism. Both PPP1R3B and FXR have been suggested to take part in glycogen metabolism, which may explain the association of PPP1R3B gene variants with altered hepatic computed tomography attenuation. We analyzed the effect of PPP1R3B rs4240624 variant on bile acid composition in individuals with obesity. The study cohort consisted of 242 individuals from the Kuopio Obesity Surgery Study (73 men, 169 women, age 47.6 ± 9.0 years, body mass index 43.2 ± 5.4 kg/m2) with PPP1R3B genotype and liver RNA sequencing (RNA-seq) data available. Fasting plasma and gallbladder bile samples were collected from 50 individuals. Bile acids in plasma did not differ based on the PPP1R3B rs4240624 genotype. However, the concentration of total bile acids (109 ± 55 vs. 35 ± 19 mM; P = 1.0 × 10-5) and all individual bile acids (also 7α-hydroxy-4-cholesten-3-one [C4]) measured from bile were significantly lower in those with the AG genotype compared to those with the AA genotype. In addition, total cholesterol (P = 0.011) and phospholipid (P = 0.001) levels were lower in individuals with the AG genotype, but cholesterol saturation index did not differ, indicating that the decrease in cholesterol and phospholipid levels was secondary to the change in bile acids. Liver RNA-seq data demonstrated that expression of PPP1R3B, tankyrase (TNKS), Homo sapiens chromosome 8 clone RP11-10A14.5 (AC022784.1 [LOC157273]), Homo sapiens chromosome 8 clone RP11-375N15.1 (AC021242.1), and Homo sapiens chromosome 8, clone RP11-10A14 (AC022784.6) associated with the PPP1R3B genotype. In addition, genes enriched in transmembrane transport and phospholipid binding pathways were associated with the genotype.
Conclusion: The rs4240624 variant in PPP1R3B has a major effect on the composition of gallbladder bile. Other transcripts in the same loci may be important mediators of the variant effect.

adenosine triphosphate binding cassette subfamily C

Homo sapiens chromosome 8 clone RP11‐375N15.1

Homo sapiens chromosome 8 clone RP11‐10A14.5

Homo sapiens chromosome 8 clone RP11‐10A14

aquaporin

body mass index

7α‐hydroxy‐4‐cholesten‐3‐one

cholesterol saturation index

computed tomography

expression quantitative trait loci

false discovery rate

farnesoid X receptor

glycodeoxycholic acid

gene ontology

hexokinase

kilobase

laparoscopic Roux‐en‐Y gastric bypass

messenger RNA

nonalcoholic fatty liver disease

nonalcoholic steatohepatitis

nuclear receptor subfamily 1 group H member 4

protein phosphatase 1 regulatory subunit 3B

RNA sequencing

solute carrier

tankyrase

The protein phosphatase 1 regulatory subunit 3B gene (PPP1R3B) is located on chromosome 8p and encodes the regulatory subunit historically known as GL. GL regulates the activity of protein phosphatase 1, which activates glycogen synthase and inactivates glycogen phosphorylase, which in turn catalyzes the rate‐limiting step in glycogenolysis.( ) PPP1R3B polymorphisms have been associated with type 2 diabetes and plasma glucose,( , ) serum lipids and lipoproteins,( , ) and plasma lactate levels.( ) The meta‐analysis of genome‐wide association identified PPP1R3B as a genome‐wide significant loci for hepatic steatosis measured by computed tomography (CT). The variant rs4240624, located 175 kilobase (kb) upstream of PPP1R3B and 229 kb downstream of tankyrase (TNKS), has been associated with decreased hepatic CT attenuation (the effect allele is a wild‐type A allele with allelic frequency of 90%( )). Similarly, the minor allele of rs4841132 variant (a proxy variant in complete linkage with rs4240624) has been associated with increased hepatic CT attenuation and with an increase in plasma levels of alanine aminotransferase but not with hepatic fat content,( , ) as also published for rs4240624.( ) Several factors can cause a change in hepatic CT attenuation, with increased hepatic fat, which decreases CT attenuation, being the most prevalent, followed by glycogen and iron content, both of which increase attenuation.( , ) Interestingly, the preprint of a recent study identifies rs4240624 as the strongest genome‐wide genetic determinant of liver volume whereas no association with liver fat or liver iron content was observed.( ) Overall, the mechanisms for how the PPP1R3B polymorphism contributes to liver CT attenuation or liver metabolism are unclear.

It has been suggested that PPP1R3B is a target of farnesoid X receptor (FXR, also known as nuclear receptor subfamily 1 group H member 4 [NR1H4]),( ) which not only is an important regulator of bile acid metabolism but also affects glycogen metabolism.( ) On the other hand, bile acid metabolism is linked to both nonalcoholic fatty liver disease (NAFLD)( ) and obesity.( ) Therefore, we analyzed the effect of the PPP1R3B rs4240624 variant on bile acid composition in individuals with obesity.

Participants and Methods

Subjects

The study cohort consisted of 242 individuals from the Kuopio Obesity Surgery Study( ) with rs4240624 (PPP1R3B) genotype and liver RNA sequencing (RNA‐seq) data available (Table 1). Of these, 50 individuals had bile acids measured from both plasma and gallbladder bile (Supporting Table S1).

TABLE 1

Characteristics of the Liver mRNA Cohort (n = 242) Based on the PPP1R3B rs4240624 Genotype

	AA n = 170	AG n = 67	GG n = 5	P Value (K‐W)
Sex (male/female)	48/122	22/45	3/2	0.267*
Age (years)	47.8 ± 9.1	46.2 ± 8.7	52.7 ± 8.1	0.205
BMI (kg/m2)	43.5 ± 5.7	42.5 ± 4.7	44.1 ± 3.2	0.430
Total cholesterol (mmol/L)	4.2 ± 0.9	4.1 ± 0.9	3.9 ± 0.5	0.067
LDL cholesterol (mmol/L)	2.3 ± 0.8	2.3 ± 0.8	1.8 ± 0.4	0.327
HDL cholesterol (mmol/L)	1.1 ± 0.3	1.1 ± 0.3	1.0 ± 0.2	0.697
Triglycerides (mmol/L)	1.6 ± 0.3	1.1 ± 0.3	1.0 ± 0.4	0.093
Type 2 diabetes (%)	36	43	0	0.133*
Fasting glucose (mmol/L)	6.5 ± 1.8	6.7 ± 2.0	5.7 ± 0.5	0.524
Fasting insulin (mU/L)	18.5 ± 11.8	21.4 ± 22.1	14.8 ± 8.0	0.439
ALT (U/L)	46 ± 31	39 ± 20	35 ± 16	0.312
Histology (n = 239)
Steatosis grade (n)				0.403*
<5%	59	32	3
5%‐33%	65	19	1
33%‐66%	24	6	1
>66%	21	8	0
Lobular inflammation grade (n)				0.503*
None	109	49	10
<2 foci per 200× field	49	12	4
2‐4 foci per 200× field	11	4	0
>4 foci per 200× field	0	0	0
Fibrosis stage (n)				0.906*
0	90	40	2
1	68	21	3
2	7	2	0
3	3	2	0
4	1	0	0

chi‐square test.

Abbreviations: ALT, alanine aminotransferase; HDL, high‐density lipoprotein; K‐W, Kruskal‐Wallis; LDL, low‐density lipoprotein.

The study protocol was approved by the Ethical Committee of Northern Savo and followed the Declaration of Helsinki. Written consent was obtained from all study subjects.

Clinical and Laboratory Measurements

Plasma glucose, insulin, serum lipids, and lipoprotein lipids were determined, as described.( )

Bile Acid Measurements

Bile acids were measured from fasting serum samples and from samples obtained from transhepatic gallbladder punctuation during elective laparoscopic Roux‐en‐Y gastric bypass (LRYGB). Measurements were carried out using a Nexera X2 ultra high‐performance liquid chromatography system (Shimadzu, Kyoto, Japan) coupled to a 5500 Qtrap mass spectrometer interfaced with an electrospray ion source (ABSciex, Toronto, Canada), as described.( ) Total, primary, secondary, tertiary, primary conjugated, and secondary conjugated bile acids were calculated based on individual bile acids. Secondary conjugated bile acids in plasma were calculated without taurolithocholate acid and secondary conjugated bile acids in bile were calculated without glycolithocholic acid because these were not measured.

Cholesterol and phospholipids from bile samples were measured enzymatically using a Clariostar analyzer (BMG Labtech, Ortenberg, Germany), as described.( ) Total lipid content and cholesterol saturation index (CSI) were calculated based on the critical tables by Carey.( )

Liver Biopsies and Histology

Liver biopsies were obtained with ultrasonic scissors during the elective LRYGB operation. Overall histologic assessment of liver biopsy samples was performed by one pathologist according to standard criteria.( , )

Genotyping

PPP1R3B rs4240624 was genotyped with TaqMan single‐nucleotide polymorphism genotyping assays (Applied Biosystems) according to the manufacturer’s protocol.

Liver RNA‐Seq Data

RNA samples were isolated using the miRNeasy (Qiagen) kit, and sequencing libraries were prepared using Ribo‐Zero Gold (Illumina). Reads were aligned against the human reference genome GRCh38 (release 29) using two‐pass STAR RNA‐aligner.( ) On average, 37.5 million (77.4%) read‐pairs mapped to a unique site of the human genome and were used for further analysis. The Rsubread R package( ) was used to count all reads mapping within exon features. To reduce random variability resulting from low counts, we required at least 80% of samples to have at least 10 reads, resulting in 15,243 genes in the final analysis. The gene level count values were then normalized using a trimmed‐mean of M values converted to count per million using edgeR( ) and inverse normal transformed. Technical factors were obtained from STAR RNA‐aligner (summary of mapping statistics) and from Picard( ) (quality metrics). Principal component analysis (PCA) was used to identify effects of confounding factors on the transcriptome. Expression data were corrected for uniquely aligned reads percentage and 3′ bias, body mass index (BMI), sex, and age for the subsequent statistical analyses.

Statistical Analysis

Data are presented as mean ± SD. Differences in categorical variables were evaluated by the χ 2 test. Differences between the study groups were compared using Kruskal‐Wallis and Mann‐Whitney independent samples test. Spearman’s rank correlation was used for correlation analyses. Analyses were conducted with SPSS version 25 (IBM Inc., Armonk, NY). P < 0.05 was considered statistically significant.

Pairwise gene correlations were measured with Pearson correlation coefficients. The Benjamini‐Hochberg (false discovery rate [FDR]) method was applied to adjust for multiple testing. An FDR value <0.05 was considered statistically significant.

Expression Quantitative Trait Loci Analysis

Expression quantitative trait loci (eQTL) analysis was used to investigate the effect of the rs4240624 genotype on gene expression levels in the liver. To improve power to map eQTL and remove the impact of cofactors, which can mask the genetic signal, we determined the hidden covariates with PCA. Matrix eQTL( ) was run in the R environment (www.r‐project.org) using the additive linear model with 23 principal components that were selected based on optimization on chromosome 21 as covariates to account for confounding variation. The rs4240624 variant was defined as cis‐eQTL if it were within 500 kb of either end of the gene and as trans‐eQTL if it were beyond that point. eGenes were defined as genes associated with eQTL (cis‐eGenes for local and trans‐eGenes for distant transcripts).

To reveal the functional relationship among key genes, gene ontology (GO) pathway analysis was performed on trans‐eGenes with nominal P < 0.01 as a cutoff, using the web‐based gene set analysis toolkit.( )

Results

PPP1RB3 rs4240624 Genotype Did Not Correlate With Liver Histology

From the initial 242 subjects (73 men, 169 women, age 47.6 ± 9.0 years, BMI 43.2 ± 5.4 kg/m2), 170 subjects had PPP1R3B rs4240624 genotype AA, 67 had genotype AG, and 5 had genotype GG. Clinical characteristics did not differ based on the genotype (Table 1). There were no significant associations between PPP1R3B rs4240624 genotype and liver histology (available from 239 subjects) in this cohort.

Bile Acids in Fasting Plasma Did Not Differ Based on the rs4240624 Genotype

From the whole study cohort, bile acids were measured from 50 subjects (16 men, 34 women, age 46.8 ± 10.0 years, BMI 41.7 ± 4.5 kg/m2); from those, 33 had PPP1RB3 rs4240624 genotype AA and 17 had genotype AG. Clinical characteristics did not differ based on the genotype (Supporting Table S1).

Total bile acids were not associated with the rs4240624 genotype (Fig. 1). However, glycodeoxycholic acid (GDCA) was higher in those with the AG genotype compared to those with the AA genotype (0.14 ± 0.11 μM vs. 0.31 ± 0.26 μM; P = 0.034). Importantly, serum 7a‐hydroxy‐4‐cholesten‐3‐one (C4), a marker of bile acid synthesis,( ) did not associate with the genotype (0.07 ± 0.08 μM vs. 0.05 ± 0.06 μM; P = 0.241) (Supporting Table S2).

FIG. 1

Total bile acid levels based on the PPP1R3Bgenotypes AA and AG at rs4841132. (A) In plasma. (B) In gallbladder bile. Kruskal‐Wallis test used for statistics. Horizontal line indicates the median, whiskers indicate the minimum and maximum values and box borders indicate lower and upper quartiles. Abbreviation: BA, bile acid.

Total Bile Acid Levels in Bile Were Significantly Lower in Those With the PPP1R3B rs4240624 AG Genotype

Next, we analyzed if the PPP1R3B rs4240624 genotype had an effect on bile acid concentration in samples taken from gallbladder. The total concentration of bile acids in bile was much lower in those with the AG genotype compared to the AA genotype (108.7 ± 55.4 vs. 35.1 ± 19.0 mM; P = 1.0 × 10–5) (Fig. 1). Similarly, all individual bile acid concentrations in bile, including C4 (0.001 ± 0.001 vs. 0.0002 ± 0.001 mM; P = 0.015), were lower in those with the AG genotype (Table 2).

TABLE 2

Different Bile Acids in Bile Based on the PPP1R3B rs4240624 Genotype

Bile Acids in Bile (mM)	AA (n = 33)	AG (n = 17)	P Value (K‐W)
Primary BAs	0.030 ± 0.035	0.009 ± 0.014	3.3 × 10−4
Cholic acid	0.017 ± 0.024	0.006 ± 0.011	0.002
Chenodeoxycholic acid	0.012 ± 0.013	0.003 ± 0.004	3.1 × 10−4
Secondary BAs	0.005 ± 0.004	0.002 ± 0.001	1.6 × 10−4
Deoxycholic acid	0.004 ± 0.004	0.001 ± 0.001	3.8 × 10−4
Lithocholic acid	0.0003 ± 0.0002	0.0001 ± 0.0001	0.006
Tertiary BAs
Ursodeoxycholic acid	0.001 ± 0.001	0.0003 ± 0.0004	0.041
Primary conjugated BAs	79.3 ± 48.5	22.2 ± 12.2	5.7 × 10−6
Glycocholic acid	30.7 ± 20.0	9.6 ± 5.4	2.8 × 10−4
Taurocholic acid	10.5 ± 7.1	3.3 ± 2.6	3.2 × 10−5
Glycochenodeoxycholic acid	31.4 ± 25.8	7.5 ± 4.6	5.8 × 10−5
Taurochenoxycholic acid	6.7 ± 5.6	1.8 ± 1.6	1.7 × 10−5
Secondary conjugated BAs*	29.4 ± 15.2	12.9 ± 8.6	1.4 × 10−4
Glycodehydrocholic acid	18.6 ± 10.6	9.0 ± 6.5	0.004
Taurodeoxycholic acid	6.9 ± 4.5	2.8 ± 2.7	0.001
Taurolithocholate	0.2 ± 0.2	0.1 ± 0.1	0.037
Glycoursodeoxycholic acid	3.3 ± 2.9	0.9 ± 0.7	1.1 × 10−4
Tauroursodeoxycholic acid	0.4 ± 0.3	0.1 ± 0.1	3.8 × 10−5
Hyodeoxycholic acid	0.0002 ± 0.0003	0.00004 ± 0.0001	0.018
C4	0.001 ± 0.001	0.0002 ± 0.001	0.015

Missing glycolithocholic acid (GLCA).

Abbreviations: BA, bile acid; K‐W, Kruskal‐Wallis.

Next, we calculated ratios of different bile acid to total bile acid concentration in bile to estimate possible changes independent of the change in concentration of total bile acids. The ratio of primary conjugated bile acids to total bile acids was lower in those with the AG genotype. In contrast, the ratio of secondary conjugated bile acids to total bile acids (glycolithocholic acid was not available for our analysis) was higher in those with the AG genotype. In addition, the ratio of GDCA to total bile acids was higher in those with the AG genotype (Supporting Table S3).

Total Lipid Content in Bile Is Decreased in Those With the PPP1R3B AG Genotype

Because bile acids are known to drive the biliary secretion of phospholipid and cholesterol, we also measured the concentration of these lipids in the bile. Similar to the level of bile acids, cholesterol (25.7 ± 22.3 vs. 13.3 ± 11.1 mM; P = 0.011) and phospholipid (74.1 ± 36.3 vs. 36.9 ± 19.9 mM; P = 0.001) concentrations were lower in subjects with the AG genotype compared to those with the AA genotype. Concordantly, total lipid content in bile was lower in those with the AG genotype (16.6 ± 6.6 vs. 10.1 ± 5.3 g/dL; P = 0.003), but the CSI did not differ (99.8% ± 60.2% vs. 101.2% ± 46.8%; P = 0.444) (Fig. 2).

FIG. 2

Lipids in gallbladder bile based on thePPP1R3Bgenotypes AA and AG at rs4240624. (A) Cholesterol, (B) phospholipid, (C) total lipid content, and (D) CSI. Kruskal‐Wallis test used for statistics. Horizontal line indicates the median, whiskers indicate the minimum and maximum values and box borders indicate lower and upper quartiles.

Liver‐Specific Impact of rs4240624 on PPP1R3B, TNKS, AC022784.1, AC021242.1, and AC022784.6 Expression

Because changes in gene expression are one of the most important determinants of phenotypic changes, we investigated the effect of the rs4240624 genotype on messenger RNA (mRNA) expression by using the liver RNA‐seq data from the whole cohort (n = 242). The genomic region surrounding rs4240624 (±500 kb) is comprised of 25 genes (four protein coding genes, five pseudogenes, and 16 RNA genes), seven of which were detected in the liver (Fig. 3). Interestingly, while all four protein coding genes were expressed both in the liver and in subcutaneous adipose tissue collected from the same individuals (data not shown), the expression of Homo sapiens chromosome 8, clone RP11‐10A14 (AC022784.6), Homo sapiens chromosome 8 clone RP11‐375N15.1 (AC021242.1), and Homo sapiens chromosome 8 clone RP11‐10A14.5 map 8 (AC022784.1) were detected only in the liver. We identified five genes (cis‐eGenes) regulated by the rs4240624 variant: PPP1R3B, AC022784.6, AC021242.1, AC022784.1, and TNKS (β = 0.787, FDR, 7.1 × 10−30; β = 0.360, FDR, 3.4 × 10−4; β = −1.037, FDR, 1.0 × 10−32; β = −0.966, FDR, 4.8 × 10−35; and β = −0.769, FDR; 3.8 × 10−17, respectively) (Fig. 3). Importantly, rs4240624 did not affect subcutaneous adipose tissue expression levels of any of the detected genes in cis (Fig. 3; data not shown), indicating that association of rs4240624 with expression levels of local genes is tissue specific. Additionally, we detected 202 trans‐eGenes with a nominal P < 0.01.

FIG. 3

Genes detected in the liver. (A) Genomic region of rs424062 (yellow boxes) depict genes detected in the liver samples. Genomic coordinates are based on the Human Genome Assembly GRCh38 (hg38). (B) NAFLD‐associated variant rs4240624 is a cis‐eQTL for five genes in the liver (AC022784.1, AC021242.1, TNKS, PPP1R3B, AC022784.6) but not in the subcutaneous adipose tissue. Violin plots display gene expression levels in liver based on the rs4240625 genotype (AA, AG, and GG); effect sizes (β) and cis‐eQTL FDR values are presented. Overlaid boxplot represent upper and lower quartiles, with horizontal line depicting the median. Whiskers depict 1.5 times interquartile range. Abbreviations: ERl1, exoribonuclease 1; MFHAS1, malignant fibrous histiocytoma amplified sequence 1.

Expression of PPP1R3B, AC021242.1, and AC022784.6 Correlated With Total Bile Acids in Bile

Although correlations were not very strong, mRNA expression of PPP1R3B and to a lesser extent AC021242.1 and AC022784.6 tended to correlate with several bile acids in bile but not in plasma (except between PPP1R3B and total bile acids) (Supporting Tables S4 and S5).

Pathway Analysis Revealed Changes in the Transmembrane Transporter Activity Library Based on the PPP1R3B Genotype

To further explore the biological functions of rs4240624‐mediated regulation of gene expression in the liver, we performed a GO analysis. Identification of trans‐eGenes is challenging due to the fact that the effect sizes of trans‐eQTLs are typically small and underpowered due to a heavy multiple testing correction burden.( ) Hence, we performed the GO analysis on nominally significant trans‐eGenes to get more insight into potentially affected pathways rather than particular genes. The overrepresented groups identified according to GO molecular function were involved in active transmembrane transporter activity (GO:0022804 enrichment ratio, 3.5; FDR, 0.028) and phospholipid binding (GO:0005543 enrichment ratio, 3.2; FDR, 0.028). In addition, we observed a significant overrepresentation of nuclear receptors (Wikipathways WP2882 enrichment ratio, 7.82; FDR, 1.1 × 10−4). Strikingly among the trans‐eGenes, we found FXR (NR1H4; nominal P = 0.007) (Table 3).

TABLE 3

GO and WikiPathways Enrichment Analysis of a Set of PPP1R3B rs4240624‐Associated trans‐eGenes

Identified eight enriched genes belonging to the molecular function active transmembrane transporter activity (GO:0022804), enrichment ratio, 3.5; FDR, 0.028
ENSG00000118777	ABCG2	ATP binding cassette subfamily G member 2 (junior blood group)
ENSG00000123191	ATP7B	ATPase copper transporting beta
ENSG00000124564	SLC17A3	solute carrier family 17 member 3
ENSG00000127249	ATP13A4	ATPase 13A4
ENSG00000128524	ATP6V1F	ATPase H+ transporting V1 subunit F
ENSG00000139370	SLC15A4	Solute carrier family 15 member 4
ENSG00000171365	CLCN5	Chloride voltage‐gated channel 5
ENSG00000173214	MFSD4B	Major facilitator superfamily domain containing 4B
Identified 10 enriched genes belonging to the molecular function phospholipid binding (GO:0005543), enrichment ratio, 3.2; FDR, 0.028
ENSG00000021762	OSBPL5	Oxysterol binding protein‐like 5
ENSG00000100003	SEC14L2	SEC14‐like lipid binding 2
ENSG00000110243	APOA5	Apolipoprotein A5
ENSG00000118960	HS1BP3	Hematopoietic cell‐specific Lyn substrate 1 binding protein 3
ENSG00000135317	SNX14	Sorting nexin 14
ENSG00000141738	GRB7	Growth factor receptor bound protein 7
ENSG00000154217	PITPNC1	Phosphatidylinositol transfer protein cytoplasmic 1
ENSG00000159023	EPB41	Erythrocyte membrane protein band 4.1
ENSG00000197858	GPAA1	Glycosylphosphatidylinositol anchor attachment 1
ENSG00000247596	TWF2	Twinfilin actin binding protein 2
Identified seven enriched genes belonging to the nuclear receptors meta‐pathway (WP2882), enrichment ratio, 7.8; FDR, 1.1 × 10−4
ENSG00000110243	APOA5	Apolipoprotein A5
ENSG00000114315	HES1	Hes family bHLH transcription factor 1
ENSG00000012504	NR1H4	Nuclear receptor subfamily 1 group H member 4
ENSG00000186350	RXRA	Retinoid X receptor alpha
ENSG00000072310	SREBF1	Sterol regulatory element binding transcription factor 1
ENSG00000118503	TNFAIP3	TNF alpha‐induced protein 3
ENSG00000136810	TXN	Thioredoxin
Identified four enriched genes belonging to the interleukin‐18 signaling pathway (WP4754), enrichment ratio, 5.4; FDR, 0.022
ENSG00000125817	CENPB	Centromere protein B
ENSG00000159023	EPB41	Erythrocyte membrane protein band 4
ENSG00000106397	PLOD3	Procollagen‐lysine,2‐oxoglutarate 5‐dioxygenase 3
ENSG00000118503	TNFAIP3	TNF alpha‐induced protein 3
Identified five enriched genes belonging to the VEGFA‐VEGFR2 signaling pathway (WP3888), enrichment ratio, 3.46; FDR, 0.031
ENSG00000175183	CSRP2	Cysteine and glycine rich protein 2
ENSG00000159023	EPB41	Erythrocyte membrane protein band 4
ENSG00000133961	NUMB	NUMB, endocytic adaptor protein
ENSG00000106397	PLOD3	Procollagen‐lysine,2‐oxoglutarate 5‐dioxygenase 3
ENSG00000136810	TXN	Thioredoxin

Abbreviations: ATP, adenosine triphosphate; ATPase, adenosine triphosphatase; TNF, tumor necrosis factor; VEGF, vascular endothelial growth factor; VEGFR, vascular endothelial growth factor receptor; WP, WikiPathways.

Transcripts in the PPP1R3B Loci Have Different Co‐Expression Patterns

Finally, to further explore how the changes in the expression of rs4240624‐associated cis‐eGenes might affect the expression of distant genes, we performed gene expression correlation analysis using Pearson’s correlation coefficient as a measure of gene co‐expression. Based on our results from bile samples and the role of PPP1R3B in glycogen metabolism, we focused on selected genes of bile acid synthesis, regulation, and transport and furthermore on genes taking part in glycogen metabolism and water transport (Fig. 4). We observed significant correlations between expression of PPP1R3B and several genes involved in glycogen metabolism, confirming its reported role in glycogenesis. Additionally, we observed correlations between PPP1R3B and genes involved in bile acid biosynthesis (cytochrome P450 family 7 subfamily B member 1 [CYP7B1], alpha‐methylacyl‐coenzyme A [CoA] racemase [AMACR], solute carrier family 27 member 5 [SLC27A5], bile acid‐CoA:amino acid N‐acyltransferase [BAAT], D‐box binding PAR bZIP transcription factor [DBP], and hydroxy‐delta‐5‐steroid dehydrogenase, 3 beta‐ and steroid delta‐isomerase 7 [HSD3B7]) and regulation (NR1H4 [i.e., FXR], ; i.e., NR1H3], and nuclear receptor subfamily 1 group I member 3 [NR1I3]). Furthermore, PPP1R3B expression correlated with canalicular multispecific organic anion transporter 2 (encoded by the adenosine triphosphate binding cassette subfamily C member 3 [ABCC] gene, more often called MRP3) involved in bile acid transport.( ) Moreover, we observed a correlation between PPP1R3B and aquaporin (AQP)3, AQP7, and AQP9. Finally, we found correlations between remaining cis‐eGenes (AC022784.6, AC021242.1, AC022784.1, and TNKS) and bile acid metabolism and transport‐related genes (Fig. 4).

FIG. 4

Pearson's correlation between identifiedcis‐eQTL genes and selected genes taking part in bile acid and glycogen metabolism and water transport. Pearson's correlation coefficients are shown with continuous gradient colors, with red representing positive and blue representing negative correlation. Asterisk symbols correspond to FDR‐adjusted *P< 0.05, **P< 0.01, ***P< 0.001. Abbreviations: ABCB11, ATP binding cassette subfamily B member 11 aka BSEP, bile salt export pump; ABCC2, ATP binding cassette subfamily C member 2; ABCC3, ATP binding cassette subfamily C member 3; ABCC4, ATP binding cassette subfamily C member 4; ACOX2, acyl‐CoA oxidase 2; AGL, amylo‐alpha‐1, 6‐glucosidase, 4‐alpha‐glucanotransferase; AKR1C4, 3α‐steroid dehydrogenase; AKR1D1, delta‐5‐3‐oxysterod 5β‐reductase; AMACR, alpha‐methylacyl‐CoA racemase; AQP1, aquaporin 1; AQP11, aquaporin 11; AQP3, aquaporin 3; AQP4, aquaporin 4; AQP7, aquaporin 7; AQP9, aquaporin 9; CYP27A1, sterol 27‐hydroxylase; CYP7A1, cholesterol 7a‐hydroxylase; CYP7B1, oxysterol 7a‐hydroxylase; CYP8B1, sterol 12α‐hydroxylase; DBP, D‐box binding PAR bZIP transcription factor; EPHX1, epoxide hydrolase 1; EPM2A, glucan phosphatase, laforin; G6PC, glucose 6‐phosphatase catalytic subunit; G6PC3, glucose 6‐phosphatase catalytic 3; GAA, glucosidase alpha, acid; GBE1, 1,4‐alpha‐glucan branching enzyme 1; GCK, glucokinase; GYG1, glycogenin 1; GYS2, glygocen synthase 2; HK1, hexokinase 1; HK2, hexokinase 2; HK3, hexokinase 3; HSD3B7, 3β‐hydroxy‐delta5‐C27‐steroid oxioreductase; LAMP2, lysosome‐associated membrane protein 2; LDHA, lactate dehydrogenase A; NHLRC1, NHL repeat containing E3 ubiquitin protein ligase 1; NR1H3, nuclear receptor subfamily 1 group H member 3; NR1H4, nuclear receptor subfamily 1 group H member 4 aka FXR, farnesoid X‐activated receptor; NR1I2, nuclear receptor subfamily 1 group I member 2; NR1I3, nuclear receptor subfamily 1 group I member 3; NROB2, nuclear receptor subfamily 0 group B member 2; PFKM, phosphofructokinase, muscle; PGM1, phosphoglycomutase 1; PHKA2, glygocen phosphorylase kinase regulatory subunit alpha 2; PHKB, phosphorylase kinase regulatory subunit beta; PHKG2, glygocen phosphorylase kinase catalytic subunit gamma 2; PPP1R3C, protein phosphatase 1 regulatory subunit 3C; PPP1R3D, protein phosphatase 1 regulatory subunit 3D; PYGL, glygocen phosphorylase L; SCP2, sterol carrier protein 2; SLC10A1, solute carrier family 10 member 1; SLC27A5, solute carrier family 27 member 5; SLC2A2, solute carrier family 2 member 2; SLC37A4, glucose 6‐phosphatase translocase; SLC51A, solute carrier family 51 alpha subunit; SLCO1A2, solute carrier organic anion transporter family member 1A2; SLCO1B1, solute carrier organic anion transporter family member 1B1; SLCO1B3, solute carrier organic anion transporter family member 1B3; SLCO2B1, solute carrier organic anion transporter family member 2B1; UAP1, UDP‐N‐acetylglucosamine pyrophosphorylase 1.

Interestingly, transcripts in the PPP1R3B loci showed different co‐expression patterns. PPP1R3B expression correlated positively with AC022784.6 (r = 0.38; FDR, 2.1 × 10−7) and negatively with TNKS and AC021242.1 (r = −0.37; FDR, 3.2 × 10−7; and r = −0.32; FDR, 1.6 × 10−5, respectively). On the other hand, TNKS showed a very strong positive correlation with AC022784.1 and AC021242.1 (r = 0.56; FDR, 2.0 × 10−17; and r = 0.80; FDR, 4.4 × 10−52, respectively) (Supporting Fig. S1). In addition, hexokinase 1 (HK1), HK2, HK3, phosphofructokinase, muscle (PFKM), glycogenin 1 (GYG1), glycogen synthase 2 (GYS2), and lysosomal associated membrane protein 2 (LAMP2) did not correlate with PPP1R3B but showed association with TNKS. Finally, we observed that TNKS correlated with ABCB11 (also known as bile salt export pump [BSEP]), the transporter responsible for canalicular transport of the major bile acid species, and other transporters, such as ABCC2, solute carrier organic anion transporter family member 1A2 (SLCO1A2), and SLCO1B3, whereas PPP1R3B did not.

Discussion

The major finding of this study is that PPP1R3B rs4240624 associated with bile acid concentration in gallbladder bile. In addition, liver mRNA expression of PPP1R3B showed positive correlation with total bile acids in plasma and negative correlation with total bile acids in gallbladder bile (Supporting Tables S4 and S5). Interestingly, the genotype had only minimal effect on bile acid concentrations in plasma or bile acid synthesis monitored by C4. Furthermore, cholesterol, phospholipid, and total lipid concentrations in gallbladder bile were lower in subjects with the AG genotype. However, the CSI did not associate with the genotype, explained by the fact that the decrease in concentrations of the biliary lipids was proportional.

Our study is the first to report the association of PPP1R3B rs4240624 genotype with altered biliary bile acid secretion. While variants of PPP1R3B at rs4240624 and rs4841132 have been linked with NAFLD and nonalcoholic steatohepatitis (NASH) in some studies,( , ) the totality of evidence in the literature does not support an association between this variant of PPP1R3B and NAFLD and NASH.( , ) Accordingly, we did not observe association between NAFLD and rs4240624 in this cohort. The link between PPP1R3B and bile acid metabolism has been suggested based on the finding that PPP1R3B was a target of FXR (NR1A4) in slices of human liver and in mouse liver samples.( ) It is also known that bile acids can increase glycogenesis and decrease gluconeogenesis through activation of FXR.( ) Thus, we analyzed FXR and several other genes involved in bile acid metabolism to reveal more information about the function of the PPP1R3B gene. Our gene expression analyses based on liver mRNA‐seq data showed that those with the G allele had higher PPP1R3B mRNA expression in the liver, which is in agreement with our previous report linking the rs4841132 variant (in complete linkage disequilibrium with rs4240624) with a lower risk of steatosis and fibrosis in a cohort of 1,388 individuals( ) but contrary to some other findings.( ) Moreover, liver mRNA expression of FXR correlated with PPP1R3B, strengthening the finding of the association of these two genes.( ) Notably, we found that rs4240624 is a trans‐eQTL for FXR (nominal P = 0.007). A direct interaction of PPP1R3B with FXR should have influenced bile acid synthesis but was not observed in this study, indicating more complex and indirect mechanisms. FXR might be an important mediator between PPP1R3B and bile acids, but we conclude that it did not explain the effect of the PPP1R3B genotype on bile acids.

On the other hand, gene expression quantitative trait locus analysis of ~16,000 liver genes followed by GO analysis revealed significant enrichment in the GO category involved in transport (active transmembrane transporter activity) (Table 3). Interestingly, this pathway included two genes from the solute carrier family, which are membrane‐bound transporters, that have substrates that include, for example, glucose, fatty acids, lipids, and bile salts; for example, SLC17A3 belongs to the vesicular glutamate transporter family.( ) It has been suggested that bile acids could link to glutamate metabolism through FXR.( )

Importantly, we observed that four other transcripts located in close proximity to rs4240624 were significantly associated with the PPP1R3B rs4240624 genotype. Three of these (TNKS, AC022784.1, and AC021242.1) were down‐regulated in subjects with the G allele and AC022784.6 was moderately up‐regulated (Fig. 3). Based on these transcriptomics results, we suggest that another possible link between PPP1R3B and altered bile composition in gallbladder bile could be TNKS. We found that TNKS mRNA expression in the liver was associated with the PPP1R3B genotype and liver mRNA expression of PPP1R3B (Figs. 3 and 4). TNKS regulates WNT/β‐catenin signaling,( ) which has been shown to play a role in bile acid synthesis and transport through regulation of FXR activation.( ) Importantly, we observed a correlation between TNKS and bile salt export pump (ABCB11), a critical player involved in the secretion of bile salts into bile.( )

We also found that the PPP1R3B genotype associated with AC022784.1, AC022784.6, and AC021242.1 mRNA expression in the liver (Fig. 3). Interestingly, LOC157273, one of the transcript variants of the AC022784.1 gene (AC022784.1‐248), is a long noncoding RNA that has been reported to down‐regulate PPP1R3B expression.( ) While we could not link AC022784.1 with PPP1R3B expression or with bile acid levels, we cannot exclude it because our analysis did not focus on the splice variant composition. On the other hand, we observed a correlation between expression of PPP1R3B and another long noncoding RNA, AC022784.6 (Fig. 4; Supporting Fig. S1). Furthermore, AC021242.1 had positive correlation with total bile acids in bile and negative correlation with the CSI (Supporting Tables S4 and S5). Therefore, our results suggest that AC021242.1 might function as a liver‐specific regulator of bile acid metabolism together with PPP1R3B.

Our results also provide further insight into earlier findings suggesting alterations in bile acid metabolism have been associated with NAFLD.( ) In our study, the secondary bile acid GDCA was increased in plasma and decreased in bile in subjects with the AG genotype compared to those with the AA genotype. However, the ratio of GDCA to total bile acids in bile was increased in those with the AG genotype. This bile acid has been shown to be decreased in the liver in patients with NASH.( ) In contrast, high serum concentrations of GDCA have been associated with liver injury.( ) We also found that the ratio of primary conjugated bile acids to total bile acids in bile was lower in AG while that to secondary conjugated bile acids to total bile acids was higher in AG than in AA phenotypes. Lower bile acid concentration probably influences the composition of the microbiota and hence may increase bile acid dehydroxylation.( )

Variants of PPP1R3B have not been associated with gallstones.( ) Intriguingly, our preliminary data with limited number of subjects with gallstone disease data available showed that 2 out of 25 (8%) subjects with genotype AA had gallstones compared to 5 out of 15 (33%) with genotype AG (P = 0.021). This observation will require further studies. Importantly, our findings that bile acids from bile associated with the polymorphism remained essentially the same when excluding those with gallstones, indicating the observed changes in bile acid levels are primarily associated with genotype not gallstones.

A major strength of this study is that we had fasting bile acid samples taken from gallbladder punctuation during the LRYGB, which in conjunction with plasma samples provided new insights into bile acid homeostasis. To our knowledge, there are only few human studies in more recent literature with bile acids measured from bile. The presence of liver mRNA sequencing data allowed us to investigate the molecular pathways regulated by the PPP1R3B genotype. A limitation of this study is the relatively small sample size, which might limit the statistical power, especially with RNA‐seq data. It is probably the reason why we could not link the genotype with liver histology, although we have done it previously in a larger cohort (rs4841132 with complete linkage disequilibrium with rs4240624).( )

In conclusion, we found that the PPP1R3B genotype at rs4240624 associates strongly with total bile acid concentration in bile but not in plasma. Because the genotype did not associate with phospholipid and cholesterol content in bile, our results strongly suggest that factors involved in biliary bile acid transport are primarily affected. The role of altered expression of TNKS and three noncoding genes associated with the rs4240624 variant in the liver needs to be investigated further.

Fig S1

Click here for additional data file.

Table S1‐S5

Click here for additional data file.