A cohort of 376 Mexican MEZ patients with childhood-onset SLE was recruited from two tertiary level institutions in Mexico City. From this cohort, a set of 225 trio families was used for transmission disequilibrium test (TDT) analysis. All patients were younger than 16 years old, and met the American College of Rheumatology (ACR) SLE classification criteria (18). The control group included 400 healthy sex-matched adults from Mexico City (MEZ). To compare the allele frequencies of the 15 SNVs from the Xq28 region of both the MEZ and Amerindian (AM) populations with those reported for other populations in the 1,000 Genomes database (1KGP) phase 3, we used imputed whole-exome sequencing from 1,074 AM individuals from the Metabolic Analysis in an Indigenous Sample (MAIS) cohort (median Amerindian ancestry of 0.96) (19, 20). To validate the imputed risk alleles, we genotyped tag SNVs in 300 AM individuals belonging to four of the most representative Mexican ethnic groups: Tarahumara from the North (N), Nahuatl from the Central East (CE), Zapoteco from the South (S), and Mayan from the Southeast (SE).

This study was conducted in accordance with the Declaration of Helsinki. All individuals signed an informed consent, and the study was approved by local ethics and research committees from all participating institutions. Parents provided consent for child participation, and all children older than 8 years old assented.

Genomic DNA was isolated from whole blood samples using the Gentra-Puregene kit (Qiagen, Valencia, CA, United States). We genotyped 15 SNVs using TaqMan Assays on Demand (Applied Biosystems, Foster City, CA, United States). The genotyped SNVs were rs2266890 and rs13397 located in TMEM187; rs3027898, rs2239673, rs763737, rs5945174, rs7061789, and rs1059702 located in IRAK1; and rs2075596, rs3027933, rs17435, rs1734787, rs1734791, rs1734792, and rs2239464 located in MECP2. The call rate was >98% for all SNVs. Genotypes were confirmed by randomly sequencing 10% of the samples using an automated ABI PRISM 310 Genetic Analyzer (Applied Biosystems, Life Technologies, CA, United States), which showed 100% reproducibility. To estimate ancestry in 320 SLE patients and 400 controls, we used the 6.0 SNV array (Affymetrix, Santa Clara, CA, United States) or GoldenGate genotyping assay (Illumina, San Diego, CA, United States), which contained 96 ancestry-informative markers (AIMs) validated in previous studies (21). The median Amerindian ancestry was 0.72 for cases and 0.66 for controls.

The allele frequencies were compared between Mexican populations and other continental populations taken from the 1KGP phase 3, using a chi-square test in R v.3.1 (http://www.r-project.org). To validate the frequencies of the SNVs in Indigenous individuals, we selected rs1059702, rs2075596, rs2239464, rs2266890, and rs3027898, which we predicted to be tag SNVs. For the SLE case-control association, logistic regression was performed. The statistical power was computed using the Genetic Power Calculator program (22), assuming a disease prevalence of 0.06% (23) and a significance level of 0.05. The association was evaluated under an additive model, and was adjusted by sex and ancestry, using the first two principal components. A p value of <0.05 after Bonferroni multiple testing correction was considered statistically significant. All analyses were performed using PLINK v.1.9 software (24). Haploview v.4.2 program was used for LD calculation through pairwise r² values, allele, and haplotype frequencies, and to construct haplotypes (25). LD plots were constructed for the Mexican Indigenous and MEZ populations, and imputation was performed using IMPUTE 2 (26), using the 1KGP phase 3 as a reference panel. Variants with an information score of >0.9 were included in the analysis.

Total RNA was obtained from PMBCs collected from eight adult healthy control carriers of the Xq28 risk haplotype (H1) and 12 non-carriers, using Trizol reagent (Invitrogen, Carlsbad, CA, United States) according to the manufacturer’s instructions. RNA integrity was assayed by agarose gel electrophoresis. For the mRNA reverse transcription, 500 ng of total RNA, and oligo dT₁₈ primer were used, along with the RevertAid enzyme (Thermo Fisher Scientific, Massachusetts, United States). We measured IRAK1 expression using the probe Hs01018347_m1, and we measured the two MECP2 isoforms (A and B) using the probes Hs00172845_m1 and Hs01598237_m1 (Thermo Fisher Scientific, Massachusetts, United States), respectively. The data were normalized using the GAPDH Endogenous Control TaqMan probe (Thermo Fisher Scientific, Massachusetts, United States). All qPCR assays were performed using a ViiA7™ Real Time system thermocycler (Thermo Fisher Scientific, Massachusetts, United States). Relative mRNA expression levels were calculated using the 2^–ΔCt method (27). Statistical analysis was performed using Student’s t-test, with GraphPad Prism software v.9.0 (GraphPad, San Diego, CA, United States). A p value < 0.05 was considered to indicate significant differences.

Based on the sample size, a statistical power of 80% was reached. The case-control analysis for each of the 15 analyzed single variants in the Xq28 region revealed that all of the tested alleles were significantly associated with childhood-onset SLE, with ORs ranging from 1.49 to 1.75 after correction for ancestry and sex (Table 1). We subsequently replicated this association using TDT analysis, by genotyping the 15 SNVs in 225 trios. We found that all risk alleles were over-transmitted, although only two SNVs located on IRAK1 (rs763737 and rs59455174) and three located on MECP2 (rs2075596, rs3027933, and rs17435) were significantly associated with childhood-onset SLE (Table 2). Furthermore, only the haplotype H1 was associated with susceptibility to childhood-onset SLE (OR = 1.97, p = 4 × 10^–6) in the Mexican population (Table 3).

We compared the allele frequencies of these SNVs among our 400 MEZ controls with the frequencies reported for other populations in the 1KGP phase 3. Only the South and East Asian populations exhibited greater allele frequencies than the MEZ population from our present study. The frequencies of all the analyzed risk alleles within our sample were higher compared with those reported for persons of Mexican Ancestry from Los Angeles, California, USA (MXL); were four times greater than in the Iberian Population in Spain (IBS), and were even greater than those reported in the African population (AFR) (Table 4 and Supplementary Table 1). To confirm whether these risk alleles were inherited by the MEZ from the Indigenous population, we used imputed whole-exome sequencing from Indigenous populations belonging to the MAIS Cohort (median Amerindian ancestry of 0.96). Genotype validation was performed using tag SNVs (Supplementary Table 2) in AM individuals belonging to four of the most representative Mexican ethnic groups: Tarahumara, Nahuatl, Zapoteco, and Mayan (which show a concordance of >99%). We found that the allele frequencies of the risk variants were all significantly enriched in the AM population, being some of the highest frequencies worldwide. Our analysis also revealed an increasing gradient from the north to the south of the Mexican territory (Supplementary Table 3).

Evaluation of the haplotype H1 structure in the MEZ and AM populations revealed that the SNVs located in IRAK1–MECP2 were in strong LD, while those located in TMEM187 exhibited lower LD values. Similar results were found in all Latin American populations and in East Asians (EAS) from the 1KGP phase 3. In contrast, in the IBS and CEU populations, the risk haplotype maintained a strong LD, including the variants in TMEM187 (r²: IBS = 85 and 92; CEU = 81 and 69). Interestingly, African populations (YRI and GWD) did not show LD among most of the risk alleles in the Xq28 region, except in five (rs3027898, rs2239673, rs763737, rs5945174, and rs7061789) of the seven analyzed SNVs in IRAK1, which were observed to be in complete LD in all populations (Figure 1).

Because of these specific population differences in the structure of the Xq28 region, we were able to observe a variable number of haplotypes in this region with a frequency ≥ 0.01. For example, five haplotypes were observed in the MEZ and MXL populations, four in AM, three in IBS, and 14 in African populations. The risk haplotype H1 was found at a frequency of 0.61 in AM, 0.49 in MEZ, and 0.13 in IBS, while it was absent in African populations, such as the Mende in Sierra Leona (MSL) and Yoruba in Ibadan Nigeria (YRI) (Figure 2).

To examine the effect of the risk haplotype on the expression of IRAK1 and the two MECP2 isoforms in healthy control H1 haplotype carriers and non-carriers, we performed real-time quantitative PCR (qPCR). Notably, the levels of both the MECP2A and B transcript isoforms were significantly lower (40%) in the risk haplotype carriers. The analysis of IRAK1 transcripts did not show significant differences between carriers and non-carriers (Figure 3).