A functional variant rs912304 for late-onset T1D risk contributes to islet dysfunction by regulating proinflammatory cytokine-responsive gene STXBP6 expression

Abstract

Background

Our previous genome‑wide association studies (GWAS) have suggested rs912304 in 14q12 as a suggestive risk variant for type 1 diabetes (T1D). However, the association between this risk region and T1D subgroups and related clinical risk features, the underlying causal functional variant(s), putative candidate gene(s), and related mechanisms are yet to be elucidated.

Methods

We assessed the association between variant rs912304 and T1D, as well as islet autoimmunity and islet function, stratified by the diagnosed age of 12. We used epigenome bioinformatics analyses, dual luciferase reporter assays, and expression quantitative trait loci (eQTL) analyses to prioritize the most likely functional variant and potential causal gene. We also performed functional experiments to evaluate the role of the causal gene on islet function and its related mechanisms.

Results

We identified rs912304 as a risk variant for T1D subgroups with diagnosed age ≥ 12 but not < 12. This variant is associated with residual islet function but not islet-specific autoantibody positivity in T1D individuals. Bioinformatics analysis indicated that rs912304 is a functional variant exhibiting spatial overlaps with enhancer active histone marks (H3K27ac and H3K4me1) and open chromatin status (ATAC-seq) in the human pancreas and islet tissues. Luciferase reporter gene assays and eQTL analyses demonstrated that the biallelic sites of rs912304 had differential allele-specific enhancer activity in beta cell lines and regulated STXBP6 expression, which was defined as the most putative causal gene based on Open Targets Genetics, GTEx v8 and Tiger database. Moreover, Stxbp6 was upregulated by T1D-related proinflammatory cytokines but not high glucose/fat. Notably, Stxbp6 over-expressed INS-1E cells exhibited decreasing insulin secretion and increasing cell apoptosis through Glut1 and Gadd45β, respectively.

Conclusions

This study expanded the genomic landscape regarding late-onset T1D risk and supported islet function mechanistically connected to T1D pathogenesis.

Background

Type 1 diabetes (T1D) is an autoimmune disease caused by the destruction of pancreatic β cells and the absolute insufficiency of insulin [1], which has been increasing worldwide for several decades [2]. T1D is a multifactorial disease, with immune, genetic, viral, and environmental factors contributing to its development [3]. As the initial trigger of T1D [4, 5], over 80 non-HLA regions have been reported to be associated with T1D risk [5,6,7] except for the HLA risk region (especially HLA-DR-DQ on 6p21) [5]. Studies also indicate that multiple non-HLA risk regions influence the risk for islet autoimmunity [8,9,10]. As an autoimmune disease, most T1D candidate genes were functionally related to immune responses, such as CTLA4, IL-2RA, PTPN2, and CD226 [5]. However, cumulative evidence supports a role for islet β-cell function in T1D pathogenesis [11]. Several candidate genes in T1D non-HLA risk regions, such as INS, GLIS3, CTSH, and CLEC16A, were reported to affect islet function [12]. Our previous study [13] also indicated that the genetic risk score (GRS) affected fasting C-peptide levels in newly diagnosed T1D individuals.

Syntaxin binding protein 6 (STXBP6), the sixth member of the syntaxin-binding protein family, is located in 14q12, which is a vital part of the SNAP receptor (SNARE) complex and is crucial in neuronal vesicle trafficking [14]. A prior study has reported that SOX4 affected extracellular insulin delivery by regulating STXBP6 expression [15, 16], which is related to glucose-stimulated insulin secretion by participating in the expansion of fusion pores during the insulin exocytosis process [15, 16]. Vujkovic et al. demonstrated that 14q12 was a T2D risk region, with the leading GWAS-significant variant rs11159347 near STXBP6 [17]. Scanning from our previous T1D GWAS [13], rs912304, another independent leading variant in STXBP6, was found to be associated with T1D risk. However, it remains to be determined whether rs912304 or its high linkage disequilibrium (LD) variant is functionally related to T1D and its related mechanism.

Therefore, this study investigated the association between the variant rs912304 and T1D subgroups and related clinical indicators, including islet autoimmunity and islet function. We further assessed the functional variant(s) and implacted gene(s) in 14q12. Finally, we evaluated the role of the putative causal gene on islet function and its related mechanism.

Methods

Study participants

A total of 2423 T1D and 5082 subjects were enrolled in the study. As previously described [13], we recruited 1045 T1D and 1308 healthy individuals for GWAS scanning, and the validation subjects consisted of 1378 patients and 3774 controls recruited from southeastern, central, and southern China. Only patients diagnosed with insulin-dependent diabetes mellitus within six months and having at least one positive autoantibody were included in the T1D subjects. Controls were outpatients from the same region with no diabetes and no family history of diabetes. All samples were collected with appropriate informed consent from all participants or guardians. The study was approved by the Ethics Committee of the First Affiliated Hospital of Nanjing Medical University and was conducted following the principles of the Declaration of Helsinki.

Laboratory measurements for glycemic traits and islet‑specific autoantibodies

A total of 2353 non-diabetic controls from the validation stage measured plasma glucose and serum insulin levels at fasting, 30 min, and 120 min based on an oral glucose tolerance test (OGTT). The plasma glucose and serum insulin levels were measured using a fully automated biochemical analyzer (Roche, Switzerland) and insulin radioimmunoassay kit (BMIBT, China). Insulinogenic Index (IGI) and corrected insulin response (CIR) were used as indicators of insulin release stimulated by oral glucose, while insulin sensitivity was evaluated using ISI Matsuda and other indices. Next, a mixed meal glucose tolerance test (MMTT) was performed on the 202 newly diagnosed T1D individuals. This test was not required for T1D patients with fasting C-peptide levels ≤ 25 pmol/L, but only the fasting C-peptide level was measured by chemiluminescence (Roche Diagnostics, Switzerland). In addition, islet-specific autoantibodies, including ZnT8A, GADA, and IA-2A, were measured when T1D subjects enrolled, which was detected separately by radioimmunoassay, as previously described [13].

Genotyping assay

Genomic DNA was extracted from peripheral blood lymphocytes using the QIAamp DNA blood extraction kit (QIAGEN, Germany). Genotyping of rs912304 in STXBP6 was from two aspects: 1045 T1D and 1308 healthy individuals were extracted for this variant in the T1D GWAS screening section [13]. The validated 1378 T1D and 3774 subjects were genotyped using TaqMan assays (Applied Biosystems). The distributions of genotypes in healthy and T1D subjects were in accordance with the Hardy–Weinberg equilibrium (P > 0.05).

In silico annotation of the functional variant and casual gene

For the potential functional variants, we performed with RegulomeDB v.2 [18], a tool incorporating genomic assays (TF ChIP-seq and DNase-seq from the ENCODE database and those overlapping the footprints and QTL data), to prioritize functional variant(s) located within 14q12. We also verified the results of functional predictions using HaploReg v4.2 [19], which utilized epigenomics data (chromatin states and histone H3 marks) primarily from Roadmap Epigenomics and ENCODE projects.

To prioritize the most causal gene in 14q12, we first annotated using our Variant-to-Gene (V2G) pipeline using Open Targets Genetics [20], which integrates molecular phenotype quantitative trait loci experiments, chromatin interaction experiments, variant effect predictor (VEP) from Ensembl, the distance between the variant and each gene’s canonical transcription start site (TSS), etc. We further identified the top-ranked eQTL gene using a meta-analysis P value of multi-tissues from GTEx v8 [21], confirmed by Tiger, Translational human pancreatic Islet Genotype tissue-Expression Resource [22].

Dual-luciferase reporter gene assay

Dual luciferase reporter gene assays were performed to assess whether the double-stranded oligonucleotide containing either G or T allele at the variant rs912304 showed different enhancer activity. DNA fragments covering the G or T allele at rs912304 (hg38 Chr14:25,006,964) were amplified, and constructed into the pGL3-promoter plasmids. The pGL3-promoter plasmid with DNA fragment containing either G or T allele, or empty pGL3-promoter vector in equal amounts was transiently co-transfected into HEK-239 T, INS-1E and MIN6 cell lines with pRL-TK plasmid using Lipofectamine 3000 transfection reagent (Invitrogen). Twenty-four hours post-transfection, cells were collected, and luciferase activity was measured by Dual-Luciferase Reporter Assay System (Promega). Firefly luciferase activity was expressed as relative luciferase units (RLU) after correction for Renilla luciferase activity to adjust for transfection efficiency. Data were normalized to those cells transfected with empty pGL3-Promoter vector. All assays were performed with a minimum of four replications in at least three independent experiments.

Cell culture and treatment

Isolation and culture of primary mouse islets, mouse β-cell lines (MIN6 and β-TC6), and HEK-293 T cells were performed as previously described [23]. The rat β-cell line INS-1E (donated by Prof. Xiao Han from Nanjing Medical University) was cultured in RPMI-1640 medium (11 mM glucose) containing 10% FBS, sodium pyruvate (1 mmol/L), HEPES (10 mmol/L), penicillin (100 U/mL), streptavidin (100 μg/mL), β-mercaptoethanol (50 μmol/L), and L-glutamine (2 mmol/L).

T1D-related cytokines consisted of 5 ng/ml IL-1β, 100 ng/ml INF-γ and 10 ng/ml TNF-α, respectively. A mixture of the three cytokine mixes was added to the complete medium of MIN6 cells, INS-1E cells, and β-TC-6 cells, respectively, according to the amount needed. The cells were spread in the culture plate at a density of about 1 × 10⁶ cells/ml and then incubated in a constant temperature and humidity incubator (containing 5% CO₂) at 37 °C. Upon reaching approximately 70% confluency, the cell culture medium was replaced with cytokines or free medium. In addition, high glucose (20 mM glucose) and high FFA (0.5 mM palmitic acid) were treated in INS-1E with a similar procedure. RNA and protein samples were accordingly collected after 48 h of incubation.

siRNA transfection and recombinant adenovirus infusion

MIN6 cells were seeded in 24-well plates and cultured until they reached 70% confluence. Transient transfection with siRNA targeting the Stxbp6 gene was done using Lipo3000 (Life Technologies, USA). All assays were performed 48 h after transfection. INS-1E cells were seeded in 6-well or 24-well plates and cultured until they reached 70% confluence. The cell culture medium was replaced the next day by one containing empty or Stxbp6 recombinant adenovirus (AdMax system), the amount of virus added per ml of cell culture medium = MOI × cell numers × 10³/virus titer), and the medium was changed after 6 h of transfection according to the cell status. After 48 h of incubation, the samples were collected, and the corresponding RNA and protein were extracted. SiRNA and recombinant adenovirus were purchased from Shanghai GenePharma CO. LTD and Shanghai OBio Tech CO. LTD, respectively.

GSIS functional assay

Transfected INS-1E were grown in 24-well plates and treated for 48 h accordingly. In brief, we added glucose to Kreb’s buffer to make 2 mM and 20 mM glucose solution and then discarded the solution after incubating them with glucose-free Kreb’s buffer for 60 min. Next, we added 500 μl 2 mM and 20 mM glucose solution per well and incubated for 60 min, collected the supernatant, centrifuged at 4 ℃, 3000 rpm, 15 min, and aspirated the supernatant, and stored at – 80 ℃. Finally, the cell culture supernatant was measured using insulin ELISA kits (MS300, Ezassay, China).

TUNEL assay

For quantifying cell apoptosis after adenovirus transfection, TUNEL staining was performed with YF594 TUNEL Assay Apoptosis Detection Kit (US Everbright, China) in INS-1E cells according to the instructions.

Bulk RNA sequencing

RNA-seq methods have been previously described [23]. Briefly, total RNA was extracted from INS-1E using TRIzol reagent (Invitrogen, CA, USA) according to the manufacturer's instructions. On-chip electrophoresis was used to measure the quality and quantity of extracted RNA using the Agilent RNA 6000 Nano Kit and Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA). Library preparation was performed on samples exhibited 1.9 ≤ A260/A280 ≤ 2.2, RNA integrity number > 8.0, and 28S/18S > 1.0. By enriching poly (A) mRNA with magnetic oligo (dT) beads, total RNA was purified. With the random hexamer (N6) primers, double-stranded cDNA fragments were synthesized and subjected to end repair and adapter ligation. Next, PCR amplification was used to construct cDNA libraries sequenced on BGISEQ-500 by BGI CO. LTD (Wuhan, China). Gene expression was quantified using RSEM 1.1.12 software and normalized using the method of fragments per kilobase of exon model per million reads mapped (FPKM). DEseq2 (logFC ≥ 0.5 and adjusted p ≤ 0.05) was used to identify differentially expressed genes (DEGs) between the Stxbp6 over-expressed and control groups.

q‑PCR

Total mRNA was extracted from tissues and cells using TRIzol reagent (Life Technologies). RNA from each sample was reverse transcribed using PrimeScript™ RT reagent kit (Takara Bio, Kusatsu, Japan). qRT-PCR was performed with SYBR Green PCR Master Mix (Takara Biotechnology) and analyzed using a real-time PCR system (StepOneplus, ABI). The relative gene expression was calculated using the comparative CT (2 ^− ΔΔCT) method with GAPDH as a reference gene. The primer sequences used are listed in Additional file 1: Table S1.

Western blotting

Cell lysates were prepared using RIPA lysis buffer on ice. The protein concentration was measured using the bicinchoninic acid (BCA) Assay Kit (Beyotime). WB was performed using the indicated antibodies (STXBP6 Polyclonal antibody, 10,976–4-AP, Proteintech; β-tubulin Polyclonal antibody, 10,068–1-AP, Proteintech; Caspase 3 Polyclonal antibody, 19,677–1-AP, Proteintech; Cleaved Caspase 3 Polyclonal antibody, 25,128–1-AP, Proteintech; ENO2 Monoclonal antibody, 66,150–1-Ig, Proteintech; GLUT1 Polyclonal antibody, 21,829–1-AP, Proteintech; GADD45B Antibody, DF2375, Affinity). Protein blots were visualized using SuperSignal West Pico Chemiluminescent Substrate (Pierce Biotechnology) in a ChemiDoc XRS + system (Bio-Rad Laboratories).

Statistical analysis

The frequencies of allele and genotype were assessed by direct counting. Associations between the variant and T1D risk in the case–control study and departures from the Hardy–Weinberg equilibrium (HWE) were analyzed by the χ² test. After testing for heterogeneity, a combined meta-analysis was performed using the Mantel–Haenszel (fixed-effects) and DerSimonian and Laird (random-effects) tests. The association of different genotypes with T1D risk was estimated using odds ratio (OR) tests with 95% CI. The heterogeneity between T1D subgroups was performed using Cochran’s Q test, considered significant with I² ≥ 75% and P ≤ 0.05. The additive model examined associations for glycemic quantitative traits using logistic regression analysis. Values of serum insulin, HOMA-B, IGI, CIR, HOMA-IR, and ISI_Matsuda were logarithmically transformed. The student's unpaired two-tailed t-test or Mann–Whitney test was used for pairwise comparisons depending on whether the data followed a Gaussian distribution. One-way ANOVA, or the Kruskal–Wallis test, was used for more than two groups. Error bars represent mean ± SD. P values were considered statistically significant with a P value of < 0.05. These statistical analyses were performed using SPSS 18.0, Stata11.0, and GraphPad Prism 9.0 software.

Results

Rs912304 in 14q12 is a risk variant for T1D subgroups with diagnosed age ≥ 12 years old

In our previous T1D GWAS, 14q12 regions were associated with T1D risk with the leading variant rs912304 (Additional file 1: Table S2). In a further case–control study, we first replicated and validated the association between rs912304 and T1D (P = 0.0177). Combined results with the GWAS scanning showed this variant was associated with T1D risk (OR (95% CI) = 1.18 (1.09–1.27), P = 1.29E − 05) (Table 1). We further analyzed T1D subjects stratified by the age of T1D diagnosis and found that this variant was a suggestive risk variant for those diagnosed age ≥ 12 years old (OR (95% CI) = 1.21 (1.12–1.31), P = 2.54E − 06), but not those diagnosed age < 12 (OR (95% CI) = 1.05 (0.94–1.17), P = 0.377). A significant heterogeneity was observed between these two subgroups (I² = 76.6%, P = 0.039). We also stratified T1D subjects according to the number of autoantibodies and found that rs912304 was associated with single (OR (95% CI) = 1.17 (1.07–1.27), P = 1.75E − 04) and multiple (OR (95% CI) = 1.14 (1.04 − 1.26), P = 1.03E − 02) antibody-positive T1D individuals without heterogeneity (I² = 0, P = 0.692) (Table 2). In addition, some T2D subjects might be misclassified as late-onset T1D [24]. T1D subjects diagnosed age ≥ 12 were stratified by autoantibody status, and we found that both the subjects with multiple autoantibodies and those with single autoantibodies showed significant association with the variant rs912304 without any heterogeneity (Additional file 1: Table S3).

Table 1 The association between rs912304 in STXBP6 and T1D risk

Table 2 The association between rs912304 in STXBP6 and T1D subgroups stratified by islet-specific autoantibody status and age at T1D diagnosis

Rs912304 affects residual islet function in T1D individuals of diagnosed age ≥ 12 years old

We next evaluated the impact of the variant rs912304 on islet function indexes in healthy individuals based on OGTT. This variant was significantly correlated with CIR and DI (P = 0.007 and 0.047, respectively) (Table 3). Subsequently, we assessed its effect on residual β-cell function in newly diagnosed T1D individuals, and it was correlated with fasting C-peptide level (P = 0.0289) (Fig. 1A) not the AUC of 2-h C-peptide level (Fig. 1D). After further stratified by the age of T1D diagnosis, this variant was significantly associated with fasting and AUC of 2-h C-peptide level (P = 0.0186, and 0.0298, respectively) (Fig. 1C, F) in T1D subjects with diagnosed age ≥ 12, which was consistent with their association with residual β-cell function (≥ 200 pmol/L) in this subgroup (P = 0.037) (Additional file 1: Table S4). However, no such correlation was found in those diagnosed age < 12 (Fig. 1B, E; Additional file 1: Table S4). In addition, we also investigated whether this variant affects the positive rate of islet-specific autoantibodies in T1D subjects, but no significant correlation was observed (Fig. 1G–I, Additional file 1: Table S5).

Table 3 The correlation between rs912304 and glycemic quantitative traits in 2353 healthy individuals based on OGTT

Fig. 1

The relationship between the variant rs912304 and residual islet function and islet autoimmunity in T1D individuals. A-C The association between rs912304 and fasting C-peptide level in total newly diagnosed T1D subjects, those diagnosed age <12 and ≥12 years old, respectively. D-F The association between rs912304 and 2-hour C-peptide AUC in total newly diagnosed T1D subjects, those diagnosed age <12 and ≥12 years old, respectively. G-I The association between rs912304 and the positive rate of ZnT8A, GADA, and IA-2A in total T1D subjects. The genotype distribution of rs912304 in newly diagnosed T1D subjects: GG = 52, TG = 102, TT = 48; T1D subjects with diagnosed age <12: GG = 10, TG = 26, TT = 8; and T1D subjects with diagnosed age ≥12: GG = 42, TG = 76, TT = 40. The genotype distribution of rs912304 in islet autoimmunity: GG=616, GT=1004, TT=420 for ZnT8A; GG=644, GT=1061, TT=443 for GADA; GG=639, GT=1059, TT=443 for IA-2A. The above correlations were analyzed by linear regression analysis and were corrected for sex and duration of T1D disease. Data are mean with 95%CI. *P<0.05 was considered as significant

Rs912304 is a functional variant regulating STXBP6 expression in islets

We evaluated the potential functional variant(s) in 14q12 risk regions. We first used RegulomeDB to prioritize functionally important variants within high linkage disequilibrium (LD, r² ≥ 0.8) to rs912304. Among the twelve high LD variants, six (rs912304 and the other five variants) were ranked 1d-1f with high posterior probability (scores > 0.5) (Additional file 1: Table S6). Next, we verified the results of functional predictions using HaploReg v4.2 [19], which utilized epigenomics data (chromatin states and histone H3 marks) primarily from Roadmap Epigenomics and ENCODE projects. In this dataset, we found that only rs912304 exhibited spatial overlap with enhancer, an active chromatin state, in tissues related to islet function (both pancreas and pancreatic islets) (Fig. 2A, Additional file 1: Table S7). It co-localized more specifically with an enhancer region marked by H3K4me1 broad peaks in the pancreas and pancreatic islets (Fig. 2A). Further functional follow-up in the published bulk islets dataset [25, 26], we also found that rs912304 mapped to ATAC-Seq peaks (Fig. 2B). Collectively, rs912304 is the most potential functional variant in pancreatic tissue and islets.

Fig. 2

Rs912304 is a functional variant regulating STXBP6 expression in islet β-cells. A The overlaps between rs912304 and chromatin state, histone H3 marks in Roadmap Epigenomics and ENCODE projects. The color red means rs912304 exhibited spatial overlaps with peaks of epigenomics data (chromatin states or histone H3 marks), blue means no spatial overlaps, and grey means the data were missing. B UCSC genome browser data showed that ATAC-Seq peaks overlapped with a genomic sequence containing rs912304 in published human islets. The location of rs912304 is highlighted with the red line. C Stxbp6 is the most significant eQTL gene regulated by rs912304 in 404 human islets derived from the Tiger database, highlighted in red. D The plasmid constructs containing different alleles of rs912304 enhancer and STXBP6 gene promoter. E-G The transcriptional activity of STXBP6 promoter after transfected to HEK-293T, islet β cell line INS-1E, and MIN6, respectively. The transcriptional activity was detected by the luciferase dual reporter gene method. Data are mean±SD. P-values were calculated by t-test, and at least three independent replicates were performed in each group. **P<0.01; ***P<0.001

As intron variants likely affect gene regulation [7], we further assessed the most putative candidate gene regulated by rs912304. Predicted by Open Targets Genetics [20], we found that among the six assigned candidate genes, STXBP6 was functionally implicated as the top-ranked potential causal gene with the highest variant to gene (V2G) score (Additional file 1: Table S8). Multi-tissue eQTL comparison from the GTEx v8 database revealed that STXBP6 is the unique eQTL gene regulated by rs912304, with a meta-analysis P-value = 1.32871e − 77 (Additional file 2: Fig. S1). Tiger, an integrative dataset with 404 human islets [22], also showed that STXBP6 is the most significant eQTL gene (P = 8.6e − 3) affected by rs912304 (Fig. 2C). Then, we performed luciferase reporter assays in HEK-293 T cells and two widely used cellular models of the β-cell (rat INS-1E and mouse MIN6). The construct containing the G allele of rs912304 exhibited higher enhancer activity of the STXBP6 gene than that of the T allele in HEK-293 T (P < 0.001), INS-1E (P < 0.01), and MIN6 (P < 0.01) (Fig. 2D–G). These data indicate that rs912304 was the most functional variant in pancreatic islets, which regulated the candidate gene STXBP6 as an enhancer.

T1D-related proinflammatory cytokine, not high glucose or FFA, up-regulates Stxbp6 expression on islet beta cells

To investigate the effects of different stress on Stxbp6 expression in islet β cells, we simulated them with three conditions in vitro: inflammatory cytokine, high glucose (20 mM glucose), and high fat (FFA). Under high glucose or FFA, Stxbp6 expression was not altered in cultured INS-1E and MIN6 cell lines (Fig. 3A, B). The results were similar to those observed under similar conditions derived from the published GEO database (Additional file 2: Fig. S2). However, we reanalyzed single-cell RNA sequencing (scRNA-seq) in the Human Pancreas Analysis Program (HPAP) [27]. The results showed that STXBP6 expression was upregulated in islet beta cells of T1D subjects (Fig. 3C). Meanwhile, after stimulating with T1D-related proinflammatory cytokine (TNF-α, IL-1β, and IFN-γ), we observed that Stxbp6 mRNA level was increased in primary mouse islets, and three β-cell lines (INS-1E, MIN6, and β-TC6) (all P < 0.05), (Fig. 3D–G). These were consistent with the fact that the Stxbp6 protein level was upregulated in the islet β-cell line stressed with cytokines (P < 0.05) (Fig. 3H, I). These data suggested that Stxbp6 was specifically stimulated in islet beta cells by T1D-related cytokines but not high glucose or lipids.

Fig. 3

Stxbp6 is significantly upregulated in pancreatic islet β-cells stimulated with T1D-related proinflammatory cytokines. A-B The alternation of Stxbp6 mRNA levels in INS-1E and MIN6 under high glucose (20mM glucose) and high lipid (0.5mM palmitic acid) stimulation. C The difference of STXBP6 mRNA expression in human islets from nPOD. AAB means islet antibody-positive normoglycemic population. D-G The mRNA level of Stxbp6 in primary mouse islets, INS-1E, MIN6, and β-TC6 cell lines after stimulation with mouse/rat cytokines (5 ng/ml IL-1β, 100 ng/ml INF-γ and 10 ng/ml TNF-α) for different periods. H-I Alternation of Stxbp6 protein level in INS-1E and MIN6 cell lines after 6 and 24 h stimulation with the above proinflammatory cytokines. Data are mean±SD. P-values were calculated by t-test, and at least three independent replicates were performed in each group. *P<0.05; **P<0.01; ***P<0.001

Stxbp6 promotes β-cell apoptosis and inhibits insulin secretion via regulating Gadd45β and Glut1 expression respectively

We first over-expressed Stxbp6 in INS-1E cells with recombinant adenovirus to assess its effect on apoptosis. Over-expression of Stxbp6 did not affect apoptosis-associated protein Caspase 3 but significantly increased the expression of cleaved Caspase 3 (Fig. 4A). The apoptosis detection by TUNEL assay also revealed a significant increase in the Stxbp6 overexpression group (Fig. 4B) in INS-1E cells, suggesting that the Stxbp6 has a role in promoting apoptosis in pancreatic β-cells. Furthermore, although no significant change in the transcript levels of Ins1 and Ins2 was observed in the Stxbp6 overexpressed INS-1E cells (Fig. 4C), insulin secretion, not intracellular insulin content, was suppressed under high glucose stimulation conditions (P < 0.05) (Fig. 4D, E).

Fig. 4

Effect of Stxbp6 overexpression on apoptosis, insulin synthesis, and secretion in INS-1E cells. A Protein expression of Stxbp6, Caspase 3 and cleaved- Caspase 3 in INS-1E cells treated with Stxbp6 recombinant adenovirus. The relative change of protein was calculated on the right. B Apoptosis detection by TUNEL assay. Scale bar 100μm. C Changes in insulinogenic gene Ins1/Ins2 mRNA levels after overexpression of Stxbp6 in INS-1E cell. D Insulin secretion at glucose concentrations of 2 mM and 20 mM, respectively; E, Intracellular insulin levels after overexpression of Stxbp6. Data are mean±SD. P-values were calculated by t-test, and at least three independent replicates were performed in each group. *P<0.05; **P<0.01; ***P<0.001

Further RNA sequencing revealed that 67 and 285 genes were significantly down-regulated or up-regulated in Stxbp6 overexpressed INS-1E cells (RPKM > 1) (Fig. 5A, B). Pathway enrichment analysis revealed that the downstream genes regulated by the Stxbp6 gene were significantly enriched in insulin secretion and the p53 signaling pathway (Fig. 5C, Additional file 1: Table S9). We validated the RNA-seq results and showed that about half of the downstream genes were associated with cell survival and apoptosis (e.g., Cdk2, Gadd45β, and Tp73; Fig. 5D), insulin secretion (e.g., Eno2 and Slc2a1; Fig. 5E), and other related genes (e.g., Cxcl12, Cygb, and Mgarp; Additional file 2: Fig. S3A, B). Knockdown of Stxbp6 with SiRNA (Si-332) in cytokine-stimulated MIN6 cells (Additional file 2: Fig. S3C), we found that mRNA levels of genes related to apoptosis (Cdk2 and Gadd45β, Fig. 5F), endocrine system (Slc2a1, Fig. 5G), immune system, and other signaling pathways (e.g., Cxcl12, Additional file 2: Fig. S3D, E). We further verified the apoptosis- and islet-function-related genes, consistent in both overexpression and knockdown experiments at the protein level. We found that the expression level of the Gadd45β was significantly increased in INS-1E cells after the overexpression of the Stxbp6 (P < 0.01), suggesting that Gadd45β was associated with Stxbp6 induced apoptosis in INS-1E cells. In addition, while Eno2 did not change significantly, the expression of Glut1 decreased significantly (P < 0.01), suggesting that the Glut1 was associated with Stxbp6-mediated reduction of insulin secretion in INS-1E cells (Fig. 5H).

Fig. 5

Gadd45β and Glut1 are two key downstream genes regulated by Stxbp6 in islet β-cell lines. A Volcano plot of recombinant adenovirus overexpressing Stxbp6 gene in INS-1E cells with downstream differential gene expression. Red dots represent genes with rising differential expression, and blue dots represent genes with falling differential expression. DEGs were screened using logFC>0.5, FPKM>1, P-adj<0.05, Q value<0.1. B Heatmap of Stxbp6 gene downstream related genes based on RNA-Seq hints. C Analysis of pathway enrichment of DEGs using the KEGG database. D Effect of overexpression of Stxbp6 gene in INS-1E cells on cell proliferation and apoptosis genes. E Effect of overexpression of Stxbp6 gene in INS-1E cells on genes related to endocrine system. F-G Effect of knocking down Stxbp6 gene in MIN6 cells on genes related to the endocrine system and apoptosis. H Western blotting detection of protein expression of Gadd45β, Glut1, and Eno2 in INS-1E cells treated with Stxbp6 recombinant adenovirus. The relative change of protein was calculated on the right. Data are mean±SD. P-values were calculated by t-test, and at least three independent replicates were performed in each group. *P<0.05; **P<0.01; ****P<0.0001

Discussion

In this study, we first reveal that rs912304 in STXBP6 is a T1D risk variant associated with fasting and 2-h C-peptide levels in newly diagnosed T1D individuals diagnosed age over 12 years old. Next, rs912304, as an enhancer, functionally regulates STXBP6 expression in islet beta cells. Finally, Stxbp6 can inhibit insulin secretion and promote islet β-cell apoptosis by regulating Glut1 and Gadd45β expression.

Studies have observed an increasing incidence of adult T1D in multiple ethnicities [28] and revealed significant heterogeneities in T1D subgroups of age at T1D diagnosis [29]. Although LD score regression analysis showed a highly correlated genetic basis in early- and late-onset T1D individuals [13], our previous study demonstrated that HLA regions showed stronger effect sizes in early-onset T1D individuals [13]. Studies also identified unique risk regions and the same risk regions with more effect size in early-onset T1D individuals [30, 31]. Steck et al. [32] further found that non-HLA risk variants, such as rs6476839 in GLIS3, and rs4900384 near VRK1 were significantly associated with progression to T1D in individuals < 12 years old, whereas rs3184504 in SH2B3 was significant in those ≥ 12. This study revealed that rs912304 in STXBP6 is a suggestive non-HLA risk variant for late-onset T1D individuals (diagnosed age ≥ 12). These suggested further exploration with GWAS to hunt for more potential non-HLA risk regions for T1D subgroups with different diagnosed ages.

Prior studies demonstrated a shared genetic basis for T1D and islet-specific autoantibody positivity, such as rs1990760 in IFIH1 for IA-2A and rs12722495 in IL2RA for GADA positive rate [33, 34]. However, we found no significant association between rs912304 and the positivity of any islet-specific autoantibodies (ZnT8A, GADA, or IA-2A), indicating the variant did not affect islet autoimmunity. Multiple T1D genetic risk variants also affected islet function, for instance, rs3825932 in CTSH affected insulin secretion [35]. For the variant rs912304, prior studies have revealed its effect on 2-h glucose challenge [36] and HbA1C levels [37]. Our study found that rs912304 was associated with glucose-stimulated insulin secretion in healthy individuals based on OGTT, and further affected residual islet function in newly diagnosed T1D individuals, especially of diagnosed age ≥ 12. Thus, we speculated that rs912304 might contribute to T1D by regulating islet function.

For complex common diseases, multi-risk variants reaching GWAS significance are usually observed in a certain risk region, and the lead variant is not always the most likely one to be functional. Moreover, most risk regions contain multiple risk variants and candidate genes, making it crucial to identify the most functional variant and causal gene within these regions. Prior studies [7, 38, 39] have discovered the most functional variants and putative implicated genes for T1D risk using fine mapping, epigenetic data, etc. Studies also found that 14q12 was the risk region for white blood cell and monocyte count [40], rs2332462 and rs2038700 in STXBP6 (in high LD to rs912304) were predicted to be the most causal variants for these two counts, respectively. However, our study indicated that rs912304 is the most casual variant in pancreatic tissues or islets by overlapping with DNase I hotspots, H3K4me1, and ATAC-Seq. These implicated that different tissues might have different functional variants. In addition, from the ChIP-Seq experiments integrating into HaploReg v4.2 [19], DNA regions spanning rs912304 were found to bind two small MAF proteins (MAFF and MAFK), which regulate islet beta cell survival [41]. Regulatory motifs prediction by Position Weight Matrix (PWM) showed a different allele-specific alternation of this variant to pdx1, which involves the early development of the pancreas and regulates glucose-dependent insulin gene expression [42]. Collectively, these data further demonstrated rs912304 as an islet-specific functional variant.

A prior study revealed that STXBP6 induced by needle-shaped hydroxyapatite (n-nHA)-treated macrophages triggers autophagy, which markedly promotes macrophage fusion into multinucleated giant cells (MNGCs) [43], which implicated STXBP6 as an immune function related gene. Moreover, 14q12 was reported to be a T2D risk region with the leading variant rs11159347 [17]. Using locus-to-gene pipeline (L2G) from Open Targets Genetics [20] prioritized STXBP6 as the most potential causal gene for T2D with the highest L2G score. This study identified STXBP6 as the most putative causal gene in islet beta cells, regulated by the functional T1D risk variant rs912304, independent of rs11159347. Prior studies have demonstrated that the proinflammatory cytokines IFN-γ, IL-1β, and TNF-α are produced mainly by islet-infiltrating immune cells and affect beta cell function and apoptosis, contributing to T1D development [44, 45]. Our study further revealed that STXBP6 was upregulated in islet β-cells under such proinflammatory cytokine stimulation. Consistently, human islets scRNA-seq analysis in the HPAP dataset [27] showed that STXBP6 was increased in islet beta cells of T1D subjects.

Despite its relatively low expression in islet beta cells, STXBP6 could participate in the expansion of fusion pores during the insulin exocytosis process [15, 16] and inhibit glucose-stimulated insulin secretion regulated by SOX4 [15, 16]. Although Glut2 is well recognized as the glucose transporter contributing to insulin secretion [46], this gene was not regulated by Stxbp6 in mouse pancreatic beta cells. However, this study indicated that Glut1, another glucose transporter isoform, was down-regulated by Stxbp6 in islet beta cells, facilitating glucose-stimulated insulin secretion [47]. Moreover, our study suggested that Stxbp6 also promoted islet β cell apoptosis, consistent with the downstream genes of Stxbp6 that were enriched in the p53 pathway. Gadd45β, a typical gene in the p53 pathway, was upregulated by Stxbp6 in islet beta cells and is involved in the positive regulation of the apoptotic process [48, 49]. Interaction networks from GeneMANIA (http://genemania.org/) and STRING (https://cn.string-db.org/) revealed that Stxbp6 might regulate Glut1 and Gadd45β indirectly. As cumulative evidence supported a critical role for β-cells contributing to T1D pathogenesis [11], we speculated that STXBP6 might contribute to T1D pathogenesis via regulating islet function and beta cell apoptosis.

Our study also has several limitations. First, due to the differences in minor allele frequency, LD patterns, and effect size in different populations, no single population is sufficient for fully uncovering the variants underlying T1D in all populations. Thus, studies should further explore the association between rs912304 and late-onset T1D risk in other populations. Second, although Stxbp6 regulated Glut1 and Gadd45β in islet beta cells, the exact contributions of these two genes in mediating the effects of Stxbp6 on beta cell function and survival still need further investigation. Finally, the function of Stxbp6 in the immune regulation related to T1D risk should be further studied, as both innate and adaptive immune abnormalities play critical roles in T1D pathogenesis [50].

Conclusions

In summary, we have identified rs912304 as a functional variant for late-onset T1D risk, affecting residual islet function in newly diagnosed T1D individuals. This variant regulates STXBP6 expression in islet beta cells, a causal gene promoting islet β-cell apoptosis and inhibiting insulin secretion. These findings expanded the genomic landscape regarding T1D risk and prompted the genetic heterogeneity in T1D subgroups with different diagnosed ages. This study also showed the evidence supported a role for β-cells in T1D pathogenesis, which needs further exploration.