We postulated that type 2 diabetes (T2D) predisposes patients to exocrine pancreatic diseases through (epi)genetic mechanisms. We explored the methylome (using MethylationEPIC arrays) of the exocrine pancreas in 141 donors, assessing the impact of T2D. An epigenome-wide association study of T2D identified hypermethylation in an enhancer of the pancreatic lipase–related protein 1 (PNLIPRP1) gene, associated with decreased PNLIPRP1 expression. PNLIPRP1 null variants (found in 191,000 participants in the UK Biobank) were associated with elevated glycemia and LDL cholesterol. Mendelian randomization using 2.5M SNP Omni arrays in 111 donors revealed that T2D was causal of PNLIPRP1 hypermethylation, which in turn was causal of LDL cholesterol. Additional AR42J rat exocrine cell analyses demonstrated that Pnliprp1 knockdown induced acinar-to-ductal metaplasia, a known prepancreatic cancer state, and increased cholesterol levels, reversible with statin. This (epi)genetic study suggests a role for PNLIPRP1 in human metabolism and exocrine pancreatic function, with potential implications for pancreatic diseases.
We performed this study to identify epigenetic changes with type 2 diabetes (T2D) in the pancreas.
This study addresses whether T2D induces epigenetic changes that could explain why individuals with T2D are more prone to pancreatic diseases.
We found hypermethylation at PNLIPRP1 associated with T2D and identified a role of this gene in cholesterol metabolism.
This study has important implications in the prevention of pancreatic diseases, because their molecular mechanisms remain largely unknown.
Introduction
The pancreas is an organ featuring distinct endocrine and exocrine compartments. The exocrine pancreas, which comprises 98% of the pancreas (1), is composed primarily of acini that secrete enzymes involved in digestion into the ducts (2). The remaining 2% constitutes the endocrine pancreas, which comprises pancreatic islets, the main function of which is to secrete insulin (3). This distinction is also seen in diseases associated with each compartment; type 2 diabetes (T2D) results from the dysregulation of the endocrine pancreas (3), whereby the dysfunction of the exocrine pancreas induces pancreatitis (i.e., inflammation) or pancreatic cancer, an increasingly leading cause of death (4).
Epidemiologic studies have shown that individuals with T2D have an increased risk of developing pancreatic cancer (5), and the reverse is also true (6). In addition, there is an association between T2D and chronic pancreatitis, leading to acinar atrophy (7–9).
A recent study exploring the effects of T2D and age on the exocrine pancreas found that compared with individuals of controls without diabetes and patients with type 1 diabetes, pancreases of patients with T2D presented with greater acinar-to-ductal metaplasia (ADM), a process whereby exocrine acinar cells transdifferentiate into ductal-like cells (10). Furthermore, ADM correlated positively with T2D duration, and this duration correlated with pancreatic cancer initiation, validating the link between T2D and pancreatic cancer. In this respect, pancreatic ADM is considered a precancerous state of the pancreas. Understanding exocrine–endocrine functional interactions has important implications in the prevention of these deadly pancreatic diseases, because as yet, their molecular mechanisms remain largely unexplored.
We hypothesized that functional alterations in the exocrine pancreas are linked to local epigenetic changes in the context of environmental factors, which may exacerbate the risk of further exocrine disease.
Research Design and Methods
Methylation in Human Samples
We selected 155 whole-pancreas samples, collected based on the Innovative Medicines Initiative for Diabetes (IMIDIA) consortium (11), according to sample availability. Of these samples, 32 had T2D, based on clinical characteristics and tests, according to the American Diabetes Association guideline (2019). Sample collection was followed by next-of-kin informed consent, with the approval of the local ethics committees in Pisa and Hannover.
DNA was extracted from whole-pancreas samples using the NucleoSpin Tissue Kit (cat. no. T740952.50; Macherey-Nagel). Bisulphite conversion was performed in 800 ng DNA using the EZ DNA Methylation Kit (cat. no. 5001; Zymo Research), and DNA was subjected to Illumina 850K MethylationEPIC array. Array data were imported using the minfi R package (12). Quality control steps removed CpG probes that were located on sex chromosomes or near single nucleotide polymorphisms (SNPs), cross-hybridizing, or non-CG or had a detection threshold P value <0.01. Samples with <99% probes with a detection P value <0.01 were excluded. Probe-design biases and batch effects were normalized using R packages Enmix and SVA (13,14). Samples were removed for having a call rate <99%, and two sex-discordant samples were removed. After quality control, 746,912 probes and 141 samples remained (Table 1), comprising 69 females and 72 males. Population structure was evaluated by principal component analysis, with 1000 Genomes as reference (Supplementary Fig. 2).
Clinical characteristics of participants
Characteristic . | No diabetes (n = 113)* . | T2D (n = 28)* . | P† . |
---|---|---|---|
Sex | 0.016 | ||
Female | 61 (54) | 8 (29) | |
Male | 52 (46) | 20 (71) | |
Age, years | 67 (53, 76) | 74 (69, 79) | 0.004 |
BMI, kg/m2 | 25.0 (23.1, 27.1) | 26.0 (24.2, 28.2) | 0.13 |
T2D duration, years | — | 9 (6, 12) | — |
T2D treatment | 0 | 25 (89) | — |
Statin treatment | 10 (17) | 4 (29) | 0.4 |
Characteristic . | No diabetes (n = 113)* . | T2D (n = 28)* . | P† . |
---|---|---|---|
Sex | 0.016 | ||
Female | 61 (54) | 8 (29) | |
Male | 52 (46) | 20 (71) | |
Age, years | 67 (53, 76) | 74 (69, 79) | 0.004 |
BMI, kg/m2 | 25.0 (23.1, 27.1) | 26.0 (24.2, 28.2) | 0.13 |
T2D duration, years | — | 9 (6, 12) | — |
T2D treatment | 0 | 25 (89) | — |
Statin treatment | 10 (17) | 4 (29) | 0.4 |
Data given as n (%) or median (interquartile range).
Pearson χ2 test or Wilcoxon rank sum test.
Comparison With Pancreatic Islet Profiles
To compare the profile of exocrine pancreatic extractions with pancreatic islets, we used pancreatic islet data that were generated using the same methods as those described above and were also collected based on the IMIDIA consortium. The pancreatic islet samples were obtained from 144 organ donors, with an age range of 22 to 96 years (mean age 69; range 22–96 years). There was no statistically significant difference in age, sex, BMI, or T2D status between the exocrine and pancreatic islet donors. Characteristics of the donors are detailed in Supplementary Table 1. Principal component analysis of methylation of exocrine and pancreatic islets was determined from β values from both groups using the flashPCA package in R.
Epigenome-Wide Association Study
Cell composition was estimated using the R package RefFreeEWAS (15), which estimated eight cell types (Supplementary Fig. 1). We adjusted for J-1 cell types (where J is the number of cell types defined) and, to avoid collinearity, removed the cell type with the lower estimated proportion. A linear regression model was applied for T2D association, corrected for age, sex, BMI, and cellular composition. Bonferroni correction was applied for multiple testing. The study design is summarized in Fig. 1.
Overview of the study design. Whole-pancreas tissue was extracted from organ donors, and DNA methylation data were generated using the MethylationEPIC 850K array to perform epigenome-wide association studies (EWAS) of T2D in organ donors. Further genetic and functional validation was performed.
Overview of the study design. Whole-pancreas tissue was extracted from organ donors, and DNA methylation data were generated using the MethylationEPIC 850K array to perform epigenome-wide association studies (EWAS) of T2D in organ donors. Further genetic and functional validation was performed.
UK Biobank Used to Identify Rare Variant Associations in PNLIPRP1
Exome sequencing data from up to 191,000 participants in the UK Biobank were used to identify rare null variants (minor allele frequency [MAF] <1%) in PNLIPRP1. This research is part of UK Biobank research application 67575. We tested the association of null variants for several traits, including BMI, glucose, and lipid traits, using the MiST method (16), which tests rare variants in a single cluster (at the gene scale). Score π represents the mean effect of the cluster; τ represents the heterogeneous effect of the cluster. The overall P value tests the association between the set of variants and the trait of interest. For each trait, we adjusted for relevant covariates. We considered a trait to be significant if the P value associated with the direct burden effect of the cluster was significant (Pπ value [i.e., P < 0.05]) and the direction of effect for the variants was consistent, revealed by the absence of heterogeneity (τ P value [i.e., P > 0.05]).
Mendelian Randomization
To assess the direction of causality between traits of interest (i.e., T2D, LDL cholesterol [LDL-C], and CpG methylation), we performed bidirectional two-sample Mendelian randomization (MR). To obtain genetic associations for both T2D and LDL-C, we consulted previously published large-scale European genome-wide association studies of T2D (17) and LDL-C (18). All MR analyses were performed using the TwoSampleMR (version 0.5.7) R software package.
To obtain genetic associations between SNPs and DNA methylation, we genotyped 111 control samples from our organ donors (Illumina HumanOmni2.5 arrays) using the Illumina iScan (method detailed in Supplementary Methods). Methylation quantitative trait loci (mQTL) analysis was performed using QTLtools software (19), adjusting for age, sex, and BMI. Forward and reverse MRs are detailed in Supplementary Methods. The inverse variance weighted (IVW) method was applied to compute the causal effect estimate. The F statistic was used to verify the strength of instrument. The MR-Egger method was applied to investigate the presence of horizontal pleiotropy. Linkage disequilibrium was assessed with the ld_clump() function from the R software package ieugwasr (version 0.1.5).
The leave-one-out analysis was used to verify whether any variant was driving our findings. Heterogeneity was assessed using the Cochran Q test.
RNA Expression of PNLIPRP1 in Organ Donors
Eleven samples were processed (six controls and five individuals with T2D), matched for age, sex, and BMI. RNA was extracted using Trizol (cat. no. 15596-026; Thermo Fisher Scientific). A detailed description of the results is provided in Supplementary Methods.
Functional Characterization
All in vitro assays were performed using the rat AR42J acinar cell line (cat. no. CRL-1492; American Type Culture Collection). Methods are further detailed in Supplementary Methods.
High Glucose and Insulin Treatment
AR42J cells were treated with 20 mmol/L glucose or 100 nmol/L insulin, or both, for 72 h, and cells were harvested for RT–quantitative (qPCR).
Akt Response and Glucose Uptake
Akt response was tested using Western blot in AR42J cells stimulated with or without 200 nmol/L insulin for 1 h. Details of all antibodies in this study are listed in Supplementary Table 3. Glucose uptake was measured using the colorimetric glucose detection kit (cat. no. EIAGLUC; Invitrogen), following the manufacturer protocol.
siRNA Knockdown
AR42J cells were transfected with PNLIPRP1 siRNA or nontargeting control siRNA. RNA and protein were harvested after 48 or 72 h (qPCR and Western blotting method detailed in Supplementary Methods).
Immunohistochemistry
Human pancreatic tissue sections, from controls and individuals with T2D, were obtained by the INSERM UMR 1190 unit (University of Lille, Lille, France). We excluded any fibrotic process or pathologic alterations in the pancreatic tissues.
RNA Sequencing
Library preparation from PNLIPRP1 knockdown (KD) RNA was performed using the KAPA mRNA HyperPrep Preparation Kit (Roche) and sequenced with NovaSeq 6000 (Illumina) (method provided in Supplementary Methods). Differential expression analysis was determined using DeSeq2. Pathway analysis was performed using Metascape or EnrichR. MTS proliferation (cat. no. 197010; Acbam) and Cholesterol Ester-Glo Assay (cat. no. J3190; Promega) methods were performed per protocol instructions in PNLIPRP1 KD.
Data and Resource Availability
The data sets generated and/or analyzed during the current study are available in the Figshare repository (https://figshare.com/s/16ad75ace2b169a1623c). No applicable resources were generated or analyzed during the current study.
Results
DNA Methylation Changes in the Exocrine Pancreas Associated With T2D
We measured DNA methylation (MethylationEPIC arrays) in whole-pancreas samples obtained from 141 organ donors (age 17–89; median 67 years) of European descent (Table 1 and Supplementary Fig. 2). To validate that the observed whole-pancreas epigenetic profile was indeed mainly exocrine tissue, we compared the methylation profile with pancreatic islet methylation profiles (Infinium MethylationEPIC arrays) from 125 individuals from the organ donor cohort, matched for age, sex, and BMI (Supplementary Table 1). The methylation profiles of pancreatic islets and whole pancreas were distinct, and this was consistent with four pancreatic islet samples extracted and handled in parallel with the exocrine preparations (Supplementary Fig. 2). Pancreatic islet preparations with low islet purity (i.e., percentage of islets compared with other pancreatic tissues, such as exocrine) did not overlap with exocrine preparations (Supplementary Fig. 3).
We performed an epigenome-wide association study to determine DNA methylation changes associated with T2D. We detected a single CpG significantly associated with T2D: hypermethylation of cg15549216 CpG (Bonferroni-corrected P = 0.025) (Fig. 2A). This CpG was located in the gene body of PNLIPRP1, encoding pancreatic lipase–related protein 1, and was 11.4% more methylated in patients with T2D compared with controls (estimate 0.6; SE 0.1) (Fig. 2B). The cg15549216 CpG was also positively correlated with glucose levels (P = 1.34 × 10−4) (Supplementary Fig. 4). We detected one differentially methylated region associated with increased T2D risk, which encompassed this cg15549216 CpG, and two flanking CpGs, cg06606475 and cg08580014 (Fig. 2C and D), that were consistent in direction of effect: the cg06606475 CpG located 921 bp upstream of the cg15549216 probe (mean difference 9.0%; P = 5.9 × 10−5; estimate 0.38; SE 0.09) and the cg08580014 CpG located 370 bp downstream of the cg15549216 probe (mean difference 6.3%; P = 4.0 × 10−4; estimate 0.27; SE 0.07). The Genehancer and dbSUPER enhancer databases showed that in human whole-pancreas tissue, the region in and surrounding the cg15549216 probe lay in a 12-kb stretch of the only known super-enhancer of the PNLIPRP1 gene (20,21) (Fig. 2E).
Epigenome-wide association study in whole pancreas in T2D. A: Volcano plot depicting the differentially methylated CpGs associated with T2D, highlighting the only significant CpG (Bonferroni-corrected P < 0.05). Colors indicate the absolute estimates (M values). B: Boxplot showing the density of cg15549216 methylation levels in individuals with T2D (orange) and controls without diabetes (purple), for each individual. The x-axis represents the percentage methylation β value. B–D: Boxplots showing the density of cg06606475 and cg08580014 methylation levels, respectively, in patients with T2D and controls without diabetes, for each individual. E: Representation of the PNLIPRP1 gene, with the differentially methylated region shaded in orange. The enhancer region is depicted in purple, and the estimates and P values of all CpGs in the PNLIPRP1 gene are shown. F: Scatterplot depicting the interaction between methylation of cg15549216 and age in organ donor individuals. Controls without diabetes are shown in purple, and patients with T2D are shown in orange.
Epigenome-wide association study in whole pancreas in T2D. A: Volcano plot depicting the differentially methylated CpGs associated with T2D, highlighting the only significant CpG (Bonferroni-corrected P < 0.05). Colors indicate the absolute estimates (M values). B: Boxplot showing the density of cg15549216 methylation levels in individuals with T2D (orange) and controls without diabetes (purple), for each individual. The x-axis represents the percentage methylation β value. B–D: Boxplots showing the density of cg06606475 and cg08580014 methylation levels, respectively, in patients with T2D and controls without diabetes, for each individual. E: Representation of the PNLIPRP1 gene, with the differentially methylated region shaded in orange. The enhancer region is depicted in purple, and the estimates and P values of all CpGs in the PNLIPRP1 gene are shown. F: Scatterplot depicting the interaction between methylation of cg15549216 and age in organ donor individuals. Controls without diabetes are shown in purple, and patients with T2D are shown in orange.
Because aging is one of the major risk factors for T2D (22), we verified whether the cg15549216 CpG was associated with age in our cohort. We found that age was nominally associated with the hypermethylation of the cg15549216 probe (unadjusted P = 0.010); however, this association was not significant after adjustment for T2D status (unadjusted P = 0.93). We considered whether both T2D and age could be contributing factors in cg15549216 methylation by including an interaction term between T2D status and age in our model. At the cg15549216 probe, we found an increase in methylation in individuals with T2D as age increased compared with controls without diabetes (M value estimate 0.01; P = 4 × 10−10) (Fig. 2F). The cg15549216 probe was not associated with BMI (P = 0.54), sex (P = 0.32), T2D duration (P = 0.39), or statin treatment (P = 0.32).
Rare PNLIPRP1 Null Variants Are Associated With Metabolic Traits
To investigate the putative role of PNLIPRP1 in glucose or lipid metabolism, we tested whether rare variants (MAF <1%) in the PNLIPRP1 gene were associated with T2D and related metabolic traits, including glucose and lipid traits. Using whole-exome sequencing data from up to 191,000 participants in the UK Biobank, we identified 44 null variants (i.e., nonsense, frameshift, canonical ±1 or 2 splice sites) (Supplementary Table 4).
We found that PNLIPRP1 null variants were associated with increased glycemia (Pπ = 1.1 × 10−3; effect size 0.13; SE 0.040). We also found that PNLIPRP1 rare variants were significantly associated with increased LDL-C (Pπ = 0.034; effect size 0.10; SE 0.049), HDL-C (Pπ = 0.026; effect size 0.05; SE 0.022), waist-to-hip ratio (Pπ = 2.7 × 10−3; effect size 0.23; SE 0.006), waist circumference (Pπ = 7.6 × 103; effect size 2.19; SE 0.820), BMI (Pπ = 3.1 × 10−3; effect size 0.039; SE 0.011), diastolic blood pressure (Pπ = 2.02 × 10−5; effect size 2.89; SE 0.677), and systolic blood pressure (Pπ = 2.9 × 10−3; effect size 3.30; SE 1.112). Null variants were not associated with T2D risk (Pπ = 0.48), obesity (Pπ = 0.29), hypertension (Pπ = 0.94), or triglyceride levels (Pπ = 0.46) (Table 2).
PNLIPRP1 null variant associations
Trait . | n of individuals . | Π . | SE . | Pπ . | n of variants . |
---|---|---|---|---|---|
Diastolic BP | 168,374 | 2.889 | 0.677 | 2.0 × 10−5 | 39 |
WHR | 190,739 | 0.023 | 0.006 | 2.7 × 10−3 | 41 |
Log BMI | 187,727 | 0.039 | 0.011 | 3.1 × 10−3 | 40 |
Glucose | 159,764 | 0.133 | 0.041 | 1.1 × 10−3 | 38 |
Systolic BP | 168,367 | 3.303 | 1.1119 | 2.9 × 10−3 | 39 |
Waist circumference | 190,751 | 2.188 | 0.820 | 7.6 × 10−3 | 41 |
HDL | 156,077 | 0.050 | 0.022 | 2.6 × 10−2 | 38 |
LDL | 169,625 | 0.104 | 0.049 | 3.4 × 10−2 | 38 |
Obesity | 107,219 | −0.146 | 0.138 | 0.29 | 31 |
Log triglycerides | 169,815 | −0.009 | 0.011 | 0.46 | 39 |
T2D | 171,651 | 0.177 | 0.252 | 0.48 | 39 |
Hypertension | 180,290 | 0.010 | 0.126 | 0.94 | 40 |
Trait . | n of individuals . | Π . | SE . | Pπ . | n of variants . |
---|---|---|---|---|---|
Diastolic BP | 168,374 | 2.889 | 0.677 | 2.0 × 10−5 | 39 |
WHR | 190,739 | 0.023 | 0.006 | 2.7 × 10−3 | 41 |
Log BMI | 187,727 | 0.039 | 0.011 | 3.1 × 10−3 | 40 |
Glucose | 159,764 | 0.133 | 0.041 | 1.1 × 10−3 | 38 |
Systolic BP | 168,367 | 3.303 | 1.1119 | 2.9 × 10−3 | 39 |
Waist circumference | 190,751 | 2.188 | 0.820 | 7.6 × 10−3 | 41 |
HDL | 156,077 | 0.050 | 0.022 | 2.6 × 10−2 | 38 |
LDL | 169,625 | 0.104 | 0.049 | 3.4 × 10−2 | 38 |
Obesity | 107,219 | −0.146 | 0.138 | 0.29 | 31 |
Log triglycerides | 169,815 | −0.009 | 0.011 | 0.46 | 39 |
T2D | 171,651 | 0.177 | 0.252 | 0.48 | 39 |
Hypertension | 180,290 | 0.010 | 0.126 | 0.94 | 40 |
BP, blood pressure; WHR, waist-to-hip ratio.
PNLIPRP1 Common Variants Are Associated With LDL-C but Not T2D- or Glucose-Related Traits
Using the Type 2 Diabetes Knowledge Portal (23), we found that common SNPs (MAF >1%) located at the PNLIPRP1 locus were strongly associated with increased LDL-C (P = 2.0 × 10−14) (Supplementary Tables 5 and 6). In addition, common variants in PNLIPRP1 were associated with non-HDL cholesterol, cholesterol levels, and apolipoprotein B after multiple testing (P ≤ 2.5 × 10−6), but not with T2D- or glucose-related traits, in up to 1.61 million participants (Supplementary Fig. 5).
Causal Relationship Between PNLIPRP1 Methylation and T2D-Related Traits
To identify whether PNLIPRP1 methylation was causally associated with T2D and LDL-C, and the reverse, we performed MR. To identify SNPs acting as proxies for PNLIPRP1 methylation, we genotyped 111 control individuals in our cohort and performed mQTL. For T2D, we used 118 proxy SNPs associated with increased T2D risk and six proxy SNPs associated with PNLIPRP1 rmethylation (Supplementary Tables 7 and 8). Using the IVW method, we found evidence that increased T2D risk was causal of cg15549216 hypermethylation, with an estimate of 0.23 (95% CI 0.029–0.43; P = 0.025) (Fig. 3A). However, we did not find a significant association using MR methods following different assumptions, such as the simple median and weighted median methods, possibly pointing to the presence of horizontal pleiotropy or heterogeneity in our genetic instruments. To ensure the validity of our results, we used another robust MR method, MR-Egger, which showed no evidence of causal association (estimate 0.102; 95% CI −0.353 to 0.557; P = 0.659), and ruled out the possibility of horizontal pleiotropy (intercept 0.011; 95% CI −0.024 to 0.045; P = 0.55). In the reverse direction, we found no evidence of a causal effect of PNLIPRP1 methylation on T2D risk (Fig. 3A).
MR to determine causality. Bidirectional causality between T2D and PNLIPRP1 methylation (A) and between LDL-C and PNLIPRP1 methylation (B) using the IVW, simple median, weighted median, and MR-Egger methods. The x-axis represents the estimates. Statistically significant associations are shown in black text, and nonsignificant values are shown in gray, with estimates and 95% CIs in parentheses.
MR to determine causality. Bidirectional causality between T2D and PNLIPRP1 methylation (A) and between LDL-C and PNLIPRP1 methylation (B) using the IVW, simple median, weighted median, and MR-Egger methods. The x-axis represents the estimates. Statistically significant associations are shown in black text, and nonsignificant values are shown in gray, with estimates and 95% CIs in parentheses.
For LDL-C, we used nine proxy SNPs to represent PNLIPRP1 methylation and 241 proxy SNPs for increased LDL-C levels (Supplementary Tables 9 and 10). We found no evidence of a causal effect of LDL-C on PNLIPRP1 methylation; however, PNLIPRP1 methylation was associated with increased LDL-C levels, with an estimate of 0.064 (95% CI 0.024–0.105; P = 0.0019) (Fig. 3B). This result was consistent with simple median (estimate 0.071; 95% CI 0.043–0.099; P = 5.9 × 10−7) and weighted median methods (estimate 0.068; 95% CI 0.044–0.091; P = 2.2 × 10−8). The MR-Egger intercept showed no evidence of horizontal pleiotropy (estimate 0.000; 95% CI −0.012 to 0.012; P = 0.31) (Fig. 3B). In all significant MR analyses, the leave-one-out analysis found that no single SNP altered the results, suggesting that the observed association was robust (Supplementary Tables 11 and 12). Altogether, our data provide some evidence that PNLIPRP1 hypermethylation increases LDL-C levels and suggest a trend that T2D status increases cg15549216 methylation.
Correlation of cg15549216 CpG With PNLIPRP1 Gene Expression in Whole Pancreas
To determine whether the expression of PNLIPRP1 was dysregulated in T2D, we performed gene expression quantification of PNLIPRP1 in extracted RNA from whole-pancreas tissue in a subset of five donors with T2D and six controls without diabetes, who were matched for age, sex, and BMI (Supplementary Table 13). We found that PNLIPRP1 expression was downregulated in donors with T2D (53%; P = 0.011) (Fig. 4A). PNLIPRP1 expression was not associated with age (P = 0.17). We performed linear regression to determine whether methylation levels were associated with PNLIPRP1 gene expression. We found that cg15549216 methylation was significantly correlated with decreased PNLIPRP1 gene expression (P = 0.042; R2 = 0.35) (Fig. 4B).
Expression of PNLIPRP1 at the RNA and protein levels. A: RNA expression of PNLIPRP1 in a subset of 11 individuals (five with T2D vs. six controls without diabetes), matched for age, sex, and BMI. t tests were performed to test for differences between T2D and nondiabetic samples for the PNLIPRP1 gene. Error bars represent SEs. B: Correlation between the cg15549216 probe and RNA expression of the PNLIPRP1 gene determined by qPCR and compared with housekeeping gene RPLP0. Dots shown in black are donors without diabetes, and those in purple are individuals with T2D. C: Immunofluorescence of healthy human whole-pancreas tissue samples stained for PNLIPRP1 and KRT19, a ductal marker (top), and insulin (INS; a marker for pancreatic islets). Nuclei were stained with DAPI.
Expression of PNLIPRP1 at the RNA and protein levels. A: RNA expression of PNLIPRP1 in a subset of 11 individuals (five with T2D vs. six controls without diabetes), matched for age, sex, and BMI. t tests were performed to test for differences between T2D and nondiabetic samples for the PNLIPRP1 gene. Error bars represent SEs. B: Correlation between the cg15549216 probe and RNA expression of the PNLIPRP1 gene determined by qPCR and compared with housekeeping gene RPLP0. Dots shown in black are donors without diabetes, and those in purple are individuals with T2D. C: Immunofluorescence of healthy human whole-pancreas tissue samples stained for PNLIPRP1 and KRT19, a ductal marker (top), and insulin (INS; a marker for pancreatic islets). Nuclei were stained with DAPI.
PNLIPRP1 Expression in Whole Pancreas Is Restricted to Acinar Cells
To identify target cell types of PNLIPRP1 dysregulation, we examined the expression of PNLIPRP1. In the GTEx database (24), PNLIPRP1 was exclusively expressed in the whole pancreas, compared with 53 tissues (Supplementary Fig. 6), and the TIGER database (25) showed that PNLIPPR1 was mainly expressed in the whole pancreas (median transcript per million [TPM] 2,581), compared with pancreatic islets (median TPM 36) (Supplementary Fig. 7). We analyzed several single-cell pancreas data sets (26,27) and found PNLIPRP1 was mainly expressed in acinar cells (Supplementary Table 14). To confirm at the protein level the acinar expression of PNLIPRP1, we performed immunofluorescence in human pancreatic tissue stained for PNLIPRP1, KRT19, a marker of ductal cells, and insulin to mark pancreatic islets. This revealed that the PNLIPRP1 protein was not expressed in pancreatic islets or ductal cells but exclusively in acinar cells (Fig. 4C).
PNLIPRP1 KD Induces Decreased Cell Cycle and Increases Cholesterol Biosynthesis Genes
Based on the acinar-specific expression pattern of PNLIPRP1, we explored the functional role of PNLIPRP1 in the rat acinar cell line AR42J. To investigate the downstream consequences of PNLIPRP1 dysregulation, we performed KD of PNLIPRP1 in the AR42J cell line. We confirmed that our PNLIPRP1 KD induced a 60% reduction in gene expression of PNLIPRP1 (Fig. 5A) and a 34% decrease at the protein level (Supplementary Fig. 8). RNA sequencing confirmed that PNLIPRP1 was one of the most significant downregulated genes (false discovery rate 0.025). We explored potential dysregulated pathways by focusing on the 1,024 genes with an unadjusted P < 0.05 (Supplementary Table 15). Using enrichR, we found that the top pathway for the downregulated genes (587 genes) was the cell-cycle pathway (adjusted P = 2.21 × 10−12) (Supplementary Table 16), and the most upregulated pathways (437 genes) were cholesterol biosynthesis and SREBP control of lipid biosynthesis (adjusted P = 0.0044) (Supplementary Table 17 and Fig. 5B). We confirmed using MTS proliferation that PNLIPRP1 KD led to a 22% reduction in cell proliferation (P = 0.0027) (Fig. 5C). We also confirmed that total cholesterol content was increased by 29% after PNLIPRP1 KD (P = 0.0017) (Fig. 5D). We then found that treatment of PNLIPRP1 KD with simvastatin (i.e., a statin drug that inhibits the cholesterol synthesis) was able to rescue the increase in cholesterol biosynthesis to levels of the control (Fig. 5E). Taken together, these findings demonstrate that the invalidation of PNLIPRP1 dysregulates in an inverse direction cholesterol metabolism and cell-cycle regulation.
PNLIPRP1 KD in AR42J rat cell line. A: Confirmation of PNLIPRP1 KD in four biologic replicates using qPCR and tested using a t test. B: Volcano plot of RNA sequencing results, depicting PNLIPRP1 (black) and downregulated cell-cycle genes (blue). The top five most significant genes are labeled. Upregulated genes from cholesterol synthesis and SREBP control of lipid biosynthesis genes (red and dark red) and ADM genes (green) are labeled. C: MTS proliferation assay in AR42J in PNLIPRP1 KD after 72 h of siRNA treatment. D: Total cholesterol measured after 48 h of PNLIPRP1 KD. All performed in three biologic replicates. E: Statin treatment after PNLIPRP1 KD, compared with untreated DMSO controls. Performed in two biologic replicates. Error bars represent SEs. F: Gene expression of acinar and ductal markers using qPCR in AR42J in PNLIPRP1 KD after 72 h of siRNA treatment. P values are from t tests performed between each gene in the siRNA of PNLIPRP1 (siPnliprp1) and untargeted control (siControl). The y-axis shows the fold change, and error bars represent SEs. Experiments were performed in four biologic replicates. *P < 0.05, ***P < 0.001. ns, not significant.
PNLIPRP1 KD in AR42J rat cell line. A: Confirmation of PNLIPRP1 KD in four biologic replicates using qPCR and tested using a t test. B: Volcano plot of RNA sequencing results, depicting PNLIPRP1 (black) and downregulated cell-cycle genes (blue). The top five most significant genes are labeled. Upregulated genes from cholesterol synthesis and SREBP control of lipid biosynthesis genes (red and dark red) and ADM genes (green) are labeled. C: MTS proliferation assay in AR42J in PNLIPRP1 KD after 72 h of siRNA treatment. D: Total cholesterol measured after 48 h of PNLIPRP1 KD. All performed in three biologic replicates. E: Statin treatment after PNLIPRP1 KD, compared with untreated DMSO controls. Performed in two biologic replicates. Error bars represent SEs. F: Gene expression of acinar and ductal markers using qPCR in AR42J in PNLIPRP1 KD after 72 h of siRNA treatment. P values are from t tests performed between each gene in the siRNA of PNLIPRP1 (siPnliprp1) and untargeted control (siControl). The y-axis shows the fold change, and error bars represent SEs. Experiments were performed in four biologic replicates. *P < 0.05, ***P < 0.001. ns, not significant.
An increase in cholesterol content and a reduction in proliferation are indicative of ADM, a response to acinar cells stress (28,29). Therefore, we assessed the expression of acinar and ductal markers in our RNA sequencing data, data subsequently confirmed by qPCR. We found that the expression of four acinar markers was downregulated: Prss1 (P = 0.0066), Amy2 (P = 0.0029), Cpa2 (P = 0.035), and Ctrl (Ctrl was downregulated in three of four replicates [P = 0.12]). Two ductal markers were upregulated: Krt19 (P = 0.0010) and Hnf1b (P = 0.020) (Fig. 5F). Taken together, these results suggest that PNLIPRP1 invalidation induces ADM.
Reduction of PNLIPRP1 Expression in Response to Diabetogenic Exposure
To address whether T2D induces dysregulation of PNLIPRP1 and mimics the phenotype observed in PNLIPRP1 KD, we assessed the impact of a diabetogenic environment (i.e., high glucose and high insulin) on AR42J cells. We confirmed that AR42J cells were responsive to insulin treatment by analyzing phosphorylated Akt (Supplementary Fig. 9). We found that high glucose and insulin induced a decrease in PNLIPRP1 expression (Fig. 6A). We tested whether the treatment induced glucose uptake and found it decreased glucose in the medium of treated cells (P < 0.0001) (Fig. 6B), indicating increased glucose uptake by the cells. We confirmed using the MTS proliferation assay that the diabetogenic milieu led to a 25% reduction in cell proliferation (P = 0.0089) (Fig. 6C). Total cholesterol content was unchanged (P = 0.1002) (Fig. 6D). We also assessed acinar and ductal markers and found that the expression of two acinar markers was downregulated (Cpa2 [P = 0.0281] and Ctrl [P = 0.0022]), whereas Prss1 and Amy2 (P > 0.05) were unchanged. Ductal markers were significantly upregulated: Krt19 (P = 0.0140) and Hnf1b (P = 0.0344) (Fig. 6E). Taken together, these results suggest that high glucose and insulin treatment induced downregulation of PNLIPRP1, a decrease in proliferation, and ADM in these cells. We performed immunofluorescence in a T2D sample, compared with a healthy control, and did not find any notable differences in PNLIPRP1 expression, but we did observe that some cells expressed both ductal marker KRT19 and PNLIPRP1 expression, suggesting an ADM phenotype (Supplementary Fig. 10).
High-glucose (HG) and insulin treatment of AR42J rat cell line. A: RNA expression of PNLIPRP1 after HG and insulin treatment, determined by qPCR; statistical significance was determined using a two-tailed t test. Data are representative of four biologic replicates. B: Glucose uptake of PNLIPRP1 after HG (20 mmol/L glucose) and insulin (100 nmol/L) treatment, compared with untreated controls, determined by qPCR; statistical significance was determined using a two-tailed t test. Data are representative of six biologic replicates. C: MTS proliferation assay in 48 h AR42J treated with HG and insulin, compared with controls. D: Total cholesterol measured after 48 h of HG and insulin, compared with controls. E: Gene expression of acinar and ductal markers in HG- and insulin-treated cells, compared with controls, in AR42J. P values are from t tests performed for each gene. The y-axis shows the fold change, and error bars represent SEs. Experiments were performed in three biologic replicates.
High-glucose (HG) and insulin treatment of AR42J rat cell line. A: RNA expression of PNLIPRP1 after HG and insulin treatment, determined by qPCR; statistical significance was determined using a two-tailed t test. Data are representative of four biologic replicates. B: Glucose uptake of PNLIPRP1 after HG (20 mmol/L glucose) and insulin (100 nmol/L) treatment, compared with untreated controls, determined by qPCR; statistical significance was determined using a two-tailed t test. Data are representative of six biologic replicates. C: MTS proliferation assay in 48 h AR42J treated with HG and insulin, compared with controls. D: Total cholesterol measured after 48 h of HG and insulin, compared with controls. E: Gene expression of acinar and ductal markers in HG- and insulin-treated cells, compared with controls, in AR42J. P values are from t tests performed for each gene. The y-axis shows the fold change, and error bars represent SEs. Experiments were performed in three biologic replicates.
Discussion
To our knowledge, this is the first epigenetic study of the impact of T2D on the whole pancreas. We found a strong association between T2D and hypermethylation of the cg15549216 probe, located in the only known enhancer of the PNLIPRP1 gene.
Additional observations revealed that PNLIPRP1 is a strong candidate for T2D-driven epigenetic changes, leading to deleterious events, particularly in acinar cell cholesterol metabolism. This is supported by several lines of evidence: 1) cg15549216 hypermethylation was causal of increased LDL-C, 2) partial extinction of PNLIPRP1 increased cholesterol metabolism in functional studies, and 3) this cholesterol impairment was rescued by statin treatment in vitro. Furthermore, our PNLIPRP1 extinction studies also showed that acinar cells started to express ductal markers and cell-cycle arrest, which both suggest ADM. Cholesterol synthesis and ADM are functionally linked; several studies have shown that cholesterol biosynthesis triggers ADM, and this is reversed in vitro through statin (28).
Although ADM should in theory be reversible, during sustained stress these cells are strongly associated with progression to pancreatic ductal adenocarcinoma (PDAC), the most common malignancy of the exocrine pancreas (30–33). A recent study followed the time course of mouse PDAC using single-cell transcriptomics and revealed that after ADM, metaplastic cells had a signature distinct from that of ductal cells and that PNLIPRP1 was significantly downregulated in these metaplastic cells compared with acinar cells (34). Additionally, PNLIPRP1 is one of the most significantly downregulated genes in PDAC tumors (35), and there are somatic deletions in PNLIPRP1 in PDAC tissue (36), pointing to a role of this gene in pancreatic reprogramming toward malignancy. Additionally, in vitro and in vivo models have shown that statins delayed the formation of PDAC in mutated mice (37,38).
This is reinforced by human studies, where statin use resulted in a significant reduction in pancreatic cancer risk and improved patient survival (39). For instance, a meta-analysis that included 14 studies showed that statin use was associated with an overall reduction in death in individuals with PDAC, particularly significant in individuals who had undergone surgery, but not in those with advanced disease (40). Indeed, in patients with advanced cancer, the late addition of statins to standard anti-cancer therapy did not show improvement in overall survival.These findings reveal the protective properties of statins, especially when administered early and before disease development, highlighting the importance of studying the early events in pancreatic injury (41).
We found that rare and common genetic variants were associated with increased LDL-C and metabolic traits. Specifically, rare variants were associated with LDL-C, HDL-C, BMI, glucose levels, waist-to-hip ratio, waist circumference, diastolic blood pressure, and systolic blood pressure, indicating a systemic effect on metabolic disease. However, we did not find an association with disease states, including T2D, obesity, or hypertension. This could be due to the involvement of PNLIPRP1 in intermediate metabolic pathways or disease progression stages. This could also be due to limitations in statistical power. Further research with larger cohorts and more detailed phenotypic data may help clarify these relationships. Our null variant data are supported by PNLIPRP1 knockout mice, which display mild hyperglycemia and impaired insulin sensitivity, exacerbated by a high-fat diet (42). Additionally, evolutionary studies have shown that PNLIPRP1 was consistently lost in herbivore or carnivore species consuming a low-fat diet (43), suggesting a role in lipid metabolism. Taken together, our study data suggest that both nature (genetic variation) and nurture (through epigenetic marks) contribute to the regulation of PNLIPRP1 and various aspects of systemic metabolic health.
This study has several important limitations to consider. First, the sample size was relatively small, underscoring the need for larger exocrine pancreas biobanks to better explore and understand the DNA methylation landscape in this tissue. Second, functional experiments were conducted using the AR42J rat cell line because of challenges associated with working with human exocrine tissue. Because the methylation region is not conserved in rat cells, further research is required to directly link DNA methylation with PNLIPRP1 gene expression in human cells. Finally, although our findings indicate an association between DNA methylation and gene expression, we lacked sufficient tissue to investigate the impact of DNA methylation on protein expression. Addressing this link is crucial for a more complete understanding.
Understanding the sequence of events leading to pancreatic disease in humans is challenging and necessitates a comprehensive and integrative approach. Our findings not only underscore the importance of understanding the underlying molecular mechanisms relating T2D and pancreatic disease but also suggest preventive strategies (e.g., statin or other lipid metabolism modifiers to protect exocrine tissue). Despite the challenges inherent in studying the exocrine pancreas, unraveling the molecular intricacies linking T2D and pancreatic diseases remains crucial for advancing our understanding of these conditions.
This article contains supplementary material online at https://doi.org/10.2337/figshare.26510815.
Article Information
Funding. This study was funded by the Société Francophone du Diabète (SFD-R21100EE/RAK21020EEA) and French National Research Agency (Agence Nationale de la Recherche [ANR]–10-LABX-46 [European Genomics Institute for Diabetes] and ANR-10-EQPX-07-01 [Lille Integrated Genomics Advanced Network for Personalized Medicine]), the EU Horizon Europe Research and Innovation Programme (OBELISK grant agreement 101080465), the France Génomique consortium (ANR-10-INBS-009). This work was supported by the National Center for Precision Diabetic Medicine – PreciDIAB (ANR-18-IBHU-0001; FEDER 20001891/NP0025517; MEL 2019_ESR_11). A.B. is supported by the European Research Council (OpiO 101043671).
Duality of Interest. No potential conflicts of interest relevant to this article were reported.
Author Contributions. L.Mau., M.B., L.N., M.C., and A.K. performed the statistical analyses and generated the figures. L.Mau., R.B., A.L., and J.F. performed the functional work. L.Mau and V.P. performed the MR analyses. L.Mar., C.D.L., A.J., S.L., M.S., and P.M. provided organ donor samples and clinical data. B.T. and S.A. performed RNA sequencing. M.D. and A.B. performed UK Biobank analyses. F.P. and J.K.-C. provided pancreatic tissue samples. P.F. and A.K. designed the study. All authors critically reviewed and edited the manuscript. A.K. is the guarantor of this work and, as such, had full access to all the data in the study and takes responsibility for the integrity of the data and the accuracy of the data analysis.
Prior Presentation. Parts of this study were presented in abstract form at the 56th Annual Meeting of the European Association for the Study of Diabetes (EASD), virtual conference, 22–25 September 2020; 8th Meeting of the Study Group on the Genetics of Diabetes. Lille, France, 11–13 May 2022; Société Francophone de Diabète (SFD), Montpellier, France, 21–24 March 2023; SFD, Toulouse, France, 19–22 March 2024; and 59th Annual Meeting of the EASD, Hamburg, Germany, 3–6 October 2023.