Phenome-Wide Association Study of Latent Autoimmune Diabetes from a Southern Mexican Population Implicates rs7305229 with Plasmatic AntiGlutamic Acid Decarboxylase Autoantibody (GADA) Levels

1. Introduction

World Health Organization estimates that 5 - 14% of patients with diabetes have latent autoimmune diabetes in adults (LADA) [1,2]. It is described that subjects with LADA have intermediate features between type 1 diabetes (T1D) and type 2 diabetes (T2D). LADA is characterized by slower autoimmune destruction of the β-cells compared to T1D, where glutamate decarboxylase autoantibodies (GADA) are the predominant serum marker [3,4]. This autoimmune process results in progressive hyperglycemia, and thus, patients with LADA do not require insulin therapy at the onset of diabetes[5]. On the other hand, the etiology of LADA involves insulin resistance. As with T2D, the risk of LADA is modified by lifestyle factors such as body mass index and physical activity [6]. Moreover, the initial clinical presentation of LADA resembles that of T2D. Consequently, patients with LADA are initially diagnosed with T2D, and they receive inadequate therapies for their diabetes until insulin dependence develops [5,7].

Buzzetti et al. (2020) suggest that the phenotype of LADA is heterogeneous with some patients similar to T1D and others indistinguishable from T2D [8]. There are variations in the phenotype of LADA between populations, such as populations from Europe, the United Arab Emirates, and Yemen have intermediate features between T1D and T2D [9,10,11]. Nonetheless, in Chinese populations, LADA has similar insulin resistance compared to T2D, and less serum GADA levels [12,13]. The variations in the phenotype of LADA suggest that differences between populations could be attributed to genetic risk factors. In this sense, LADA has genetic traits that interact with the environmental factors of T1D and T2D [6]. The genetic risk factors together with non-genetic factors of LADA could influence the rate of β-cell destruction and thus lead to a heterogenous clinical presentation [14]. Our knowledge is that only one GWAS of LADA has been conducted in a European population; in this study, they found 5 polymorphisms in LADA (genes HLA-DQB1, PTPN22, INS, and SH2B3), these SNPs are related to T1D [15]. In the Chinese Han population, Song et al (2024) associated LADA and T1D with an RPS26 variant in a case-control study [16]. On the other hand, Mendelian randomization studies in the European population associated lifestyle risk factors in LADA with risk genotypes related to T2D (FTO and TCF7L2) [17,18]. In this sense, more studies are necessary to comprehend the genetics of LADA better. However, information about risk genotypes of LADA in non-European populations, including Latin American populations, is scarce. Genetic studies in the Latin American population could reveal LADA risk factors unique to this population. Therefore, this study aims to perform a GWAS with phenome-wide association study (PheWAS) of LADA of a population of patients with LADA from a hospital in southeastern Mexico and a biobank of the Mexican population.

2. Results

A total of 51 phenotypes were evaluated in our PheWAS analysis, which identified 981 SNPs associated with these phenotypes (p < 5x10^-5). The associations are provided in Supplementary Table 1.

In the case-control GWAS analysis for LADA and population-based controls, no genome-wide significant association was found. Still, we observed 8 single-nucleotide polymorphisms associated at the nominal level (p-value < 1x10^-5) (Supplementary Figure 1). Five SNPs are intronic in the LIFR and HDAC9 genes, two are intergenic variables, and one is located downstream of the TMEM20 gene (Supplementary Table 2).

In LADA cases (n=50), we associated GADA levels and one polymorphism with genome-wide significance (p < 5x10^-8). Likewise, we found 25 more variants associated with GADA levels at nominal significance (Supplementary Table 1). Table 1 presents significant polymorphisms associated with phenotypic features in patients with LADA (diabetes onset, height, and hypertension onset).

The SNP is rs7305229, a downstream variant of the FAIM2 gene (p = 1.84x10^-8, β = -1.27, SE = 0.18 and MAF = 0.22). Figure 1 includes the Manhattan plot for the GADA genome-wide associations. After, we performed a fine mapping of rs7305229 locus, including all SNPs within 1Mb. Figure 2 shows that the polymorphism rs7305229 had a 94.8% posterior inclusion probability of being causal for the GADA levels phenotype, the greatest posterior probability.

We realized a PheWAS with the GWAS Atlas database, which found 41 traits of the metabolic domain associated with rs7305229 (p < 5x10^-8). Bioelectrical impedance-related measures and body mass index were frequently associated with rs7305229. Other traits related to this SNP were salt consumption, age at menarche, and bone mineral density. Notably, we did not find an association between rs7305229 and the immunological domain (p > 1x10^-5) (Supplementary Table 3).

The rs7305229 was associated with GADA levels in our population and with BMI in the GWAS Atlas database. Thus, we performed a conditional analysis to identify possible dependent associations. The linear model between the normalized values of both traits showed no association (R²=0.026, p=0.262). Besides, the second genome-wide analysis had changes in comparison to the first analysis (p = 6.97x10^-8, β = -1.26, SE = 0.19).

3. Discussion

We explored GWAS associations in multiple phenotypes in individuals diagnosed with LADA from Southern Mexico, we showed that rs7305229 is associated with GADA antibody levels.

In the comparison between LADA cases and population-based controls, we identified associations with the LIFR (leukemia inhibitory factor receptor) gene. These associated SNPs of LIFR are intronic variants, and they probably change expression levels rather than alter their protein structure. In addition, LIFR is a transmembrane receptor that binds to IL-6 family members, including leukemia inhibitory factor and oncostatin-M [19]. LIFR is involved in lipolysis, adipose inflammation, hepatic triacylglycerol accumulation, and insulin resistance in mouse models [20,21]. Literature remarks that 56% of patients with T2D have NAFLD [22]. In the context of the changes in LIFR expression could be related to susceptibility to non-alcoholic fatty liver disease (NAFLD) and insulin resistance in patients with LADA, and the LIFR/JAK/STAT3 inflammatory pathway enhances adipocyte lipolysis [21]. Although, the evidence on subjects with LADA is scarce; a study of hospitalized patients in China reported that 24.6% of individuals with LADA had NAFLD [23]. In our population, we found fewer patients with LADA and NAFLD than with T2D and NAFLD; nonetheless, the risk for NAFLD in patients with LADA could be greater than in healthy subjects [24].

We associated a total of 981 SNP with nominal associations in different phenotypes, the phenotypes with more associations were age of hypertension onset (n=244), BMI (n=88), insulin doses (n=74), and cholesterol levels (n=67). Notably, these phenotypes resemble the metabolic syndrome criteria [25]. In addition, the LADA genetic constitution has polymorphisms associated with metabolic syndrome [14,15].

Our PheWAS identified a GWAS significant association with a 3’UTR genetic variant (rs7305229) located downstream of the FAIM2 gene (Fas apoptotic inhibitory molecule 2) with plasmatic GADA antibody levels in patients diagnosed with LADA. Polymorphisms of FAIM2 are associated with childhood obesity, autoimmune thyroid disease, and obsessive-compulsive disorder [26,27,28]. Besides, FAIM2 expression is associated with lymphocyte CD8+ antitumor activity. The FAIM2 protein protects from cell death induced by Fas and plays an important role in calcium balance in the endoplasmic reticulum [29]. Concerning that, FAIM2 KO mouse models showed increased apoptosis of newly activated T lymphocytes; this evidence suggests that FAIM2 protects the T cells from apoptosis, as the first days of a cytotoxic response are hostile for memory T cells [30].

We employed the GWAS Atlas data to perform a second PheWAS to compare our observations with previous results from other populations. Notably, this second PheWAS showed that SNP rs7305229 has not been associated with immunologic traits. Instead, this SNP is mainly associated with metabolic traits. In this sense, FAIM2 is mainly expressed in nervous system tissues [29]. Remarkably, Littleton et al. (2024) reported that variants of FAIM2 in 3’UTR change expression levels and influence neuronal populations in mouse models. They propose that expression of FAIM2 is required for the normal development of anorexigenic neurons in the hypothalamus; however, individuals with less expression would have fewer anorexigenic neurons and thus increased appetite [27].

There is limited evidence on the effect of genetic variations of FAIM2 on GADA antibodies, even though GADA is present in other autoimmune disorders. In this sense, GADA antibodies are detectable in some individuals with Graves' disease; furthermore, these subjects with GADA can develop T1D [31,32]. Notably, a study associated a SNP of FAIM2 with antibody levels in patients with Graves’ disease [28]. Nonetheless, other variants of FAIM2 in the 3’UTR region are associated with obesity and body adiposity in other GWAS [33,34]. Nevertheless, we did not find associations between BMI with GADA and rs7305229. Littleton et al. (2024) propose that subjects with fewer anorexigenic POMC neurons would have reduced expression of FAIM2, induce increased appetite, and higher risk of becoming overweight in childhood [27]. Interestingly, evidence shows that childhood obesity genetics are a risk factor for diabetes regardless of autoimmunity or insulin resistance [35]. In addition, high BMI is a risk factor for childhood type 2 diabetes and latent autoimmune diabetes in youth [36]. Furthermore, obesity in childhood also is a risk factor in adult age for T2D and LADA [37].

We hypothesize that our subjects with LADA and rs7305229 would have lower GADA antibody levels. Likewise, other SNPs of the FAIM2 gene could produce increased BMI and thus, a high risk of diabetes. BMI and GADA antibody levels are heterogeneous features among patients with LADA: some individuals have high antibody levels and low BMI, while others have low antibody levels and high BMI [14]. In this regard, Caucasian patients with LADA have a greater proportion of high antibody levels than Chinese subjects with LADA. Conversely, metabolic syndrome is more frequent in the Caucasian population than in the Chinese individuals. Another difference is insulin resistance; LADA has similar resistance to T2D in the Chinese population, but in Caucasian patients, it is less pronounced than in T2D [12]. Our population showed low antibody levels, high frequency of metabolic syndrome, and insulin resistance like T2D. The variability of clinical features of LADA would reflect the differences in genetic and environmental risk factors among individuals and populations [14].

We know, this is the first GWAS of LADA performed in the Mexican population; nevertheless, it has some limitations. First, all subjects with LADA were from a single care center in Southeastern Mexico. In addition, the healthy subjects’ genotypes were from a nationwide database, which limited the comparison between phenotypical features. The sample size of our study is small, and we do not have a replication cohort, reducing our statistical power. Although using the subpopulation from Southeast Mexico limited the sample size, our GWAS was well controlled, as reflected in our lambda values.

4. Materials and Methods

4.1. Sample Population

We included a subsample of 59 deeply phenotyped individuals diagnosed with LADA from a previous study [24]. These subjects attend the Diabetes Clinic of the Hospital Regional de Alta Especialidad “Dr. Gustavo A. Rovirosa Pérez” from January 2020 to May 2021. The hospital cares for patients from Chiapas, Veracruz, and Tabasco (Southeast Mexico). A diabetologist evaluated all recruited subjects. We included individuals who met the following criteria: Mexican descent, older than 30 years at the time of diabetes diagnosis, more than six months without insulin treatment after diagnosis, and positivity to GADA determined by ELISA (≥ 5 UI/mL). We followed the ethical principles from the Declaration of Helsinki. This study was approved by the Ethics Committee of the High Specialization Regional Hospital “Dr. Gustavo A. Rovirosa Pérez” (HR/ENS/ARM/8073/2023; July 5th, 2023). They received verbal and written information about this study. Every patient participated voluntarily without receiving remuneration and signed an informed consent.

We also included a population-based control sample from the Mexican Genomic Database for Addiction Research (MxGDAR/Encodat) [38]. Creation of the MxGDAR/ENCODAT database was approved by the Research Ethics Committees of the Instituto Nacional de Psiquiatría Ramón de la Fuente Muñiz (CEI/C/083/2015) and the Instituto Nacional de Medicina Genómica (01/2017/I). This study was derived from the Mexican National Survey of Tobacco, Alcohol, and Drug Use. Subjects from the 32 states in Mexico were sampled, and regions where more than half of the population speaks a Native American language were excluded. We included 3121 individuals who met the following inclusion criteria: Mexican descent, older than 18 years, without self-reported diabetes or other comorbidities.

4.2. Phenotypic Characterization

A structured questionnaire was applied to recruited patients with LADA, we collected sociodemographic data (age, gender, education, marital status, socioeconomic status). Family history, age at diabetes diagnosis, diabetes complications, comorbidities, and pharmacotherapy were retrieved from clinical records. We measured weight, height, BMI, waist circumference, and hip circumference. Fasting glucose, total cholesterol, HDL cholesterol, LDL cholesterol, and triacylglycerols were determined from blood serum. GADA antibody levels were detected with an ELISA kit (MyBioSource, USA). A total of 51 phenotypes were evaluated.

4.3. Microarrays Analysis

We extracted DNA from whole blood samples using a salting-out procedure with the Puregene Blood kit (Qiagen, USA) following the manufacturer’s instructions. DNA was hybridized with the Infinium Psycharray Beadchip (Illumina, USA) following the automated protocol. All genotyping procedures were performed in the Unidad de Alta Tecnología para Expresión y Microarreglos of the Instituto Nacional de Medicina Genómica. The microarray analysis was similar to the population-based controls from the MxGDAR database, as previously described [38].

4.4. Genotype Quality Control

Quality control of genotypes was performed with PLINK 1.9 version b7.1 [39]. We excluded genetic variants with call rate < 95%, minor allele frequency (MAF) < 05, p-value <1e-10 for a chi-squared test for Hardy-Weinberg equilibrium for cases and p < 1e^-06 for controls, and chimeric alleles (AT and CG). Individuals were excluded if they had a call rate < 95% for all the genetic variants and ± 3 standard deviations from the mean heterozygosity rate. We calculate the identity-by-state of all sample pairs to identify genetic relatedness. Pairs of individuals with pi-hat > 0.2 were selected, and related individuals with lower genotyping rates were excluded. Using this threshold, we kept 59 individuals diagnosed with LADA and 3086 population-based controls for further analysis.

4.5. Population Stratification

We used the pcair function of the GENESIS package version 2.30 for R version 4.2 to estimate principal components analysis (PCA) to evaluate population stratification. We included independent SNPs filtered by linkage disequilibrium pruning (threshold: 0.1, slide size: 106). Nonetheless, the population analysis showed the previously reported differences between the Southeast and other regions of Mexico [40]. As previously described, the population from the South of Mexico could have a genetic structure different from that of the country's center [38]. To avoid effects in our statistical tests for population stratification, we cluster cases and controls. We visualize the potential clustering of subpopulations inside our sample with Uniform Manifold Approximation and Projection (UMAP) using the first 5 principal components (PCs). Two clusters were observed: one group with individuals from various regions of Mexico (cases n=12, controls n=2592), and the other with subjects mainly from Southeastern Mexico (cases n=50, controls n=491). The subpopulation that had the most cases from our LADA sample was selected. (Supplementary Figure 2)

4.6. Genotype Imputation

We performed the imputation of autosomal genetic variants using the TOPMed Imputation Server with the TOPMed-r2 reference panel [41]. We used the following parameters as post-imputation quality control, keeping those genetic variants with MAF > 0.05, quality score (R²) > 0.7, and chi-squared p > 1x10^-6 for Hardy-Weinberg equilibrium. After quality control, we kept 5,797,363 genetic variants for the following analysis.

4.7. Phomene-Wide Association Study (PheWAS)

GENESIS package for R was employed to perform a case-control GWAS of LADA (n=50) versus controls (n=491), including the first 5 PCs, sex, and age as covariates. The statistic inflation factor (λ) was 1.082. Also, we performed genome-wide linear regressions between SNPs and phenotypic features of LADA cases with PLINK 1.9; covariates for this analysis were the first 2 PCs, sex, and age. Included phenotypes were continuous variables (anthropometric measures, biochemical measures, antibody levels, and years of disease evolution). Categorical variables were excluded from PheWAS as they had poor quality control results. A p-value < 5x10^-8 was considered genome-wide significant and a p-value < 1x10^-5 was nominally significant for both analyses. QQ plots and Manhattan plots were realized with ggplot2 and fastman libraries for R.

4.8. Fine Mapping

We used the susieR package version 0.12.35 for R version 4.2 was employed to perform a fine mapping for the associated loci using the summary statistics that showed GWAS significance. A 1000 kb window around the top genetic variant was selected as a significant locus. The LDlinkR package for R was employed to acquire the LD correlation matrix from the LDlink web application, we used the Mexican Ancestry from Los Angeles (MXL) population as a reference [42]. Genetic variants with 95% posterior inclusion probability (PIP) were considered part of the credible set.

4.9. PheWAS Analysis

To find associations of SNPs with phenotypes in previous GWAS, the GWAS Atlas website was employed for PheWAS analysis of the nominally significant SNP. Associations with p-value <5x10^-8 were considered significant.

4.10. Conditional Analysis

We realized a conditional analysis to identify phenotypes dependently associated with significant GWAS loci. For patients with LADA, a second genome-wide analysis was performed with PLINK 1.9 [43]. This analysis included body mass index as a covariant for the significant GWAS loci. Also, a linear model for body mass index and GADA levels was calculated with R.

5. Conclusions

We performed a phenome-wide association study of LADA in subjects from southeastern Mexico, and the results indicated that rs7305229 was associated with GADA antibody levels. These findings contribute to understanding the genetic mechanisms underlying LADA, as its heterogeneity among individuals may be reflected in their genetic constitutions. Because the genetic risk factors for LADA may vary among populations, it is necessary to expand research in this area.