Rare, Protein-Altering Variants in AS3MT and Arsenic Metabolism Efficiency: A Multi-Population Association Study.

BACKGROUND: Common genetic variation in the arsenic methyltransferase (AS3MT) gene region is known to be associated with arsenic metabolism efficiency (AME), measured as the percentage of dimethylarsinic acid (DMA%) in the urine. Rare, protein-altering variants in AS3MT could have even larger effects on AME, but their contribution to AME has not been investigated.
OBJECTIVES: We estimated the impact of rare, protein-coding variation in AS3MT on AME using a multi-population approach to facilitate the discovery of population-specific and shared causal rare variants.
METHODS: We generated targeted DNA sequencing data for the coding regions of AS3MT for three arsenic-exposed cohorts with existing data on arsenic species measured in urine: Health Effects of Arsenic Longitudinal Study (HEALS, n=2,434), Strong Heart Study (SHS, n=868), and New Hampshire Skin Cancer Study (NHSCS, n=666). We assessed the collective effects of rare (allele frequency <1%), protein-altering AS3MT variants on DMA%, using multiple approaches, including a test of the association between rare allele carrier status (yes/no) and DMA% using linear regression (adjusted for common variants in 10q24.32 region, age, sex, and population structure).
RESULTS: We identified 23 carriers of rare-protein-altering AS3MT variant across all cohorts (13 in HEALS and 5 in both SHS and NHSCS), including 6 carriers of predicted loss-of-function variants. DMA% was 6-10% lower in carriers compared with noncarriers in HEALS [β=-9.4 (95% CI: -13.9, -4.8)], SHS [β=-6.9 (95% CI: -13.6, -0.2)], and NHSCS [β=-8.7 (95% CI: -15.6, -2.2)]. In meta-analyses across cohorts, DMA% was 8.7% lower in carriers [β=-8.7 (95% CI: -11.9, -5.4)]. DISCUSSION: Rare, protein-altering variants in AS3MT were associated with lower mean DMA%, an indicator of reduced AME. Although a small percentage of the population (0.5-0.7%) carry these variants, they are associated with a 6-10% decrease in DMA% that is consistent across multiple ancestral and environmental backgrounds. https://doi.org/10.1289/EHP8152.

Introduction

More than 200 million people are exposed to inorganic arsenic (iAs) through drinking water worldwide (Naujokas et al. 2013). Dietary exposure to iAs (primarily through rice and grain products) is an emerging concern that has received less regulatory focus (Nachman et al. 2017). Chronic exposure to levels of iAs above the World Health Organization (WHO) safety standard for drinking water () has been recognized to pose a significant risk of adverse outcomes across multiple organ systems (Mohammed Abdul et al. 2015). Epidemiological studies in arsenic-affected areas of South America, Asia, and North America have demonstrated that chronic exposure to iAs is associated with adverse health effects, including increased risk for cardiovascular disease (Moon et al. 2017), diabetes (Sung et al. 2015), cognitive dysfunction (Karim et al. 2019; Tyler and Allan 2014), adverse birth outcomes (Milton et al. 2017), and overall mortality (Argos et al. 2010). In addition, iAs exposure increases the risk for cancers of the skin (Karagas et al. 2015), lung (Kuo et al. 2017; Lamm et al. 2015), bladder (Gamboa-Loira et al. 2017; Kuo et al. 2017), and kidney (Ferreccio et al. 2013).

The metabolism of iAs in humans involves a series of reduction and methylation reactions. iAs enters the body as (trivalent arsenite) or (pentavalent arsenate). Sequential reduction and methylation reactions are catalyzed by glutathione and arsenic ( oxidation state) methyltransferase (AS3MT), respectively, producing monomethylated (MMAIII and MMAV) and dimethylated (DMAIII and DMAV) forms of arsenic. Consumed arsenic is eliminated in the urine, primarily as a DMA, although smaller percentages are eliminated as MMA and iAs. Arsenic methylation facilitates excretion of arsenic in urine given that DMA is more rapidly expelled from the body compared with MMA or iAs (Gamble et al. 2006, 2007; Peters et al. 2015). There are a number of factors believed to impact individuals’ ability to metabolize arsenic, including genetic differences, age, sex, body mass index (BMI), smoking status, nutritional status, and arsenic exposure level (Jansen et al. 2016; Kordas et al. 2016; Shen et al. 2016). Arsenic metabolism efficiency (AME) is often represented by the percentage of arsenic species in the urine that are DMA (DMA%) (Hopenhayn-Rich et al. 1996; Vahter 1999).

Common inherited genetic variation in the 10q24.32 region (containing AS3MT) is known to impact AME and the risk for arsenic-induced skin lesions. Two independent association signals for DMA% have been reported in the 10q24.32 region in a Bangladeshi cohort (Pierce et al. 2012, 2013). These DMA%-associated variants, best represented by single nucleotide polymorphisms (SNPs) rs9527 and rs11191527, were also associated with MMA%, iAs%, and skin lesion risk, highlighting the 10q24.32 region as a key source of individual variability in AME and susceptibility to arsenic toxicity. Studies in other arsenic-exposed populations have also reported associations between AS3MT SNPs and AME phenotypes (Agusa et al. 2009; Engström et al. 2011, 2015; García-Alvarado et al. 2018), including a study of American Indian populations in the United States with low-to-moderate levels of iAs exposure (Balakrishnan et al. 2017).

Rare genetic variation [i.e., minor allele frequency (MAF) ] may also impact AME. Gao et al. (2015) estimated the heritability of DMA% due to common SNPs to be [95% confidence interval (CI) (, 39.5%) based on a standard error (SE) of 12%] in a sample of unrelated Bangladeshi individuals; however, when restricting to only close relatives, the heritability estimate increased to 63% [95% CI: (31.6%, 94.4%), based on a SE of 16%], potentially reflecting the contributions of rare variants. Family studies in American Indian communities also reported DMA% heritability estimates (Tellez-Plaza et al. 2013), also suggesting that rare variants may contribute to the heritability of AME. However, the contribution of rare variants to AME has not been assessed for any specific genes, including AS3MT.

In this study, our primary goal was to characterize the impact of rare, protein-altering variants in AS3MT on AME in three different arsenic-exposed populations. This multicohort approach allows us to assess the generalizability of our findings across cohorts of different ancestral and environmental backgrounds. Understanding the role of rare variation is critical for identifying individuals at high risk for arsenic toxicities and understanding the biological mechanisms underlying interindividual differences in susceptibility to toxicity.

Methods

Study Participants

This study uses data from three different studies of arsenic-exposed populations: the Health Effects of Arsenic Longitudinal Study (HEALS), the Strong Heart Study (SHS), and the New Hampshire Skin Cancer Study (NHSCS). We selected these three cohorts because of their known arsenic exposure, existing data on arsenic species, and available DNA for sequencing.

HEALS is a prospective cohort study designed to investigate health outcomes associated with chronic arsenic exposure through drinking water in Araihazar, Bangladesh (Ahsan et al. 2006). Arsenic concentrations have been measured in wells in the study area. A total of 11,746 participants (5,042 males and 6,704 females, 18–75 years of age) were enrolled between 1999 and 2001 after providing informed consent. Participants completed a written questionnaire and participated in clinical exams at baseline. Blood and urine samples were also collected from participants at baseline (Ahsan et al. 2006). Participants were simultaneously assessed for arsenical skin lesions at baseline and every 2 y thereafter by trained physicians using a structured protocol (Argos et al. 2011). Arsenic species in baseline urine samples were measured previously for 4,794 HEALS participants (Jansen et al. 2016). For the present study, we selected 2,719 of these participants who had available DNA for sequencing, of which 273 were skin lesion cases at baseline. Although many of these participants were randomly selected for arsenic species measurement in baseline urine, a sizable fraction () were selected at subsequent follow-up visits based on outcomes they experienced (i.e., skin lesions, respiratory symptoms, and cardiovascular conditions) in a case–cohort fashion. So although this group was relatively healthy at baseline, the participants were not randomly selected.

The SHS is the largest population-based cohort study of cardiovascular disease in American Indian men and women. The SHS includes 12 American Indian tribes and communities and was designed to estimate cardiometabolic disease morbidity and mortality and the prevalence of its risk factors (Moon et al. 2013). Between 1989 and 1991, the SHS recruited 4,549 American Indian men and women between the ages of 45 to 74. These participants have been exposed to low-to-moderate levels of iAs primarily through drinking water (Navas-Acien et al. 2009) and to a lower extent, diet (i.e., rice) (Nigra et al. 2019). Arsenic species were measured in 3,973 participants (Moon et al. 2013). These individuals are members of large families, so we restricted participation to a group of 997 unrelated individuals with existing measures of arsenic species and available DNA for sequencing.

The NHSCS is a population-based case–control study of basal cell carcinoma (BCC) and squamous cell carcinoma (SCC) of the skin (Gilbert-Diamond et al. 2013). Invasive SCC incident cases (), 25–74 years of age, were recruited from of practicing dermatologist and pathologist clinics in New Hampshire and bordering areas. Controls () frequency-matched on sex and age were selected from the New Hampshire Center for Medicare and Medicaid Services and driver’s license records. The enrollment period for cases and controls was between 2003 and 2009. Participants interviewers masked to the study hypothesis ascertained sociodemographic, lifestyle, medical, and sun-exposure information. In addition, home water and spot urine samples were collected and used to measure total urinary arsenic concentration and arsenic species (post-diagnosis measurements for SCC cases) (Gilbert-Diamond et al. 2013). Of the 993 cases and controls, 288 individuals did not have sufficient DNA for sequencing and were excluded from the present study. This resulted in 706 individuals with existing data on arsenic species selected for the present sequencing study (349 SCC cases and 357 controls).

Measurement of Arsenic Metabolites in Urine

All participants selected for this study had existing data on arsenic species in urine. In all three cohorts, speciation analysis of arsenic metabolites was performed using high-performance liquid chromatography (HPLC) (Scheer et al. 2012) followed by detection using ICPMS (Ahsan et al. 2006; Gilbert-Diamond et al. 2011; Moon et al. 2013). Details regarding the limit of detection (LOD) for each metabolite and percentage of samples below the LOD have been described previously (Ahsan et al. 2006; Gilbert-Diamond et al. 2011; Navas-Acien et al. 2009). Briefly, in HEALS the LOD was for , , MMA, and DMA. In SHS, was first oxidized to in the urine, following the methods described by Scheer et al. (2012), this method minimizes undetectable or in urine samples from populations with low levels of exposure. The LOD was for and for MMA and DMA in SHS. In NHSCS, the LOD was for , for , for MMA, and for DMA.

For this analysis, and were summed to obtain total iAs, and each arsenic species (iAs, MMA, and DMA) is expressed as a percentage of the sum of each of these three species (). Arsenocholine (AsC) and arsenobetaine (AsB) are nontoxic forms of (organic) arsenic and were excluded from our analyses. We used DMA% as our primary measure of AME; MMA% and iAs% were used as secondary measures of AME. In the SHS and NHSCS, measures of AME were not computed if species were undetectable for an individual. If only one metabolite was undetectable, this missing metabolite was estimated as the LOD divided by the square root of 2 and was used in downstream analyses. For the present study, of the 2,719 HEALS participants, 523, 155, and 4 participants had values below the LOD for , , and MMA, respectively. We kept these participants in our analyses and set arsenic species values that were to zero. In SHS, we excluded 6 (of 997) participants with arsenic species with undetectable values and calculated imputed values for and MMA in 49 and 2 participants, respectively. In NHSCS, we excluded 39 (of 706) participants with arsenic species with undetectable values and calculate imputed values for , , and MMA in 225, 598, and 29 participants, respectively.

DNA Extraction and Targeted Sequencing

The details of sample collection and DNA extraction for HEALS and SHS participants have been described previously (Lee et al. 1990; Pierce et al. 2012). In brief, genomic DNA for HEALS samples were extracted from clotted blood using the Flexigene DNA Kit (catalog no. 51,204). For SHS participants, buffy coats were extracted from fasting blood samples using organic solvents and were stored at MedStar Health Research Institute (Lee et al. 1990). For NHSCS participants, frozen buffy coats retrieved for this study were thawed at room temperature and mixed by tube inversion and 5-s vortex. DNA was extracted using the Agencourt Genfind version 2 Kit.

For all DNA samples, we conducted targeted sequencing focusing on the exons of AS3MT. Sequencing was carried out on an Illumina MiSeq instrument using Illumina’s TruSeq Custom Amplicon Assay Kit. The custom kit was designed to sequence 781 small regions of in length (which includes additional content that is not relevant to the aims of the present study, including common variation in the 10q24.32 region).

Read Alignment, Variant Calling, and Quality Control

Targeted sequencing data was processed at the University of Chicago Bioinformatics Core using the GenomeAnalysisToolkit (GATK) (Poplin et al. 2018) and the GATK Best Practices Workflow (https://gatk.broadinstitute.org/hc/en-us/articles/360035535932-Germline-short-variant-discovery-SNPs-Indels-). For all samples, raw paired-end reads were aligned to the reference genome (hg19) using the Novoalign software. HaplotypeCaller [in genomic variant call format (GVCF) mode] was used to call variants [bi-allelic and insertion–deletions (indels)] per sample, generating intermediate GVCF files. Next, we used the GenomicsDBImport tool to consolidate the contents of GVCF files across samples into a single GVCF file, which was then used for joint calling across all samples, producing a set of joint-called bi-allelic variants and indels in a single VCF. We applied two sets of filters across all cohorts to the 414 raw variants and indels detected, using more stringent thresholds for indels because they are more prone to heterozygous false positives than bi-allelic variants. Variants matching any one of the following quality metric conditions were removed: QualbyDepth (), FisherStrand (), RMSMappingQuality (), StrandOddsRatio (), MappingQualityRankSumTest (), or . We removed indels that met at least one of the following conditions: , , , , or . The first filter excluded 141 bi-allelic variants and 17 indels, resulting in 256 variants passing all the quality metrics. In addition, we excluded all variants with call rates of , resulting 135 for HEALS, 92 for SHS, and 120 for NHSCS. Overall, we were able to sequence 9 of the 11 exons in AS3MT (Illumina was unable to design flanking probes for Exons 5 and 10) at a depth of coverage ranging between 173 and 729 per exon, across all cohorts (Figure 1).

Figure 1.

Average aligned read depth (depth of coverage) of AS3MT exons. Nine of 11 exons (excluding Exons 5 and 10) in AS3MT were sequenced with an average depth of coverage across all cohorts ranging from ( is deemed high quality). Rare, protein-altering variants were identified across all cohorts: five variants in HEALS in Exons 4 and 6, four variants in SHS in Exons 4, 8, 9, and 11, and five in NHSCS in Exons 6, 7, 9, and 11. Note: CCDS, Consensus Coding Sequence project; chr, chromosome; Human Feb. 2009, February 2009 Human reference sequence; HEALS, Health Effects of Arsenic Longitudinal Study; NHSCS, New Hampshire Skin Cancer Study; Rfam, RNA families; SHS, Strong Heart Study; UCSC, University of California, Santa Cruz.

Among the 2,719 HEALS participants selected for this study, we excluded 260 participants because of a low depth of coverage () and 25 participants who did not have genome-wide SNP data available for the generation of a kinship matrix, resulting in 2,434 remaining HEALS participants. Among the 997 SHS participants, we excluded 123 participants whose samples had extremely low coverage due to a very low number of reads and 6 participants with missing metabolite data, leaving 868 SHS participants. Among the 705 NHSCS participants selected for this study, we removed 39 samples owing to missing metabolite data and one sample due to a low number of reads, resulting in 666 NHSCS participants.

In an effort to validate our sequenced variants, we used existing HEALS SNP array and exome array data (Pierce et al. 2019) for 2,363 of the HEALS participants sequenced for this study. We identified 721 overlapping variants in both data sets (sequencing data and imputed SNP array data). We compared genotypes for each variant and found that 645 of 721 variants had consistent genotypes across 2,363 individuals. Only 3 variants had consistency. In addition, we used the Genome Aggregation Database (gnomAD) version 2.1.1 (Karczewski et al. 2020) to validate observed variants and to obtain an estimate of the rare variants missed in our study (specifically in Exons 5 and 10). GnomAD provides exome sequencing data for unrelated individuals, most of which are classified into six ancestral populations (African American, Latino, East Asian, Finnish, non-Finish European, and South Asian).

Variant Annotation

All variants were annotated using the Annotate Variation (ANNOVAR) software (Wang et al. 2010), which provided information on the impact of each exonic variant on the amino acid sequence (e.g., nonsynonymous, synonymous, frameshift indel). A variant was annotated as a predicted loss of function variant (pLoF) if they were categorized by ANNOVAR as frameshift, premature stop codon, loss of start or stop codon, or splice acceptor or donor. For this analysis, we focused on variants with a that altered the protein sequence. Thus, we excluded all synonymous, intronic, and intergenic variants detected. These exclusions resulted in five, four, and five AS3MT rare, protein-altering variants in the HEALS, SHS, and NHSCS cohorts, respectively.

We used the sorting intolerant from tolerant (SIFT) (Ng and Henikoff 2003) and PolyPhen (Adzhubei et al. 2010) tools to predict how variants coding for amino acid changes impact protein function. SIFT uses information on sequence homology and physical properties of amino acids to predict the impact on protein function, whereas PolyPhen uses sequence, phylogenetic, and structural information to empirically predict the effect of the substitution. We also report Combined Annotation Dependent Depletion (CADD) scores (Rentzsch et al. 2019) for the variants we analyzed, which are based on evolutionary information, conversion metrics, functional genomic data, and transcription information. CADD scores are Phred-like C-scores that rank variants relative to all possible substitutions of the human genome (). The CADD score represents the potential level of deleteriousness. For instance, variants with CADD scores of 0–10 are in the top 10% most deleterious, CADD scores 10–20 are in the top 1%, CADD scores 20–30 are in the top 0.1%, and so on.

Identification of Common Variants in 10q24.32 Region Associated with DMA%

In light of the well-established associations between common variants in the 10q24.32 region and AME, we used our sequencing data to genotype common variants in the 10q24.32 region and identify variants showing independent associations with DMA% in each population. We performed linear conditional, forward stepwise regression tests, stratified by cohort. We restricted analyses to 4,990 variants in the 10q24.32 target region and excluded variants with a Hardy-Weinberg () and variants with a (). This quality control resulted in 383 variants across all cohorts (HEALS , SHS , and NHSCC ).

Linear conditional, forward stepwise regression tests consisted of a series of association analyses. To identify the primary association signal, we individually tested each of the 383 common variants for an association with DMA%, controlling for age, sex, and population structure in HEALS and SHS. We used PLINK (Purcell et al. 2007) for SHS and NHSCS and genome-wide complex trait analysis (GCTA) (Yang et al. 2011) for HEALS (for adjustment of kinship matrix). To identify secondary independent association signals, we included the primary signal (identified in the previous analyses) as a covariate and repeated our association analyses. We repeated these analyses, adjusting for both the primary and secondary signals, to identify any remaining independent signals in each population. These analyses resulted in the identification of two independent lead variants in HEALS, rs12573221 () and rs14553735 (), three independent lead variants in SHS, rs10786722 (), rs4919687 (), and rs10883846 (), and a single lead variant in NHSCS, rs76255497 (). We conducted analyses of MMA% and iAs% and found that the DMA% associated lead variants were consistently among the top five lead independent variants in analyses of MMA% and iAs%, across all cohorts. Given this agreement in results across metabolites, we adjusted for the DMA% associated individual SNP alleles (0, 1, or 2) in our linear regressions. This adjustment ensured association estimates were not biased due to linkage disequilibrium (LD) between common variants known to impact AME and the rare variants studied in this work.

Statistical Analysis

The classical single-variant–based association test is not typically extended to analyses of rare variants because of a lack of statistical power (Asimit and Zeggini 2010). Instead, investigators have developed methods that aggregate variants in a biologically relevant region (i.e., a gene) and evaluate their cumulative effects. This approach is now the standard for rare variant studies and can provide reasonable power to detect association between a gene-based set of rare variants and a human trait (Lee et al. 2014). In the present study, the primary hypothesis was that carrying a rare, protein-altering variant in AS3MT reduces AME (represented by DMA%). To test this hypothesis, we conducted a burden test, where we assigned each individual a binary carrier status and tested whether the mean of DMA% differs between carriers and noncarriers (using linear regression). Individuals carrying at least one rare, protein-altering variant in AS3MT were assigned a carrier status of 1, and 0 otherwise. We did not observe any individuals carrying more than one rare, protein-altering variant in AS3MT. The burden test makes the strong assumption that all variants analyzed in the gene are causal and affect the trait in the same direction and with the same magnitude. Violation of these assumptions may result in loss of power (Neale et al. 2011).

We fit linear regression models for SHS and NHSCS, adjusted for age and sex, with carrier status as the predictor and arsenic metabolites (DMA%, MMA%, and iAs%) as outcomes. For HEALS, we conducted the burden test using a mixed-linear model based association test (–mlma) and incorporated a kinship matrix as implemented in the GCTA software (Yang et al. 2011) to control for cryptic relatedness among participants. For SHS, we accounted for population structure by adjusting for the first five principle components (PCs) derived from existing genome-wide SNP data described previously (PCs provided by SHS) (Matise et al. 2011). We were not able to control for population structure in NHSCS because of the lack of existing genome-wide SNP data; however, 99.5% of NHSCS participants included in this study self-reported as non-Hispanic white (only three individuals self-reported as Hispanic or Latino). To ensure minimal bias due to population structure, we conducted the same burden test excluding the three individuals who self-reported as Hispanic or Latino. In addition, we performed the burden test with adjustment and without adjustment for common SNPs in the 10q24.32 region.

To address the possibility of bias in the association between rare variant carrier status and AME due to prevalent skin lesion and SCC case status, we repeated the burden tests excluding 293 prevalent skin lesion cases in HEALS (1 skin lesion case was a carrier of the rs3523887 rare variant) and 349 SCC cases in NHSCS. We analyzed 2,076 HEALS individuals without prevalent skin lesions and 357 NHSCS controls and adjusted for age, sex, relatedness (in HEALS only), and common SNPs in the 10q24.32 region. To estimate the association between carrier status and DMA% across all cohorts, we meta-analyzed the association estimates for each cohort using the metafor package in R (version 1.3, R Development Core Team). We chose the meta-analyses approach instead of a pooled analysis because each analysis required a unique adjustment for population structure (i.e., in HEALS we adjusted for a kinship matrix, and in SHS we adjusted for genotype principle components). Tests of heterogeneity were used to determine the appropriate selection of a fixed-effects vs. a random-effects models. We used a random-effects model for a Cochran’s -statistic with a and a fixed-effects model otherwise. We conducted similar analyses using MMA% and iAs% as outcomes to ensure consistency with the analyses of DMA%.

In secondary analyses, we estimated the association between the collective effect of rare, protein-altering variants in AS3MT on DMA% using the sequence kernel association test (SKAT) (Wu et al. 2011), within each cohort. The SKAT is a score-based variance-component test that allows the variants within a gene to have different directions and magnitudes of effect and the SKAT statistic is primarily a representation of the variance of effect sizes within the set of variants. The SKAT software does not generate an effect estimate for the individual variants. Under the null hypothesis, none of the variants in the gene set have an effect on the arsenic metabolites. Conversely, under the alternative hypothesis, at least one variant in the gene set must have an , and a larger variance in phenotypic variation explained by the set of genetic variants would result in a larger SKAT statistic and a smaller -value. The SKAT software restricts all analyses to variants with missingness.

In SKAT analyses, we made adjustments for age, sex, population structure, and dosage scores of common SNPs associated with DMA% in the 10q24.32 region. HEALS was analyzed using the family-based SKAT (famSKAT) (Wang et al. 2017), which allows adjustment for a kinship matrix. We used both burden and nonburden SKAT gene-based tests because it was unknown which method would provide superior power and whether the assumptions of the burden test were reasonable. The meta-analysis version of SKAT, MetaSKAT (Lee et al. 2013), was used to derive summary SKAT statistics of rare protein-altering variants across the three cohorts. The MetaSKAT framework is currently unable to implement kinship matrices and, thus, HEALS estimates in the meta-analyses were not adjusted for relatedness. To assess any bias in the SKAT meta-analyses due to population structure in HEALS, we conducted sensitivity analyses of the burden and SKAT test where we excluded related individuals (kinship coefficient ) resulting in 1,285 unrelated individuals.

Last, we attempted to extrapolate our results to approximate the carriers’ increased risk of developing arsenic-related toxicities. We used previously reported odds ratios (ORs) for the risk of lung and bladder cancer in Chile (Melak et al. 2014) for a 1-unit percentage increase in MMA% (for lung and bladder cancer). We then converted the reported ORs to a logORs and multiplied them by the estimated meta-analysis effect of rare variant carrier status on MMA%. These values were then exponentiated to obtain ORs for the risk of developing lung and bladder cancer for rare variant carriers. In addition, we used a SNP-based Mendelian randomization (MR) approach to estimate carriers’ increased risk of skin lesions, using data from Bangladesh. We first estimated the impact of DMA% on skin lesion risk, using previously reported association estimates for DMA%-associated AS3MT SNPs (rs9527 and rs11191527) and FTCD SNP (rs61735836) (Pierce et al. 2013, 2019). We used a likelihood-based MR method that combines the summary statistics of multiple genetic variants into a single causal estimate for DMA% (Burgess et al. 2013). We multiplied this estimate by the association of carrier status with DMA%. Finally, we exponentiated this value to obtain the OR describing the risk of developing skin lesions given a difference in DMA% due to rare variant carrier status.

Results

Characteristics of the 2,434 HEALS, 868 SHS, and 666 NHSCS participants included in our analyses are described in Table 1. Total urinary arsenic (in micrograms per liter), represented by the sum of urinary concentrations of , , MMA, and DMA, varied across cohorts (Table 1), with substantially higher total urinary arsenic in HEALS () compared with SHS () and NHSCS (). Urinary DMA%, on average, was the lowest in HEALS (71.2%), higher in SHS (76.9%), and highest in NHSCS (81.1%) (Table 1). The distribution of arsenic metabolites by cohort are shown in Figure S1. The difference in metabolite percentages across cohorts is likely driven by the difference in overall level of iAs exposure. HEALS participants exhibited the lowest percentage of urinary DMA and experienced the highest level of iAs exposure. This is consistent with the reported inverse association between DMA% and iAs exposure, suggesting that the methylation machinery can become saturated at high doses of arsenic (Ahsan et al. 2007).

Table 1

Participant characteristics stratified by cohort.

Characteristics	HEALS (n=2,434)	SHS (n=868)	NHSCS (n=666)
Sex [n (%), male]	1,433 (58.9)	348 (40)	395 (59.3)
Age [y (mean±SD)]	41.7±10.5	56±8.2	64.4±7.6
BMI {kg/m^2 [n (%)]}
<18.5 (underweight)	1,043 (43.3)	7 (0.8)	7 (1)
18.5–24.9 (normal)	1,212 (50.3)	143 (16.6)	194 (29.6)
25–29.9 (overweight)	137 (5.7)	300 (34.7)	246 (37.5)
>30 (obese)	17 (0.7)	414 (47.9)	209 (31.9)
Total urinary arsenic [μg/L (mean±SD)]	137.4±162.8	14.8±16.9	7.3±8.9
25th	38.2	5.6	3.3
50th	82.6	9.8	5.0
75th	175.6	16.9	8.4
Min, max	3.7, 1,528.9	0.4, 161.9	0.7, 111.6
iAs% (mean±SD)	15.2±6.5	8.5±5.2	7.8±5.2
25th	11.0	5.2	4.3
50th	14.1	7.5	6.8
75th	18.3	10.7	9.8
Min, max	0, 70.3	0.6, 59.7	0.5, 35.7
MMA% (mean±SD)	13.7±5.2	14.5±5.5	10.4±4.5
25th	9.8	10.7	7.3
50th	13.1	13.9	10.1
75th	16.7	17.7	13.2
Min, max	0, 34.7	0.5, 45.9	0.6, 36.7
DMA% (mean±SD)	71.2±8.7	76.9±8.8	81.8±8.1
25th	66.1	72.1	77.0
50th	71.8	77.8	82.5
75th	77.2	83.3	87.2
Min, max	27.4, 92.9	32.4, 94.5	35.1, 97.9

Note: BMI, body mass index; DMA%, percentage of dimethylarsinic acid; HEALS, Health Effects of Arsenic Longitudinal Study; iAs%, percentage of inorganic arsenic; max, maximum; min, minimum; MMA%, percentage of monomethylated arsenic; NHSCS, New Hampshire Skin Cancer Study; SD, standard deviation; SHS, Strong Heart Study.

We identified 13 rare, protein-altering variants in AS3MT across all cohorts (Table 2). All but one of these variants (rs35232887 in HEALS and NHSCS) was specific to a single cohort, with five, four, and five variant sites detected in HEALS, SHS, and NHSCS, respectively. These variants were primarily observed once or twice, with 23 total individuals carrying a rare allele at one of these sites (13 carriers in HEALS, 5 in SHS, and 5 in NHSCS). Individual-level percentage of arsenic metabolites for rare-variant carriers are shown in Table S1. All five variants detected in HEALS were also present in gnomAD, and four of five variants detected in NHSCS were present in gnomAD (Table S2). Only one of the four variants detected in SHS was present in gnomAD, which may be due to the fact that gnomAD does not include data from individuals of Native American ancestry. Overall, 0.5% of HEALS and SHS participants were carriers of an AS3MT rare, protein-altering allele, and 0.7% were carriers in NHSCS. This is similar to, but slightly below the 0.9% of carriers of rare, protein-altering AS3MT variants in gnomAD (Karczewski et al. 2019) (variants in Exons 5 and 10 were excluded from the percentage of gnomAD carriers for consistency purposes).

Table 2

Rare variants observed in AS3MT across cohorts.

Study/rs ID	Position (BP)	Major allele/minor allele	Amino acid change	Number of carriers	Exon	Mutation	PolyPhena	SIFTb	CADD scorec	Carrier(s) DMA% distribution (decile)d
HEALSe
rs144713814	104632301	C/A	Ser/Arg	1	4	Missense	Probably damaging	Deleterious	18.9	2nd
rs759484385	104632326	G/A	Val/Met	2	4	Missense	Probably damaging	Deleterious	22	1st
rs563943815	104634388	A/C	Gln/His	7	6	Missense	Possibly damaging	Deleterious	14.4	5th
rs35232887h	104634407	C/T	Arg/Trp	2	6	Missense	Probably damaging	Deleterious	21.8	1st
rs528680133	104634419	G/A	NA	1	6	Splice donor	—	—	26	1st
SHSf
NA	104632249	G/A	Cys/Tyr	1	4	Missense	Possibly damaging	Deleterious	27.6	2nd
NA	104638178	T/TC	Leu/—	1	8	Frameshift insertion	—	—	—	2nd
rs139656545	104638615	G/A	Arg/His	1	9	Missense	Benign	Deleterious	23.8	4th
NA	104660404	A/T	Arg/—	2	11	Stop gained	—	—	35	1st
NHSCSg
rs770395949	104634377	G/GAT	Lys/—	1	6	Frameshift insertion	Probably damaging	Deleterious	24.7	1st
rs35232887h	104634407	C/T	Arg/Trp	1	6	Missense	Probably damaging	Deleterious	21.8	6th
rs775931783	104636754	A/T	Glu/Val	1	7	Missense	Benign	Tolerated	22.8	5th
NA	104638627	C/T	Ala/Val	1	9	Missense	Possibly damaging	Deleterious	27.0	2nd
rs182365639	104660417	A/T	Asp/Val	Singleton	11	Missense	Benign	Tolerated	11.5	2nd

Note: —, not applicable; A, adenine; Ala, alanine; Arg, arginine; Asp, asparagine; BP, base pair; C, cytosine; CADD, combined annotation dependent depletion; Cys, cysteine; DMA%, percentage of dimethylarsinic acid; G, guanine; Gln, glutamine; Glu, glutamic acid; HEALS, Health Effects of Arsenic Longitudinal Study; His, histidine; Leu, leucine; Lys, lysine; Met, methionine; NA, not available; NHSCS, New Hampshire Skin Cancer Study; rs ID, reference single nucleotide polymorphism cluster dentification number; Ser, serine; SHS, Strong Heart Study; SIFT, sorting intolerant from tolerant; T, thymine; Trp, tryptophane; Tyr, tyrosine; Val, valine.

PolyPhen predicts the impact of amino acid changes.

SIFT predicts the impact of amino acid changes.

CADD scores of 0 to 10 are in the top 10% most deleterious, scores 10 to 20 are in the top 1%, CADD scores 20 to 30 are in the top 0.1%, and so on.

This value is based on the distribution of DMA% across carriers and noncarriers and corresponds to the average decile when a variant has carrier.

The proportion of rare variant carriers in HEALS was 13 of 2,434 (0.5%).

The proportion of rare variant carriers in SHS was 5 of 865 (0.5%).

The proportion of rare variant carriers in NHSCS was 5 of 666 (0.7%).

Only one variant (rs35232887) was observed in multiple cohorts.

Most variants were missense changes, but each cohort had at least one pLoF variant (i.e., stop gain, frameshift, or splice donor). In HEALS and NHSCS, the two pLoF variants were each observed once, rs528680133 (splice donor) and 10:104,634,377 (frameshift), with high CADD scores of 26 and 24.7, respectively. In SHS, we identified two pLoF variants, 10:104,660,404 coding for a stop-gain with a CADD score of 35 (observed in two individuals) and 10:104,638,178, coding for a frameshift insertion (observed in one individual).

Carrier status for rare, protein-altering variants in AS3MT was associated with lower mean DMA% on average compared with noncarriers, across all cohorts in analyses unadjusted for common variants in the 10q24.32 region [HEALS: (95% CI: , ), ; SHS: (95% CI: , ), ; NHSCS: (95% CI: , ), ]. These results were consistent in adjusted analyses [HEALS: (95% CI: , ), ; SHS: (95% CI: , ), ; NHSCS: (95% CI: , ), ] (Table 3 and Figure 2). In SHS, at least one rare variant carrier had missing genotypes for the common variants rs10786722 and rs4919687, so we instead adjusted for dosage of proxy SNPs, 10:104,723,620 and 10:104,685,493, respectively, with no missing genotypes among carriers. In addition, sensitivity analyses in NHSCS where we excluded three individuals who self-reported as Hispanic or Latino showed the same results (, ) and thus confirming that the inclusion of these individuals did not bias our results. Carrier status showed evidence of a positive association with MMA% and iAs% (Table 3) across all cohorts, but the association exceeded the threshold for SHS in analyses unadjusted and adjusted for common variants.

Table 3

Association between carrier status of rare, protein-altering AS3MT variants and arsenic metabolism phenotypes.

Cohort	N	Carriers (n)	DMA%		MMA%		iAs%
			β (95% CI)	p-Value	β (95% CI)	p-Value	β (95% CI)	p-Value
Unadjusteda
HEALSc	2,369	13	−9.7 (−14.4, −5.11)	3.9×10^−5	4.9 (2.2, 7.5)	0.0003	4.8 (1.3, 8.3)	0.007
SHSd	865	5	−10.2 (−17.6, −2.9)	0.006	3.9 (−0.8, 8.5)	0.11	6.4 (2.0, 10.7)	0.004
NHSCS	666	5	−9.3 (−16.3, −2.4)	0.009	5.8 (1.9, 9.8)	0.004	3.5 (−0.9, 8.0)	0.12
Meta-analysise	3,900	23	−9.8 (−13.2, −6.3)	2.2×10^−8	4.9 (2.9, 6.9)	1.3×10^−6	4.9 (2.6, 7.2)	3.6×10^−5
Adjustedb
HEALSc	2,369	13	−9.4 (−13.9, −4.8)	5.2×10^−5	4.7 (2.1, 7.4)	0.0005	4.6 (1.2, 8.0)	0.008
SHSd	865	5	−6.87 (−13.5, −0.21)	0.04	2.0 (−2.39, 6.45)	0.37	4.8 (0.63, 9.05)	0.02
NHSCS	666	5	−8.9 (−15.6, −2.16)	0.01	5.6 (1.76, 9.44)	0.004	3.3 (−1.09, 7.69)	0.15
Meta-analysise	3,900	23	−8.7 (−11.9, −5.4)	1.9×10^−7	4.5 (2.5, 6.4)	9.7×10^−6	4.3 (2.0, 6.6)	0.0002

Note: CI, confidence interval; DMA%, percentage of dimethylarsinic acid; GCTA, genome-wide complex trait analysis; HEALS, Health Effects of Arsenic Longitudinal Study; iAs%, percentage of inorganic arsenic; MMA%, percentage of monomethylarsonic acid; NHSCS, New Hampshire Skin Cancer Study; SHS, Strong Heart Study; SIFT, sorting intolerant from tolerant; SNP, single nucleotide polymorphism.

Models adjusted for sex and age (continuous).

Models adjusted for sex, age (continuous), and common variants in the 10q24.32 region (variant dosage score, e.g., 0, 1, or 2) (HEALS: rs12573221 and 145537350; SHS: 10:104723620, rs10883846; and 10:104685493, and NHSCS: rs76255497).

Burden test was conducted on a subset of 2,369 individuals with genome-wide SNP data available to adjust for kinship using the mixed-effects model in GCTA.

Included additional adjustment for first five principal components derived from genome-wide SNP data.

Performed using fixed-effect models (test of heterogeneity -values across metabolites ).

Figure 2.

Percentage of arsenic metabolites by carrier status of rare, protein-altering variants in AS3MT across cohorts. Measures of arsenic metabolites within each box range from the 25th to the 75th percentile. The 50th percentile (or median) is depicted by the horizontal line within each box. Data in the upper and lower whiskers of the box represent measures below and above the 25th and 75th percentiles, respectively. Individual data points for all noncarriers are shown to the left of the box plot and to the right of the box plot for carriers. There were a total of 13, 5, and 5 rare variant carriers in HEALS, NHSCS, and SHS, respectively. (A) Carriers of rare, protein-altering variants have reduced DMA% (on y-axis) compared with noncarriers, across all cohorts. Carriers of loss of function (pLoF) variants, labeled accordingly by cohort and shown by arrows, fall in the bottom 10th and 20th percentiles of DMA%. (B) MMA% (on y-axis) is higher among carriers of rare, protein-altering variants across all cohorts. (C) iAs% (on y-axis) is higher among carriers of rare, protein-altering variants across all cohorts. Note: DMA%, percentage of dimethylarsinic acid; HEALS, Health Effects of Arsenic Longitudinal Study; MMA%, percentage of monomethylated arsenic; NHSCS, New Hampshire Skin Cancer Study; SHS, Strong Heart Study.

We conducted meta-analyses across the three cohorts (3,900 participants) for unadjusted and adjusted models (for common variants in 10q24.32) across all metabolites. Tests of heterogeneity indicated fixed-effects models () were appropriate for all metabolites in adjusted and unadjusted models (unadjusted: , , , , and adjusted: , , ). These meta-analyses indicated that, relative to noncarriers, being a carrier of a rare protein-altering variant in AS3MT was associated with 9.8% lower mean DMA% [ (95% CI: , )], 4.9% higher mean MMA% [ (95% CI: 2.9, 6.9)], and 4.9% higher mean iAs% [ (95% CI: 2.6, 7.2)] with no adjustment for common variants in 10q24.32. Meta-analyses adjusted for common variants showed similar results, where carriers had 8.7% [ (95% CI: , )] lower mean DMA%, 4.5% [ (95% CI: 2.5, 6.4)] higher mean MMA%, and 4.3% [ (95% CI: 2.0, 6.6)] higher mean iAs%.

Sensitivity analyses in HEALS and NHSCS restricted to individuals without prevalent skin lesions or SCC showed a small attenuation in the association of rare variant carrier status with DMA% in 2,076 HEALS participants [ (95% CI: , )] and an increased strength of association in 357 NHSCS controls [ (95% CI: , )] (Table S3) compared with analyses of all individuals (Table 3). These results remained consistent with the results prior to case status exclusion.

Using SKAT, we observed for the association between the collective effects of rare, protein-altering variants in AS3MT and DMA% in SHS () and NHSCS (), but not in HEALS () (Table 4). Meta-analyses across cohorts using SKAT provided evidence of association across all cohorts (). Analyses using MMA% and iAs% as outcomes in SKAT produced results similar to the burden test (Table 4).

Table 4

Association between rare, protein-altering variants in AS3MT and arsenic metabolism phenotypes across cohorts using a non-burden (i.e., SKAT) testing method.

Cohort	N	Carriers (n)	DMA%	MMA%	iAs%
			p-Value	p-Value	p-Value
HEALSa	2,369	13	0.11	0.09	0.09
SHS	865	5	0.03	0.27	4.3×10^−6
NHSCS	666	5	0.005	0.08	0.003
Meta-analysisb	3,900	23	0.002	0.1	1.4×10^−7

Note: SKAT is a score-based variance component test; under the null hypothesis, all rare variants in the gene set have an . DMA%, percentage of dimethylarsinic acid; HEALS, Health Effects of Arsenic Longitudinal Study; MMA%, percentage of monomethylarsonic acid; iAs%, percentage of inorganic arsenic; NHSCS, New Hampshire Skin Cancer Study; SHS, Strong Heart Study; SIFT, sorting intolerant from tolerant; SKAT, sequence kernel association test.

SKAT analyses were implemented using the linear mixed model (EMMAX) built-in to the software to adjust for relatedness in the HEALS cohort.

MetaSKAT software does not allow for adjust of a kinship matrix, thus SKAT meta-analyses were not adjusted for relatedness in the HEALS cohort.

It is important to note that the meta-analyses performed in the MetaSKAT software required individual-level data and did not allow adjustment for a kinship matrix in the HEALS cohort. Thus, we performed burden and SKAT analyses of 1,285 unrelated individuals in the HEALS cohort across all metabolites without kinship adjustment. These analyses also excluded two carriers of one rare, protein altering variant in AS3MT (rs563943815). Overall, we found that excluding related individuals in HEALS slightly increased the effect of carrier status on arsenic metabolites in the burden test and decreased the -value in SKAT (Table S4). These results also highlight the possibility that not all nonsynonymous variants included in our rare variant test will impact protein function. These results suggest that accounting for kinship matrix in the meta-analyses has little impact on the estimated effect sizes from the burden test and the -values from SKAT.

Using exome sequencing data from gnomAD, we assessed the extent of missed rare, protein-altering variants in Exons 5 and 10 of AS3MT, which were not captured in the present study. In Exon 5, among all participants, there were 19 rare (), protein-altering variants (15 missense and 3 pLoF) (Table S5). In Exon 10, there were 14 rare, protein-altering variants (12 missense and 2 pLoF). The carrier frequency in gnomAD was 0.09% (52 of 60,706) for Exon 5 and 0.4% (242 of 60,706) for Exon 10. Across all populations and all AS3MT exons in gnomAD, the percentage of carriers of rare, protein-altering variants in AS3MT was 1.4% (859 of 60,706) (Table S6), approximately double what we observed in our three cohorts (0.5–0.7%).

We used results from prior studies to estimate the association of carrier status with disease outcomes (lung cancer, bladder cancer, and skin lesions) in exposed populations. A Chilean study (Melak et al. 2014) of lung and bladder cancer report logORs of 0.10 and 0.04 for a 1-unit increase in MMA%, respectively. These logORs and a 4.5% increase in MMA% among carriers (from meta-analysis, Table 3) correspond to a 59% () and 19% () estimated increase risk of developing lung and bladder cancer compared with noncarriers of AS3MT rare variants, respectively (Table S7). For skin lesions, we used an MR approach to estimate the association of DMA% with skin lesion outcomes. We used previously reported association estimates and SEs (from the HEALS cohort) for the a) of AS3MT SNPs and FTCD SNP with DMA%; and b) associations of AS3MT SNPs and FTCD SNP with skin lesions (Table S8) (Pierce et al. 2013, 2019). Given an estimated logOR of 0.069 for a 1-unit decrease in DMA% and an 8.7% reduction in DMA% among carriers (from meta-analysis, Table 3), the estimated the risk of developing skin lesions for carriers in the HEALS cohort was 82% higher () compared with noncarriers of AS3MT rare variants (Table S9).

Discussion

In this study, we conducted targeted sequencing of the coding regions of the AS3MT gene across three arsenic-exposed cohorts of different genetic ancestry. Our results suggest that on average, rare, protein-altering variants in AS3MT are associated with decreased AME (lower urinary DMA%), and these associations were independent of the effects of known common variants in this region. Our burden test results provide evidence that rare variants in AS3MT, which are predicted to affect protein structure and function, may reduce the efficiency of arsenic metabolism. These variants likely result in increased retention of arsenic in the body [as seen in AS3MT KO mice (Drobna et al. 2009)] and increased risk of arsenic-related toxicities and health effects (as seen in humans carrying common alleles associated with AME). This is the same pattern of association observed among the HEALS study participants by Pierce et al. (2012, 2013, 2019) for common variants that are associated with AME, namely, the risk allele decreases DMA% and increases MMA% and iAs%. This study represents one of the first attempts to assess the effects of rare inherited variation on AME in humans, and we address this question in multiple population groups of different ancestries and arsenic exposure levels.

The impact of common variation in the AS3MT (10q24.32) region on AME is well established. Epidemiologic studies report consistent and reproducible associations between common SNPs in and around AS3MT and arsenic metabolites in urine across multiple populations (Agusa et al. 2009; Balakrishnan et al. 2017; Engström et al. 2011, 2015; García-Alvarado et al. 2018; Pierce et al. 2012). The association of these variants with risk for arsenic-induced skin lesions (Pierce et al. 2013) highlights the relevance of these metabolism-related variants to arsenic toxicity risks. At least one of the causal variants in this region appears to have effects that are regulatory in nature, impacting expression of local genes, including AS3MT (Chernoff et al. 2020; Engström et al. 2013). Other than AS3MT, only one other gene, the gene for formimidoyltransferase cyclodeaminase (FTCD), has been shown to contain inherited genetic variation that affects AME. In a Bangladeshi population, the minor allele (A) of a missense variant (rs61735836) in Exon 3 of FTCD was shown to be associated with decreased DMA% and increased skin lesion risk in Bangladeshi individuals (Pierce et al. 2019). Common SNPs in AS3MT and FTCD account for of the variation in DMA% (Pierce et al. 2019). The FTCD association had not been identified at the time targeted sequencing was conducted, so we could not examine rare variants in FTCD in the present study.

Prior heritability studies suggest a role for rare variation in AME. For example, Gao et al. (2015) estimated the full narrow-sense heritability () of AME to be 41% [95% CI: (7.7%, 74.3%), calculated based on reported SE of 17%], including the additive effects of both common and rare variants but excluding (adjusting for) the effects of common SNPs in the AS3MT region with known associations with AME. They also estimated the heritability due to common SNPs to be only after similar adjustments. The difference in these two estimates suggests that rare variants (and unmeasured shared-environmental factors) may make an important contribution to interindividual variation in AME. Similarly, in a family study conducted within SHS, the heritability for DMA% was as high as 63% (Tellez-Plaza et al. 2013). Together, these findings suggest that rare variation contributes to familial similarity in AME phenotypes.

AS3MT knock out (KO) mice can be used to understand the potential effects of AS3MT pLoF variants in humans. AS3MT KO studies have demonstrated the critical role of AS3MT in the elimination of arsenic from tissues, clearance of arsenic in urine, and prevention of arsenic toxicities. KO mice exhibit dramatically higher arsenic concentrations and higher proportions of iAs in numerous tissue types—including liver, bladder, and kidney—several hours post arsenic dosing compared with wild-type (WT) mice with a similar dose (Chen et al. 2011; Currier et al. 2016; Drobna et al. 2009; Hughes et al. 2010). Whole-body clearance of iAs is substantially slower in KO compared with WT mice. For instance, Drobna et al. (2009) showed that at 24 h post dosing, KO mice retained 50% of the dose, whereas WT mice retained only 6%. In addition, KO mice are at an increased risk for arsenic toxicities (Douillet et al. 2017; Negro Silva et al. 2017; Yokohira et al. 2010). The AS3MT KO mice models developed to date are homozygous for the deleted gene (). In contrast, the mutation carriers in this study are heterozygous for a AS3MT rare variant likely to be damaging. Thus, we do not expect the impaired arsenic metabolism in our participants to be as extreme as that observed in AS3MT KO mice. However, the pLoF variants we observe (i.e., rs528680133 splice donor in HEALS, 10:104660404 stop gain variant in SHS, and 10:104634377 frameshift in NHSCS) are likely to produce human phenotypes that most closely resemble those of the KO mice. In human studies, we are unable to obtain tissue-specific measures of arsenic and therefore cannot determine how the elimination of arsenic varies among tissue types. However, we did observe that carriers of pLOF variants exhibit strikingly low metabolism efficiency given that these individuals were consistently in the bottom first decile of the DMA% distribution (Table 2). The effect of missense variants on AME is not as striking, and this is expected because some of these amino acid substitutions may have small or negligible effects on protein function.

In the present study, we performed two types of gene-based tests to assess the association between rare AS3MT variants and AME phenotypes. The first was a burden test where we summed each individual’s rare allele count (equivalent to carrier status in our case where no individual carries rare protein-coding allele in AS3MT), and the second was a nonburden test (SKAT) where we tested the variance component explained by the genetic variants among the total phenotypic variations. The burden test showed that on average, rare variants in AS3MT were inversely associated with DMA% across all three populations. Analyses using SKAT also showed clear association between AS3MT variants and DMA%. Both the burden test and SKAT methods also showed clear associations for AS3MT variants with both iAs% and MMA%. Meta-analyses across all three cohorts suggested that on average, carriers of rare, protein-altering variants in AS3MT have lower DMA% compared with noncarriers, a substantially larger effect compared with any single common variant associated with DMA% (i.e., AS3MT or FTCD SNPs). Furthermore, we expect that this estimated association of rare, protein-altering AS3MT variants with AME may be attenuated as compared with the association of known deleterious variants with AME because some individuals may be carrying missense variants that have a negligible impact on AS3MT function.

Our study was limited by our lack of data on variants in Exons 5 and 10 of AS3MT, due to limitations of the Illumina TruSeq custom amplicon kit. We assessed the extent of missingness in our study using exome sequencing data from gnomAD, and we estimated that because of our lack of data in Exons 5 and 10 we were unable to capture of rare variant carriers. It is worth noting that Exon 10 contains a rare variant with 205 carriers in gnomAD (frequency of 0.02%) (Table S5). Thus, it is likely that sequencing of Exons 5 and 10 would identify a substantial number of additional variants not captured in the present study. Consequently, we expect the AS3MT protein-altering variant carrier proportion to be higher than what we report in this study.

Additional research is needed to further our understanding of the role of rare, protein-altering AS3MT variants in AME and in the etiology of arsenic-related health conditions. A broader representation of arsenic-exposed populations is necessary to replicate these findings, identify additional rare AS3MT variants, and assess their effects on the risk for arsenic toxicities. We are unable to examine the association of rare AS3MT variants with skin lesion status (in HEALS) and SCC (in NHSCS) directly because of the small number of cases who are carriers of rare variants (three cases of prevalent and incident skin lesions and no cases of SCC were carriers). However, using previously reported risk estimates for skin lesions (Pierce et al. 2013, 2019) and lung and bladder cancer (Melak et al. 2014), we estimated that the risk of developing skin lesions was 82% higher given an 8.7% reduction in DMA% and the risks of developing lung and bladder cancer were 59% and 19% higher, respectively, given a 4.5% increase in MMA%. A larger cohort study is needed to assess the impact of rare AS3MT variants on risk of arsenic toxicities and explore the interactions of rare-variant carrier status with sex, lifestyle factors, nutritional status, and overall iAs exposure. Additional experimental research is also needed to fully understand the effects of these variants on AS3MT protein structure and function.

In summary, we used data from multiple arsenic-exposed cohorts to assess the association between rare, protein-altering variants in AS3MT and arsenic metabolism efficiency. We provide evidence that rare variants in AS3MT decrease arsenic metabolism efficiency, a finding consistent across multiple populations with distinct ancestral backgrounds. Although we were unable to assess the effect of these rare, protein altering variants on arsenic-related disease outcomes (due to the small number of carriers), our findings suggest that at least some of the variants may affect the internal dose of arsenic (through effects on arsenic metabolism and clearance) and, in turn, influence the risk of arsenic-related toxicities. Our results, together with our knowledge of common AS3MT variants, highlight AS3MT as a major susceptibility locus for arsenic metabolism efficiency with implications for arsenic toxicities.

Click here for additional data file.

Click here for additional data file.