Populations from Indian Himalayan and adjoining regions share their associated lactase persistence variants with the Tibetans and South Asians

Observed genotype count of putative LP-associated variants in the HAAPs

A genetic screening of the 100 bp region located 14 kb upstream of the LCT gene led to the observation of four polymorphic sites, namely, rs4988235 (-13910C >T), rs1575359915 (-13838G>A), rs1679771596 (-13906T>A), and rs4988236 (-13908C >T) (Supplementary Table 4). It was interesting to note that previously reported LP-associated variants in Sub-Saharan Africa and Saudi Arabia, respectively, rs41525747 (-13907C>G) and rs41380347 (-13915T>G), were absent in the HAAPs. Analysis of -13910C>T, attributed to the LP majorly among Northern Europeans, Americans, and Indians revealed the absence of homozygous derived genotype (TT) in all studied populations. The heterozygous (CT) genotype for the same was observed in 16 individuals among 9 HAAPs. This includes four individuals from Bhutia2, three from the Darjeeling mixed population (a mixture of ethnic groups from Darjeeling Hills) (Supplementary Table 1), two individuals each from Khasi and Sherpa populations, and a single individual each from Limbu, Rai, Tamang, Ao Naga, and Mizo populations (Table 1). The -13838G>A variant, which was earlier reported in the Tibetan population from the Tibetan Plateau as the predominant LP-associated variant, was present in 19 individuals from 7 HAAPs in heterozygous condition (GA). This included seven individuals from the Tibetan2, four individuals from Bhutia2, three from the Darjeeling mixed, two from the Sherpa, and one each from Tamang, Nyishi, and Khasi populations (Table 1). None of the individuals carried the homozygous derived genotype (AA). The other two Tibetan-specific putative LP variants -13906T>A and -13908C>T, also earlier reported in Tibetans³⁷ were observed in heterozygous conditions (TA and CT, respectively) in two individuals (one from Ao Naga and one from Tibetan2) and a single individual from the Darjeeling mixed population, respectively, among the HAAPs (Table 1).

Derived allele frequency of putative LP-associated variants in HAAPs

Although the nationwide frequency of the -13910C>T variant was found to be 10.3%¹⁶the HAAPs had a lower frequency of 1.18%, including 1.13% in TB and 1.61% in AA-speaking groups (Supplementary Table 2). It was observed to be highest among Bhutia2 (4.76%), followed by Tamang (4.17%), Sherpa (3.13%), Darjeeling mixed (3%), Limbu (1.92%), Rai (2.78%), Khasi (2.78%), Mizo (1.43%), and Ao Naga (0.81%) populations (Table 1). On the other hand, for -13838G>A, the derived allele frequency was highest in Bhutia2 (4.76%), followed by Tamang (4.17%), Tibetan2 (4.07%), Sherpa (3.13%), Darjeeling mixed (3%), Nyishi (1.47%), and Khasi (1.39%) populations (Table 1). Overall, the derived allele frequency for -13838G>A was 1.4% in HAAPs, which constituted 1.46% in TB and 0.81% in AA-speaking groups (Supplementary Table 2). For -13906T>A, the derived allele was found only in two TB-speaking populations, namely Ao Naga (0.81%) and Tibetan2 (0.58%), altogether accounting for 0.15% in the total HAAPs and 0.16% in TB-speaking HAAPs (Table 1Supplementary Table 2). The derived allele for the variant -13908C>T was observed only in the Darjeeling mixed population at a frequency of 1%, accounting for 0.07% of total HAAPs and 0.08% in the TB-speaking HAAPs (Table 1Supplementary Table 2). We detected no HAAPs deviating from the Hardy-Weinberg equilibrium (HWE) for observed variants (Table 1). Sensitivity analysis removing small populations (N < 10 and Tamang; N = 12) revealed, the derived-allele frequency f(T) for -13910C>T 1.18% in overall HAAPs, compared to 1.2% estimated after excluding populations with N < 10 and 1.14% after successive exclusion of Tamang (N = 12) along with populations with N < 10. Similarly, the frequency of -13838G>A was 1.4%, 1.42%, and 1.37%, respectively. The two rare variants, -13906T>A and -13908C>T, also showed minimal changes across the sensitivity datasets (Supplementary Table 6). The median-joining network tree suggested that the two variants -13838G>A and -13910C>T have been found separately in different individuals of the HAAPs, and no single individual was found to carry both variants together. Bhutia2, Tamang, Sherpa, Darjeeling mixed, and Khasi populations showed the presence of both variants (Supplementary Fig. 1). The frequency estimates of 95% CIs overlap substantially, indicating that the observed difference is small relative to the uncertainty around the frequency estimates. However, because overlapping confidence intervals alone do not constitute a formal test of equality, we additionally performed a paired comparison of -13910C>T and -13838G>A using an exact McNemar test that revealed no statistically significant differences between the tested pair (exact McNemar test, P = 0.736; OR = 1.19, 95% CI: 0.58–2.47).

Predicted frequency of candidate LP-associated genotypes in the HAAPs

The present study does not assess the LP phenotypes coupled with genetic variants. Rather, it calculates the predicted frequency of candidate LP-associated genotypes in HAAPs using the q² + 2pq formula¹⁶. As none of the individuals with the homozygous derived genotype was observed for the studied variants, the predicted frequency of candidate LP-associated genotypes in the HAAPs was contributed only by heterozygous genotypes of -13838G>A (2.77%) and -13910C>T (2.34%). The expected frequency of derived allele carriers in the TB-speaking HAAPs for the above two loci was 2.89% and 2.25%, respectively, and for AA-speakers 1.61% and 3.2%, respectively (Supplementary Table 2). We observed -13838G>A is responsible for 9.3% predicted frequency of candidate LP-associated genotypes in Bhutia2, followed by 8.16% in Tamang, 7.97% in Tibetan2, 6.15% in Sherpa, 5.91% in Darjeeling mixed, 2.92% in Nyishi, and 2.76% in Khasi population. On the other hand, variant -13910C>T was responsible for 9.3% of predicted frequency of candidate LP-associated genotypes in Bhutia2, followed by 8.16% in Tamang, 6.15% in Sherpa, 5.91% in Darjeeling mixed, 5.48% in both Khasi and Rai, 3.81% in Limbu, 2.84% in Mizo and 1.61% in Ao Naga population (Table 1). 11 populations out of 19 HAAPs showed the presence of predicted candidate LP-associated genotypes of -13838G>A and -13910C>T variants in them (Figs. 1, 2), In Bhutia2, Tamang, Darjeeling mixed, and Sherpa, predicted candidate LP-associated genotypes were present in equal frequency for both -13838G>A and -13910C>T. Whereas, the Khasi population showed a higher proportion of their predicted LP-associated genotypes contributed by -13910C>T compared to -13838G>A. The predicted frequency of candidate LP-associated genotypes in Tibetan2 and Nyishi was also solely contributed by -13838G>A. On the other hand, in Rai, Limbu, Ao Naga, and Mizo it was exclusively contributed by -13910C>T (Fig. 2). It is noteworthy that Banravat, Lepcha, Bhutia1, Tibetan1, Naga Sema, Poumai Naga, Chakhesang Naga, and Santhal populations had ancestral alleles for both loci.

The predicted frequency of candidate LP-associated genotypes with 95%CI in the Himalayan and adjoining populations (HAAPs), estimated following Romero et al. 2012.

Pan-India comparison of HAAPs based on linguistic affiliation

A comparison of genotype frequency across India for -13910C>T revealed the highest frequency of heterozygotes (CT) (26.23%) in IE, followed by DR (11.96%), AA (2.29%), and TB (2.17%). The homozygous-derived genotype (TT) was found only in IE (5.25%) and DR (1.64%). On the contrary, the heterozygous and homozygous derived genotypes for -13838G>A (GA, AA), -13906T>A (TA, AA), and -13908C>T (CT, TT) were absent in the IE and DR. The frequency of heterozygotes for the -13838G>A was found to be 0.25% for AA speakers. The TB speakers were heterozygotes for all the above-mentioned three variants with frequencies of 2.44%, 0.27%, and 0.14%, respectively (Table 2Supplementary Table 3).

A similar trend was noted in the distribution of the derived allele frequency of -13910C>T; the highest frequency in IE (18.36%) followed by DR (7.62%), AA (1.15%), and TB (1.08%). IE and DR lacked the derived alleles for the other three SNPs (Table 2). Whereas AA speakers had 0.13% derived allele frequency for -13838G>A, the TB speakers possessed the derived alleles for all three loci, i.e., -13838G>A (1.22%), -13906T>A (0.14%), and 13908C>T (0.07%) (Table 2Supplementary Table 3). The predicted frequency of candidate LP-associated genotypes in the IE and DR was 33.35% and 14.66%, respectively, and was contributed by the homozygous derived and heterozygous genotypes of -13910C>T. Whereas in AA speakers, it was contributed by heterozygous genotypes of -13910C>T (2.28%) and -13838G>A (0.25%). Interestingly, in Tibetans, heterozygous genotypes of both -13910C>T and -13838G>A contributed to 2.16% and 2.42%, respectively (Table 2). None of the individual Indian populations within each linguistic group showed deviation from the HWE (Supplementary Table 3). However, upon merging data based on their linguistic affiliation, the IE (p-value = 0.0003) and DR (p-value = 0.0001) deviated from HWE (Table 2). Sensitivity analysis removing the small populations with fewer than 10 individuals revealed that all of the linguistic groups have relatively consistent 95%CIs for both derived allele frequency and predicted frequency of candidate LP-associated genotypes in them (Supplementary Table 7). The Kruskal-Wallis tests revealed that IE and TB are significantly different from the other linguistic families for -13910C>T and -13838G>A variants, respectively. Comparison of TB-speaking HAAPs with the other three linguistic groups also revealed similar results to the pan-Indian TB group (Table 5a). Beta-binomial analysis also replicated the similar trend for -13910C>T (Supplementary Table 5b). The derived allele frequency for the -13910C>T (Fig. 3a) variant in the Indian subcontinent revealed a general North-West to South-East declining trend. The derived allele for variant -13838G>A was found only in the HAAPs and the Tibetan refugees (Fig. 3b).

Frequency distribution of -13910T and -13838 A in India.

As Indian HAAPs share their close genetic affinity with the Tibetans and East Asians, we compared the derived allele frequencies of the four variants obtained in the present study with the CHB and JPT data retrieved from the 1000 Genomes Project’s from Ensembl. However, for -13838G>A and -13906T>A there was no data available for the CHB or JPT populations in the 1000 Genome. For -13910C>T, the derived allele frequency, f(T), was nil in both CHB and JPT populations. However, the derived allele frequency of this variant was 0.0118 in HAAPs and 0.0113 in TB-speaking HAAPs. For -13908C>T, the derived allele frequency, f(T), was 0.015 in CHB and 0.00 in the JPT population. The frequency of this variant was 0.0007 in HAAPs and 0.0008 in the TB-speaking HAAPs which is quite low.