Genetic Variants Associated with Neuropeptide Y Autoantibody Levels in Newly Diagnosed Individuals with Type 1 Diabetes

(1) Autoantibodies to the leucine variant of neuropeptide Y (NPY-LA) have been found in individuals with type 1 diabetes (T1D). We investigated the association between the levels of NPY-LA and single nucleotide polymorphisms (SNP) to better understand the genetic regulatory mechanisms of autoimmunity in T1D and the functional impacts of increased NPY-LA levels. (2) NPY-LA measurements from serum and SNP genotyping were done on 560 newly diagnosed individuals with T1D. SNP imputation with the 1000 Genomes reference panel was followed by an association analysis between the SNPs and measured NPY-LA levels. Additionally, functional enrichment and pathway analyses were done. (3) Three loci (DGKH, DCAF5, and LINC02261) were associated with NPY-LA levels (p-value < 1.5 × 10−6), which indicates an association with neurologic and vascular disorders. SNPs associated with variations in expression levels were found in six genes (including DCAF5). The pathway analysis showed that NPY-LA was associated with changes in gene transcription, protein modification, immunological functions, and the MAPK pathway. (4) Conclusively, we found NPY-LA to be significantly associated with three loci (DGKH, DCAF5, and LINC02261), and based on our findings we hypothesize that the presence of NPY-LA is associated with the regulation of the immune system and possibly neurologic and vascular disorders.


Introduction
Type 1 diabetes (T1D) is an autoimmune disease caused by the destruction of pancreatic islet β cells [1]. The immune-mediated attack on β cells is caused by complex interactions between environmental factors and multiple risk-conferring genes [2]. More than 75 susceptibility loci for T1D have been identified with genome-wide significance [3]. The human leukocyte antigen (HLA) region on chromosome 6p21 contributes to approximately 50% of the familial genetic risk [4], and the HLA region confers the strongest risk for islet autoantibodies which are biomarkers of T1D autoimmunity [2,5,6]. HLA-DR3-DQ2 is associated with autoantibodies to glutamic acid decarboxylase (GADA) as the first appearing islet autoantibody [7]. Additionally, the presence of autoantibodies to islet antigen-2 (IA-2A) is associated with HLA-DR4-DQ8 [2,5,6,8] but is negatively associated with HLA-DQ2 [9].
Participant demographics of all the 560 individuals with T1D are shown in Table 1. This study, including the existing protocols from DanDiabKids, Steno 2007, the CPH remission phase 2015 cohort, and the Danish Remission Phase Study, has been carried out per the Declaration of Helsinki, and all were approved by the Danish Ethical Committee and by The Danish Data Protection Agency. Autoantibody measurements: Serum samples from the 560 individuals, taken within 3 months from clinical onset of T1D, have been stored at −80 • C and were analyzed for NPY-LA using a radioligand binding assay (RBA), described elsewhere [22,26,30]. In brief, NPY-L was subjected to coupled in vitro transcription-translation to generate 35 S-methionine radiolabeled autoantigens incorporated into an overnight immunoprecipitation/filtration assay. A polyclonal rabbit IgG antibody to NPY (Abcam, Cambridge Science Park, Cambridge, UK), in seven-doubling dilution steps, was used for converting Protein A Sepharose-bound radioactivity into in-house units. Analyses for NPY-LA were done in duplicates. The mean value of the intra-assay coefficient of variation (CV) for duplicates was 13%. The NPY-LA assay has not been part of the international standardization workshops [31,32].
GADA and IA-2A were analyzed within 3 months from the clinical onset of T1D and were analyzed in a standard RBA, previously described in detail elsewhere [28,32]. Levels were used instead of dichotomous values since GADA quartiles above Q1 show increasing NPY-LA values [26]. Additionally, individuals with IA-2A in the upper quartile compared to the lower quartile have a fourfold risk of developing T1D [6]. Titration above 500 U/mL for GADA and 250 U/mL for IA-2A have not been performed [28]. However, undiluted values above the standard curves were used [6,31,33]. Levels above the standard curve are underestimated as end-point titrations tend to attain higher levels.
SNP genotyping of the three sub-cohorts: A total of 560 individuals had genotyping data available, of which 486 were previously genotyped on the ImmunoChip array described elsewhere [12], and 74 were genotyped on the Illumina Infinium Global Screening Array (GSA version 24.v2.0). Before quality control (QC), the ImmunoChip array dataset and GSA array dataset had 177,022 and 730,059 SNPs, respectively. After QC and frequency pruning (-geno 0.1 -maf 0.05 -hwe 0.001 -mind 0.05) using PLINK v1.90, the ImmunoChip and GSA datasets had a genotyping rate of 0.99 and had 109,021 SNPs and 306,455 SNPs, respectively.
The two datasets were merged using PLINK v1.90 before proceeding with the imputation. A total of 13,066 SNP markers overlapped between the two platforms. Before merging, we made sure that both datasets had the same strand and build concordance. After QC and strand alignment, 402,410 SNPs remained in the merged dataset. After removing duplicate variants, 397,929 SNPs remained for imputation analysis.
Imputation with 1000 genomes: To explore regions not well covered by the Im-munoChip, imputation was performed on the merged dataset using the 1000 Genomes haplotype reference panel. Chromosome strand and allele order in the merged dataset were adjusted to match the reference 1000 Genome Phase 3 dataset (excluding non-European samples) using the conform-gt program. Pre-phasing was performed using SHAPEIT v2 [34] (with European population (-effective-size = 11,418)). Imputation was then performed using IMPUTE2 v2.3.2 [35] and 1000 Genome phase3 panel [36] with default parameters. For imputation to run efficiently, each chromosome was split into 5 MB chunks for analysis, as the IMPUTE2 program produces higher accuracy over short genomic regions. The IMPUTE2 results were obtained in a gene/sample output format. The chromosome chunks were concatenated for each chromosome and converted to VCF format using BCFtools v1.8. The annotation of the VCF file was done with the annotate function in BCFtools using the dbsnp version150 (build GRCh37p13). PLINK was used to filter insertion-deletion mutations (INDELS), all variants with 2+ alternate alleles, and SNPs with no-call alleles. The files were then converted to PLINK bed format. To assess imputation quality, the IMPUTE2 internal quality metrics (INFO score) were obtained. The INFO values range from 0 to 1, where a higher value indicates the increased quality of an imputed SNP. A total of 3,119,616 SNPs remained after post imputation QC and after filtering low imputation score SNPs (INFO < 0.4).
Association and pathway analysis: Linear regression was used to evaluate the association between the imputed genetic variants and NPY-LA. Out of the 560 samples, 11 Non-European samples were removed before proceeding with the association analysis. A quantitative association test, using a linear function in PLINK v 1.90, was performed on the NPY-LA measurements in 546 out of 549 samples (3 samples had missing values) and the imputed SNPs (n = 3,119,616). The linear model was adjusted for the following covariates: sex, age at clinical onset of T1D, ethnicity, log2 of GADA, and IA-2A measurements. Ethnicity was coded as a factor with two levels (1 = Danish, 2 = European/Turkish). The age at clinical onset of T1D was grouped into four levels: below 5, 5-10 or 10-15 years, and above 15 years. The value of genomic inflation-estimate (lambda), a measure of inflated p-values (based on median chisq), was 1. Lambda statistic should be close to 1 if the points fall within the expected range or greater than 1 if the observed p-values are more significant than expected. Manhattan plot and QQ-plots before and after adjustment for genomic control were created using the qqman package in R version 4.1.0. A p-value cutoff of p < 1.0 × 10 −4 was used to identify the top signals from the association analysis.
Variant Effect Predictor (VEP) in Genome Reference Consortium Human Build 38 Ensembl release 104 (GRCh38) was used to annotate the top variants. Expression quantitative trait loci (eQTLs) in whole blood were extracted for the top variants using the GTEx portal (The Genotype-Tissue Expression Project), hereby retrieving p-values and normalized effect sizes (NES). The NES of eQTLs is defined as the slope of the linear regression. It is computed as the effect of the alternative allele relative to the reference allele in the human genome reference GRCh38/hg38. A p-value cutoff of p < 0.05 was used as the threshold for the eQTL analysis.
Pathway analysis and network establishment were performed on the annotated top variants using Ingenuity Pathway Analysis (QIAGEN IPA, RRID: SCR_008653).

Results
Of the 560 individuals with T1D, 0.5% had missing NPY-LA values (3/560). For GADA, 53% and 21% were positive and negative, respectively. IA-2A was detected in 51% of individuals, while 23% were negative. The measurements were missing in 26% for both GADA and IA-2A, and the missing values of the two variables were set to zero. The ethnic distribution of the genotyped samples is shown in Table S1.
The association analysis between the NPY-LA measurements and the genotyped score of each SNP is illustrated in a QQ-plot ( Figure S1) that demonstrates observed and expected −log10 p-values. A Manhattan plot illustrates chromosomes with significant SNPs (Figure 1).
Genes 2022, 13,869 Of the 560 individuals with T1D, 0.5% had missing NPY-LA values (3/5 GADA, 53% and 21% were positive and negative, respectively. IA-2A was detecte of individuals, while 23% were negative. The measurements were missing in 26% GADA and IA-2A, and the missing values of the two variables were set to zero. Th distribution of the genotyped samples is shown in Table S1. The association analysis between the NPY-LA measurements and the ge score of each SNP is illustrated in a QQ-plot ( Figure S1) that demonstrates obser expected −log10 p-values. A Manhattan plot illustrates chromosomes with sig SNPs (Figure 1).
Expression quantitative trait loci (eQTLs) for the selected 348 SNPs in whole blood were retrieved from GTEx, and the analysis identified several eQTLs of DCAF5 (Table 3). Only 24 out of the 348 selected variants had eQTL effects on the host gene, and the majority of these (19 variants) mapped to the DCAF5 locus. Other loci with eQTL variants included ABCC2, CCDC6, PABPC1L, DENND2C, and TTF2 (Table 3, p-value < 0.05).  We identified a set of top molecules based on pathway analysis of the 348 SNPs (using IPA), mapping to 52 genes (Table S3). We highlight the presence of SNP rs79095830 from host gene ABCC2 and SNP rs3746583 from gene DCAF5, of which eQTL was identified ( Table 3 and Table S3).
Finally, the IPA analysis generated two networks that illustrate the influence of the NPY-LA associated genes on the p38 mitogen-activated protein kinases (MAPK) pathway, where one network includes JNK, Akt, and PI3K ( Figure 2a) and the other includes JNK (Figure 2b). We identified a set of top molecules based on pathway analysis of the 348 SNPs (using IPA), mapping to 52 genes (Table S3). We highlight the presence of SNP rs79095830 from host gene ABCC2 and SNP rs3746583 from gene DCAF5, of which eQTL was identified (Tables 3 and S3).
Finally, the IPA analysis generated two networks that illustrate the influence of the NPY-LA associated genes on the p38 mitogen-activated protein kinases (MAPK) pathway, where one network includes JNK, Akt, and PI3K ( Figure 2a) and the other includes JNK (Figure 2b).

Discussion
This observational study of 569 children with T1D was designed to examine genetic loci associated with levels of NPY-L autoantibodies at the clinical onset of T1D. Association analysis identified three loci DGKH, DCAF5, and LINC02261, which had SNPs

Discussion
This observational study of 569 children with T1D was designed to examine genetic loci associated with levels of NPY-L autoantibodies at the clinical onset of T1D. Association analysis identified three loci DGKH, DCAF5, and LINC02261, which had SNPs significantly associated with NPY-LA levels. When including GADA and/or IA-2A levels in the model, the associations did not change, indicating that the three loci are not important for the known association between NPY-LA and the islet autoantibodies [26]. The functional enrichment analysis (eQTL) identified SNPs associated with variations in expression levels in whole blood of DCAF5, ABCC2, CCDC6, PABPC1L, DENND2C, and TTF2, suggesting potential interactions between these genetic factors and NPY-LA levels. The top molecules of the pathway analysis included DGKH, DCAF5, PABPC1L, and ABCC2.
The gene, DGKH, encodes a member of the diacylglycerol kinase enzyme family [37], while DCAF5 encodes DDB1 and CUL4-associated factor 5 [38], and LINC02261 is a nonprotein-coding lncRNA. None of the genes are known susceptibility loci for T1D or islet autoimmunity. However, DCAF5 has previously been linked with T2D [39], and NPY-LA has been found in individuals with T2D [22], suggesting NPY-LA influences diabetic processes shared between T1D and T2D. ABCC2 is part of the adenosine triphosphatebinding cassette family (ABC transporters) located in the plasma membrane that enables an active transport of various compounds across the cell membrane [40,41]. Both the pathway top molecules and the pathway networks indicate that NPY-LA is associated with changes in the mitogen-activated protein kinase (MAPK)-signaling pathway, gene transcription, protein modification, and immunological functions.
MAPK-signaling pathway: Several identified genes affect cell transcription and apoptosis signaling. DGKH is known to regulate intracellular concentrations of phosphatidic acid and diacylglycerol and has a key role in promoting cell growth by activating the MAPK-signaling pathway [37]. MAPKs are implicated in altered transcriptions of genes and, thus, numerous cellular behaviors such as survival, proliferation, differentiation, and apoptosis [42,43]. The network identifies DGKH as a hub node directly or indirectly activating several kinase signaling pathways, e.g., Akt, JNK, PI3K, and LYN. Akt affects glucose uptake and utilization, glycogen, fatty acid, and protein syntheses, and has been implicated in several diseases, including diabetes, autoimmunity, and insulin resistance [44,45]. The c-Jun N-terminal kinase (JNK) pathway is one of the major signaling cassettes of the MAPK signaling pathway and plays an important role in apoptosis, inflammation, cytokine production, and metabolism [46]. The Phosphoinositide 3-kinase (PI3K) complex is linked to an extraordinarily diverse group of cellular functions but has an important role in glucose regulation and may be involved in the clinical onset of T1D [47]. LYN is a member of the Src family of tyrosine kinases and plays an important role in regulating innate and adaptive immune responses, responses to growth factors, and cytokines [48]. Interestingly, LYN was recently identified as a novel T1D susceptibility locus [3]. Also, DCAF5 indirectly interacts with Akt. Other genes from the pathway analysis and eQTLs support an association with the MAPK-signaling pathway (CCDC6, PRKG1, TTF2, PABPC1L, SND1, and ABCC2). This suggests that NPY may impact the T1D phenotype by complex regulation of pro-and anti-apoptotic signaling pathways.
Furthermore, the MAPK pathway is regulated by enzymatic post-translational modification of ubiquitination [49,50], a proteasomal degradation to regulate the activity of targeted proteins [41]. DCAF5 is assumed to function as a substrate receptor for the CUL4-DDB1 E3 ubiquitin-protein ligase complex helping ubiquitination [38]. Ubiquitin chains have been shown to play important roles in inflammatory signaling pathways [51], and interestingly, both DGKH [38] and DCAF5 influence neutrophiles of the innate immune system. Studies have shown that increasing NPY levels in mice lead to hyperglycemia, impaired glucose tolerance, hyperphagia, and obesity. On the contrary, NPY knockdown in rats leads to the development of brown adipose tissue and elevated thermogenic activity [21,[52][53][54][55]. Additionally, injection of NPY autoantibodies in mice decreases food consumption [56].
The findings suggest that NPY-LA plays a role in post-translational modification and regulation of the MAPK signaling pathway via ubiquitination, and through that influence survival of β cells or activation of inflammatory cells, thus influencing the autoimmune response in T1D and glucose regulation. Interestingly, NPY acts through the Y1 receptor to increase neuronal proliferation, mediated through activation of the MAPK pathway [57].
Neurological disorders: DCAF5 has been associated with neurogenic bowel disease [58], which is interesting since 20% of individuals with T1D in time will experience neurogenic bowel disease due to irreversible autonomic neuropathy [59]. Furthermore, prior studies have shown that defects in the expression level or stability of ABC transporters result in neuropathies [41]. Moreover, individuals with long-duration T1D are complicated with autonomic neuropathy when they have markedly reduced or absent responses of NPY to insulin-induced hypoglycemia [60]. Since we assume that NPY levels are diminished in the presence of NPY-LA [61], we hypothesize there is a possible association between NPY-LA levels, neuropathy, gene variants of DCAF5, and the function of ABCC2. Although, a previous study including few participants with T1D and T2D found no relation between peripheral or cardiac autonomic neuropathy and NPY autoantibodies [22]. Studies of NPY-LA with more participants having long-duration T1D and T2D and including identifications of genetic variants are needed to further elaborate associations between DCAF5, ABCC2, NPY-LA, and neuropathies.
The present study identified a significant association between NPY-LA and DGKH, a gene that has previously been reported in bipolar disorders [62]. Also, the lncRNA LINC02261 has previously been found in schizophrenia and autism spectrum disorder [63,64]. Interestingly, several genes that were previously found to affect the MAPK-signaling pathway are known to be associated with schizophrenia, schizoaffective, and bipolar disorders [65]. Moreover, prior studies have shown that numerous lncRNAs and their target genes modulate signaling pathways known to play a role in T1D-related neuronal dysfunctions [66].
Additionally, our pathway analysis identified genes affecting calcium channels (PRKG1 [67]). NPY acts to modulate voltage-gated Ca 2+ -channels in the neuronal system [68], and, interestingly, inhibition of glucose-stimulated insulin secretion by NPY is assumed to be dependent on Ca 2+ in β cells [24]. Voltage-gated calcium channels are present in the membrane of neurons, β cells, and muscle cells (skeletal, cardiac, and smooth). They regulate intracellular processes such as neurotransmission, secretion, and contraction. Besides, Ca 2+ is involved in the vesicular trafficking of ABC transporters to the membrane [41].
In T1D, persistent hyperglycemia has previously been found to promote neurovascular dysfunction, leading to endothelial cell dysfunction, oxidative stress, increased neuronal cell apoptosis, and inflammation [59,66]. The association with endothelial dysfunction is supported by the MAPK-signaling pathway, an important regulator of inflammation on endothelial cells and the progression of atherosclerosis [69]. Interestingly, both mutations in the ABC transporter and the minor allele of NPY (NPY-P) have been associated with atherosclerosis [70,71]. Future genetic analyses of NPY-LA associations in individuals with T2D, who have neuropathy and vascular disorders, may elaborate on possible associations.
The strength of the present study involves the large number of samples at the clinical onset of T1D with genetic information to identify loci associated with NPY-LA levels in individuals with T1D. Though the lack of replication sample/data is a limitation, as well as the use of whole blood for functional annotation in the absence of follow-up experiments in relevant cells/models, an important strength is the combination of analyses identifying associations with SNPs, expressions, and pathway analyses. Finally, there may be limitations to using publicly available databases as study design and demographic data vary.

Conclusions
In conclusion, NPY-LA is significantly associated with three loci: DGKH, DCAF5, and LINC02261. It is hypothesized that the presence of NPY-LA is associated with activation of the innate immune system, where NPY-LA may interact with the MAPK pathway in the regulation of glucose homeostasis, inflammation, cell survival, ABC transporter, and Ca 2+ levels in β cells, leukocytes, endothelial cells, and neurons in the development of T1D and possibly in T2D neuropathies or vascular disorders.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.