Potential role of regulatory DNA variants in modifying the risk of severe cutaneous reactions induced by aromatic anti-seizure medications.

OBJECTIVE: Stevens-Johnson syndrome (SJS) and toxic epidermal necrolysis (TEN) are severe cutaneous adverse drug reactions. Antiseizure medications (ASMs) with aromatic ring structure, including carbamazepine, are among the most common culprits. Screening for human leukocyte antigen (HLA) allele HLA-B*15:02 is recommended prior to initiating treatment with carbamazepine in Asians, but this allele has low positive predictive value.
METHODS: We performed whole genome sequencing and analyzed 6 199 696 common variants among 113 aromatic ASM-induced SJS/TEN cases and 84 tolerant controls of Han Chinese ethnicity.
RESULTS: In the primary analysis, nine variants reached genome-wide significance (p < 5e-08), one in the carbamazepine subanalysis (85 cases vs. 77 controls) and a further eight identified in HLA-B*15:02-negative subanalysis (35 cases and 53 controls). Interaction analysis between each novel variant from the primary analysis found that five increased risk irrespective of HLA-B*15:02 status or zygosity. HLA-B*15:02-positive individuals were found to have reduced risk if they also carried a chromosome 12 variant, chr12.9426934 (heterozygotes: relative risk = .71, p = .001; homozygotes: relative risk = .23, p < .001). All significant variants lie within intronic or intergenic regions with poorly understood functional consequence. In silico functional analysis of suggestive variants (p < 5e-6) identified through the primary and subanalyses (stratified by HLA-B*15:02 status and drug exposure) suggests that genetic variation within regulatory DNA may contribute to risk indirectly by disrupting the regulation of pathology-related genes. The genes implicated were specific either to the primary analysis (CD9), HLA-B*15:02 carriers (DOCK10), noncarriers (ABCA1), carbamazepine exposure (HLA-E), or phenytoin exposure (CD24). SIGNIFICANCE: We identified variants that could explain why some carriers of HLA-B*15:02 tolerate treatment, and why some noncarriers develop ASM-induced SJS/TEN. Additionally, this analysis suggests that the mixing of HLA-B*15:02 carrier status in previous studies might have masked variants contributing to susceptibility, and that inheritance of risk for ASM-induced SJS/TEN is complex, likely involving multiple risk variants.

We identified novel genome‐wide significant variants for ASM, carbamazepine, and HLA‐B*15:02‐negative cohorts

Interaction analysis highlighted why some HLA‐B*15:02 carriers were less likely to develop SJS/TEN

Identified suggestive variants from the various analyses were associated with regulatory regions

Suggestive variants were predicted to modify gene expression of CD9, DOCK10, HLA‐E, ABCA1, and CD24

INTRODUCTION

Epilepsy affects 7.6 per 1000 people worldwide. The mainstay of treatment is antiseizure medications (ASMs). Although effective in most patients, ASMs are associated with a range of adverse effects. Idiosyncratic reactions, most commonly cutaneous reactions, are among the most problematic adverse effects, because they are largely dose‐independent. Up to 10% of new users of ASMs develop some form of cutaneous adverse drug reactions (cADRs), typically between 1 and 4 weeks after commencement of ASMs, particularly for drugs with an aromatic ring structure. Although the majority of cADRs are mild self‐limiting maculopapular exanthema, a minority of individuals develop more severe reactions, including Stevens–Johnson syndrome (SJS) and toxic epidermal necrolysis (TEN). These reactions are characterized by blistering exanthema of macules and targetlike lesions accompanied by, to different extents, cutaneous and mucosal detachment. SJS/TEN has a mortality rate of 15%, and survivors have long‐term disabilities including vision loss and scarring. There is a pressing clinical need to identify high‐risk patients and to understand the pathophysiology of these reactions.

Asians have increased incidence of SJS/TEN after commencing aromatic ASMs (e.g., 1 in 400 in Han Chinese) compared to Caucasians (1 in 10 000). , In Asians, several common genetic variants of the human leukocyte antigen (HLA) have been identified to be associated with ASM‐induced SJS/TEN. The most reported is the HLA‐B*15:02 allele, which is strongly predictive of SJS/TEN induced by carbamazepine across broad Asian populations, and for which pretreatment screening has been recommended by the US Food and Drug Administration (FDA) since 2007. Due to this recommendation, fewer HLA‐B*15:02‐positive patients have received carbamazepine, reducing incidence of SJS/TEN arising from this drug. However, due to changes in prescribing habits of physicians, there has been an increased incidence of phenytoin‐SJS/TEN. Other risk alleles for carbamazepine‐induced SJS/TEN reported in Han Chinese include HLA‐A*24:02 and possibly HLA‐A*31:01 (also reported in European Caucasians). However, these HLA alleles are neither necessary nor specific and have low positive predictive values (<10%). In Asians, HLA‐B*15:02 has also been suggested to be associated with SJS/TEN induced by phenytoin and oxcarbazepine, although the genetic predisposition to SJS/TEN induced by ASMs other than carbamazepine is largely unknown.

We aimed to identify additional common risk variants for SJS/TEN induced by aromatic ASMs. We compared the frequencies of common variants in patients who developed ASM‐SJS/TEN (cases) and matched drug‐tolerant and population controls, followed by subanalysis of patients based on HLA‐B*15:02 carrier status and drug‐specific exposure. We modeled interaction between each significant variant identified and HLA‐B*15:02 status to determine whether carrying the variant modified the likelihood of developing SJS/TEN. Lastly, we investigated the potential functional impact of significant and suggestive variants.

MATERIALS AND METHODS

Cases and controls

Han Chinese patients with SJS/TEN induced by aromatic ASMs (carbamazepine, oxcarbazepine, lamotrigine, and phenytoin) were recruited as cases from southern China and Malaysia. As described in our previous studies, , the diagnosis of SJS/TEN was based on the criteria by Roujeau and Stern defined by skin detachment in two or more mucosal sites, and was confirmed by dermatologists. Patients who had been exposed to the matched ASMs for at least 3 months without any form of cADR were recruited as drug‐tolerant controls. Patients who had tolerated more than one ASM acted as drug‐matched controls for multiple cases. To identify potential variants that modify the effect of HLA‐B*15:02 on the risk of carbamazepine‐induced SJS/TEN, we preferentially recruited carbamazepine‐tolerant controls who were known to carry this allele. All cases and controls were of Han Chinese descent. A venous blood sample was collected from each patient for DNA extraction. Due to the 2007 FDA recommendation, few patients who harbor HLA‐B*15:02 have been prescribed carbamazepine or other aromatic ASMs. For this reason, recruitment of a replication cohort was not feasible.

The study was approved by the following institution and ethics committees: Joint Chinese University of Hong Kong‐New Territories East Cluster Clinical Research Ethics Committee (no. 2004.068), Ethics Committee of the Second Affiliated Hospital of Guangzhou Medical University (no. 18‐8071‐BO), University of Malaya Medical Center Medical Ethics Committee (no. 950.49), and the Ministry of Health Medical Research and Ethic Committee (no. nmrr‐13‐1157‐16170). All methods were performed in accordance with the relevant guidelines and regulations. Patients provided written informed consent and were subsequently deidentified.

Whole genome sequencing and bioinformatics pipeline

Whole genome sequencing was performed on the extracted DNA with an average of 30× coverage in batches using either the Illumina HiSeq X Ten platform (Kinghorn Centre for Clinical Genomics, Garvan Institute of Medical Research, Darlinghurst, NSW, Australia) or the BGI‐500 platform (BGI). The Picard/Genome Analysis Toolkit (GATK) data processing pipeline was used to align raw sequence reads to the human reference genome (build GRCh38/hg38) and to call variants (details provided in Supplementary Methods).

Quality control

Variants were restricted to autosomal regions and filtered to remove likely false positives according to the following criteria: (1) GATK variant quality score recalibration truth tranche scores < 99.95, (2) not in Hardy–Weinberg equilibrium, and (3) missing rate of >2%. Within each sequencing batch, alleles with a frequency of <3% were removed to mitigate batch effects, resulting in a common variant call set with minor allele frequency of >2.5% (Figure S1). Individuals identified as outliers by multidimensional scaling were removed together with those with a missing rate of >2% or a heterozygosity rate deviating by >3 SD from the mean (Table S1). Identity by descent testing found no cryptically related individuals, and principal component analysis using 2548 individuals from the 1000 Genomes resource confirmed East Asian descent for all cases and controls (Figure S2). Anonymized data will be shared by request from any qualified investigator under compliance with institutional policy.

Variant annotation

Variants were annotated using Ensembl's Variant Effect Predictor tool. The Ensembl database was also used to determine whether significant or suggestive variants fell within regulatory features—regions of the DNA that regulate gene expression. These features include promoter regions, enhancers, and transcription factor binding sites. Their location in the DNA is based on evidence from multiple sources including the NIH Roadmap Epigenomics Mapping Consortium, the Encyclopedia of DNA Elements (ENCODE) Consortium, and the BLUEPRINT analysis portal. Regulatory features are generated using genome‐wide assays including chromatin assays (DNase‐seq), histone modification assays (CHiP‐seq) and transcription factor binding assays. Expression quantitative trait loci (eQTL) and splicing quantitative trait loci (sQTL) data were obtained from the Genotype‐Tissue Expression (GTEx) portal, and information on regulatory elements was sourced from the Ensembl website (Release 103, February 2021, http://asia.ensembl.org/index.html).

Statistical analysis

PLINK version 1.9 was used to conduct association analyses, using the allelic model, and to calculate linkage disequilibrium (LD) between variants.

In the primary analysis, we compared allele frequencies between all ASM‐SJS/TEN cases and ASM‐tolerant controls. We performed subanalyses by stratifying cases and controls based on their HLA‐B*15:02 carrier status, and by comparing cases and tolerant controls exposed only to carbamazepine, or phenytoin. We did not perform subanalysis for oxcarbazepine or lamotrigine due to limited sample sizes. We used the gnomAD database, viewed in March 2021, and Fisher exact test to compare significant variant allele frequencies to frequencies observed in East Asian population controls. The interaction between genome‐wide significant variants and HLA‐B*15:02 status was investigated using generalized linear model with Poisson distribution, log‐link, and robust variance in Stata version 16 (StataCorp).

The standard genome‐wide p‐value cutoff for statistical significance of 5e‐08 was applied to identify variants that modified risk. To interrogate variant enrichment in subthreshold regions (suggestive variants) observed in the quantile–quantile plots, three criteria were applied: (1) a p‐value between 5e‐06 and 5e‐08, (2) five or more variants in LD within a region that harbors a regulatory feature, and (3) evidence from GTEx data that a variant alters gene expression. The three most significant findings for each variant are reported in the supplementary tables.

For the genome‐wide significant ASM variants, we used Fisher exact test to compare their frequencies between cases and population controls obtained from the gnomAD (version 3) East Asian genomes. Applying Bonferroni correction for multiple testing, p < .01 was considered statistically significant. Figures were generated using R Studio packages ggplot2 and qqman, and locusZoom version 1.4.

RESULTS

Patient characteristics

A total of 116 aromatic ASM‐SJS/TEN cases and 85 aromatic ASM‐tolerant controls were recruited and underwent whole genome sequencing. Data from four individuals (three cases and one control) failed quality control (Supplementary Methods), leaving 113 cases and 84 controls available for analysis (Table 1). In the majority of cases, SJS/TEN was induced by carbamazepine (85/113, 75.2%), followed by lamotrigine (14/113, 12.4%) and phenytoin (11/113, 9.7%). The power of the study was calculated using the online GCTA‐GREML Power Calculator (https://shiny.cnsgenomics.com/gctaPower/), which can estimate power relative to heritability based on our numbers of cases and controls, disease prevalence in the general population (1 in 1 000 000), and estimated type I error (.05). Due to the rarity of SJS/TEN in the population, lower samples sizes are needed to detect differences based on the parameters of the calculator. For this sample size, ~80% power can detect variants with heritability of ~75%. Therefore, this dataset had the capacity to detect variants despite its small sample size (Figure S3).

TABLE 1

Characteristics of cases and controls

Characteristic	Cases, n = 113	Tolerant, n = 84
Female sex	56	42
Drugs
Carbamazepine	85	77
Lamotrigine	14	15
Phenytoin	11	40
Oxcarbazepine	3	4
HLA‐B*15:02
Carrier	78	31
Noncarrier	35	53

Through retrospective screening, we identified 69% (78/113) of cases as HLA‐B*15:02+ and 40% (31/84) of controls to be HLA‐B*15:02+. The Allele Frequency Net Database identified 9.38% (712/7595) of people from the Hong Kong Chinese Bone Marrow Registry (http://www.allelefrequencies.net/) as being HLA‐B*15:02 carriers. We identified that both ASM cases (Fisher exact odds ratio [OR] = 21.5, p = 1.27e‐51) and ASM‐tolerant controls (OR = 5.65, p = 1.27e‐11) significantly differed from the population controls in HLA‐B*15:02 carrier frequency. This confirmed our biased selection of HLA‐B*15:02+ ASM‐tolerant controls. As far as the authors are aware, this is the largest known cohort of HLA‐B*15:02+ ASM‐tolerant controls.

ASM‐SJS/TEN cases versus all ASM‐tolerant subjects

After quality control, a total of 6 198 332 common variants were available for analyses. The genomic control value (λGC) for comparison between all ASM‐induced SJS/TEN and all ASM‐tolerant controls was 1.014, confirming minimal inflation (Figure 1A.1). Nine variants reached genome‐wide significance, with three being intergenic (chr3:73919920 [OR = 5.2, p = 3.6e‐08], chr8:21314026 [OR = 3.2, p = 3.3e‐08], chr12:131623246 [OR = 5.7, p = 2.0e‐09]) and the remainder falling within introns (chr4:820728 [OR = .3, p = 2.5e‐08], chr8:98263823 [OR = 4.0, p = 6.2e‐09], chr9:63832271 [OR = 6.5, p = 2.8e‐08], chr10:38643144 [OR = 5.4, p = 8.8e‐08], chr12:9426934 [OR = .3, p = 1.2e‐08]) of various genes (Figure 1A.2, Table 2). For two of these variants (chr12:9426934 and chr4:820728), the minor allele was underrepresented in ASM‐SJS/TEN. These variants were not in high LD with other significant or suggestive variants (Figure S4).

FIGURE 1

Quantile–quantile and Manhattan plots of the common variant association analysis testing 6 198 332 variants per group with at least one genome‐wide significant variant. (A) All ASMs (λ = 1.015). (B) HLA‐B*15:02‐negative cohort (λ = 1.017). (C) Carbamazepine cohort (λ = 1.014). (1) Quantile–quantile plots. Red line: expected = observed. (2) Manhattan plots. Red solid horizontal line = 5e‐08 (genome‐wide significance); Blue dotted horizontal line = 1e‐05 (suggestive variants). Lambda (λ) was calculated using the PLINK association function

TABLE 2

Variants with genome‐wide significance

Variant location	HGVS name	Variant reference cluster ID	Gene ID	A1	A2	Controls, homozygous A1/heterozygous/homozygous A2 (MAF)	Cases, homozygous A1/heterozygous/homozygous A2 (MAF)	OR (95% CI)	p
All ASMs, 113 cases vs. 84 controls
chr12:131623246	NC_000012.12:g.131623246T>C	rs4471527	Intergenic	T	C	0/14/70 (.08)	5/67/41 (.34)	5.7 (3.1–10.5)	2.04e‐09
chr8:98263823	NC_000008.11:g.98263831del	rs199755581	NIPAL2 (intronic)	CA	C	3/22/59 (.17)	36/26/48 (.45)	4 (2.5–6.5)	6.23e‐09
chr10:38643144	NC_000010.11:g.38643144A>G	rs1297852527	SLC9B1P3 (intronic)	G	A	0/14/70 (.08)	0/74/39 (.33)	5.4 (2.9–9.9)	8.75e‐09
	NC_000010.10:g.38936275A>G
	NG_025429.2:g.46033T>C
chr12:9426934	NC_000012.12:g.9426934G>C	rs77491650	DDX12P (intronic)	C	G	14/58/12 (.51)	2/49/62 (.23)	.3 (.2–.4)	1.21e‐08
chr4:820728	NC_000004.12:g.820783_820819ins	rs59567505 and rs143960439	CPLX1 (intronic)	G	GCACCCCTCACCAGCCTCACGTGAACCCCCAAGGTGGA	24/35/25 (.49)	14/22/77 (.22)	.3 (.2–.5)	1.49e‐08
chr9:63832271	NC_000009.12:g.63832271G>T	rs77542827	FRG1JP (intronic)	T	G	0/9/75 (.05)	0/61/52 (.27)	6.5 (3.1–13.6)	2.76e‐08
chr8:21314026	NC_000008.11:g.21314045_21314046insAGCTGGGAGTCAGTGAGAAAGAACAACACTGGGATCCAGTCCGG	rs778096762	Intergenic	C	CACACTGGGATCCAGTCCGGAGCTGGGAGTCAGTGAGAAAGAACA	14/23/47 (.3)	50/31/31 (.58)	3.2 (2.1–4.9)	3.29e‐08
chr3:73919920	NC_000003.12:g.73919922_73919923del	rs374138762	Intergenic	T	TTA	0/13/71 (.08)	0/69/44 (.31)	5.2 (2.8–9.9)	3.56e‐08
chr10:38643136	NC_000010.11:g.38643136C>T	rs879656274	SLC9B1P3 (intronic)	T	C	0/15/69 (.09)	0/73/40 (.32)	4.9 (2.7–8.9)	3.61e‐08
Carbamazepine, 85 cases vs. 77 controls
chr12:131623246	NC_000012.12:g.131623246T>C	rs4471527	Intergenic	T	C	0/13/64 (.08)	3/53/29 (.35)	5.8 (3.0–11.0)	1.36e‐08
HLA‐B*15:02‐negative, 35 cases vs. 53 controls
chr6:149444135	NC_000006.12:g.149444146_149444267del	rs1562468327	Intergenic	TCAGCCAGTGTGTCAGTCAGCCAGTGTTAGTCAGCCAGTGTGTCAGTCAGCCAGTGTCAGCCACCCAGTGTCAGTCAGCCAGTGTGTCAGCCACTGTCAGCCAATGTCAGCCAGTGTGTCAGC	T	2/7/44 (.10)	13/11/11 (.53)	9.7 (4.4–21.1)	5.90e‐10
chr8:98263823	NC_000008.11:g.98263831 CA>C a	rs199755581	NIPAL2 (intronic)	CA	C	0/11/42 (.10)	14/8/12 (.53)	9.7 (4.4–21.3)	6.84e‐10
chr21:9790175		Not reported in dbSNP	LINC01667 (intronic)	CCTCTCTCCAGGCTCACACATTGAAGAGAA	C	0/17/36 (.16)	16/9/10 (.59)	7.4 (3.7–15)	4.22e‐09
chr9:63832271	NC_000009.12:g.63832271G>T	rs77542827	FRG1JP (intronic)	T	G	0/3/50 (.03)	0/24/11 (.34)	17.9 (5.1–62.5)	1.45e‐08
chr17:26783109	NC_000017.11:g.26783109G>A	rs1286845082	Intergenic	A	G	0/2/51 (.02)	0/22/13 (.31)	23.8 (5.4–105.5)	2.28e‐08
chr17:26783113	NC_000017.11:g.26783113T>C	rs1597607761	Intergenic	C	T	0/1/52 (.01)	0/20/15 (.29)	42 (5.5–321.8)	3.13e‐08
chr3:183730061	NC_000003.12:g.183730061T>G	rs1391213386	YEATS2 (intronic)	G	T	1/3/49 (.05)	5/16/14 (.37)	11.9 (4.3–33.1)	3.26e‐08
chr4:39690137	NC_000004.12:g.39690141_39690142del	rs1211926109	Intergenic	AAT	A	16/18/19 (.47)	0/5/29 (.07)	.1 (.03–.2)	3.55e‐08

Gene ID indicates single nucleotide polymorphism identification.

Abbreviations: A1, PLINK assigned minor allele; A2, PLINK assigned major allele; ASM, antiseizure medication; CI, confidence interval; CPLX1, complexin 1; DDX12P, DEAD/H‐box helicase 12, pseudogene; FRG1JP, FSHD Region Gene 1 Family Member J, Pseudogene; HGVS, Human Genome Variation Society; LINC01667, Long Intergenic Non‐Protein Coding RNA 1667; MAF, minor allele frequency; NIPAL2, NIPA‐like domain containing 2; OR, odds ratio; SLC9B1P3, solute carrier family 9 member B1 pseudogene 3, YEATS2, YEATS Domain Containing 2.

The anchor position for this variant includes all nucleotides (repeats) potentially affected (HGVS is right‐shifted).

HLA‐B*15:02 stratification analysis

Among individuals who did not carry the HLA‐B*15:02 allele (35 cases and 53 controls), eight variants reached genome‐wide significance, including two, chr8:98263823 (OR = 9.7, p = 6.8e‐10) and chr9:63832271 (OR = 17.9, p = 1.5e‐08), that were identified in the primary analysis (Table 2). No variants reached genome‐wide significance among HLA‐B*15:02‐positive individuals.

Drug‐specific subanalyses irrespective of HLA‐B*15:02 status

Subanalysis of the carbamazepine‐induced SJS/TEN cases (n = 85) and tolerant controls (n = 77) identified a single genome‐wide significant variant, chr12:131623246 (OR = 5.8, p = 1.4e‐08), that had been identified in the primary analysis (Figure 1B, Table 2). No genome‐wide significant variants were identified in the subanalysis specific to phenytoin (11 phenytoin‐SJS/TEN cases and 40 phenytoin‐tolerant controls). This was likely due to the small sample sizes, as indicated in the quantile–quantile plot (Figure S5).

Significant variants differ from population controls

Comparison of significant variant allele frequencies in cases in our primary and subanalyses with those in East Asian population controls, as reported in the gnomAD genome aggregation database, found all to be significantly different (Table S2).

Interaction analysis of significant variants and HLA‐B*15:02 carrier status

The interaction analysis showed that carrying one of the significant variants identified in the primary analysis changed the likelihood of having developed SJS/TEN (Table 3, Table S3). HLA‐B*15:02 carriers who were homozygous for chr12.9426934 had .23 times the risk of those who carry only the wild‐type allele, whereas heterozygotes had .71 times the risk. For six variants, risk was increased irrespective of HLA‐B*15:02 status or zygosity. In those carrying HLA‐B*15:02, the relative risk ranged from 1.39 (heterozygous carriers of chr9:63832271) to 1.66 (heterozygous carriers of chr10:38643144). In individuals negative for HLA‐B*15:02, the relative risk ranged from 3.90 (heterozygous carriers of chr3:73919920) to 4.93 (heterozygous carriers of chr9:63832271). Presence of the chr8:98263823 variant, which reached genome‐wide significance in both the primary and HLA‐B*15:02‐negative subanalyses, increased risk for HLA‐B*15:02‐heterozygous carriers and for noncarriers of HLA‐B*15:02, irrespective of zygosity.

TABLE 3

ASM genome‐wide significant regulator variants alter probability and RR of ASM‐related Stevens–Johnson syndrome/toxic epidermal necrolysis

Variant	Neg PP, %	PP Neg + A1, %	RR a (p), Neg vs. Neg + A1	Pos PP, %	Pos PP + A1, %	RR a (p), Pos vs. Pos + A1
chr3:73919920	20.3	het = 79.3	het: 3.90 (<.001)	57.1	het = 86.8	het: 1.52 (<.001)
chr4:820728	68.4	het = 20.7 hom = 14.3	het: .30 (<.001) hom: .21 (<.001)	79.7	het = 57.1 hom = 64.7	het: .72 (.034) hom: .81 (.237)
chr8:21314026	18.4	het = 39.3 hom = 77.3	het: 2.13 (.062) hom: 4.19 (<.001)	60.0	het = 76.9 hom = 78.6	het: 1.28 (.136) hom: 1.31 (.064)
chr8:98263823	22.2	het = 42.1 hom = 100	het: 1.89 (.) hom: 4.50 (.)	67.9	het = 62.1 hom = 88.0	het: .91 (.) hom: 1.30 (.)
chr9:63832271	18.0	het = 88.9	het: 4.93 (<.001)	62.1	het = 86.0	het: 1.39 (.003)
chr10:38643136	17.0	het = 74.3	het: 4.37 (<.001)	55.4	het = 88.7	het: 1.60 (<.001)
chr10:38643144	15.4	het = 75.0	het: 4.87 (<.001)	54.4	het = 90.4	het: 1.66 (<.001)
chr12:9426934	75.7	het = 25.0 hom = 9.1	het: .33 (<.001) hom: .12 (<.001)	88.9	het = 62.7 hom = 20.0	het: .71 (.001) hom: .23 (<.001)
chr12:131623246	21.1	het = 71.4 hom = 100	het: 3.39 (<.001) hom: 4.75 (<.001)	53.7	het = 88.7 hom = 100	het: 1.65 (<.001) hom: 1.86 (<.001)

(.) indicates that value could not be estimated because all observations of chr8.98263823 in patients were Pos.

Abbreviations: A1, minor allele (see Table 1 for details); ASM, antiseizure medication; het, heterozygous; hom, homozygous; Neg, HLA‐B*15:02 negative; Pos, HLA‐B*15:02 positive; PP, predicted probability; RR, relative risk.

Calculated using the exponent of the coefficients; see Table S2 for details.

Predicted regulatory function of noncoding significant and suggestive variants

Among the 15 unique genome‐wide significant variants of the various analyses, only one had in silico predicted functional consequence. The chr4:820728 variant falls within an intronic region of the complexin 1 (CPLX1) gene. eQTL data associate the variant with expression of solute carrier family 26 member 1 (SLC26A1) in multiple tissues and cell types, including skin and cultured fibroblasts. sQTL data show potential to modify the splicing of two further genes in skin, alpha‐L‐iduronidase (IDUA) and transmembrane protein 175 (TMEM175; Table S4).

We then further explored the potential functional consequence of suggestive variants. In the primary analysis, the Manhattan plot identified a cluster of eight suggestive variants on chromosome 12 (Figure 1A.2). These variants fall within a 2.6‐kb region 94.6 kb downstream of the chr12:9426934 variant and 125 Mb downstream of chr12:131623246. They lie between the pleckstrin homology and RhoGEF domain containing G6 (PLEKHG6) gene and the CD9 molecule (CD9) gene (Figure 2A). This region is enriched for chromatin marks that are often found near regulatory elements, and six of the variants fall within promoter flanking regions at this locus. One also overlaps a region where the CCCTC‐binding factor (CTCF) binds to DNA. CTCF is a highly conserved transcription factor that can function as an activator, repressor, or insulator (blocking communication between enhancers and promoters). Tissue‐ and cell type‐specific eQTL data indicate that seven of the eight variants alter the expression levels of both genes in cultured fibroblast cells and testis tissue (Table S5).

FIGURE 2

Distinct clusters of suggestive variants. (A) Antiseizure medication cohort. (B) HLA‐B*15:02‐negative cohort. (C) HLA‐B*15:02‐positive cohort. (D) Carbamazepine‐exposed cohort. (E, F) Phenytoin‐exposed cohort. The reference single nucleotide polymorphism (SNP) is depicted as a pink diamond. The color coding for all other SNPs indicates linkage disequilibrium relative to the top SNP, as follows: red, r 2 ≥ .8; orange, .6 ≤ r 2 < .8; green, .4 ≤ r 2 < .6; blue, .2 ≤ r 2 < .4; light blue, r 2 < .2. Gray represents unknown n r 2 from the 1000 Genomes database. The bottom row represents the gene annotations.

In addition to the eight genome‐wide significant variants identified in the subanalysis of non‐HLA‐B*15:02 carriers, there were five suggestive variants that fall within an intronic region of the ATP binding cassette subfamily A member 1 (ABCA1) gene on chromosome 9 (Figure 2B). These five were predicted to modify expression of an RNA gene: ENSG00000226334.1 (Table S6).

In the analysis of HLA‐B*15:02 carriers, 17 suggestive variants were found to cluster on chromosome 2. These variants fall within an intronic region of dedicator of cytokinesis 10 (DOCK10; Figure 2C), and seven are predicted to alter the splicing ratio of DOCK10. In addition, 14 lie within a promoter flanking region of microRNA gene ENSR00000638436 (MIR4439), and four of these also overlap a CTCF binding site (Table S7).

In the subanalysis of individuals exposed to carbamazepine, seven suggestive variants overlap the HLA‐C and HLA‐B gene region on chromosome 6 (Figure 2D). Five of these variants are predicted to modify expression of HLA‐related molecules, including HLA‐C, MHC class I polypeptide‐related sequence B (MICB), and an HLA complex group 4 (HCG4) pseudogene (ENSG00000271581.1). All five variants were also predicted to alter the splicing of HLA‐E (Table S8).

In the analysis of phenytoin‐treated individuals, we identified suggestive variants on chromosome 6 (n = 36) and chromosome 11 (n = 16; Table S9). Variants on chromosome 6 are intergenic and nearby the CD24 gene. Nine of these variants overlap regulatory DNA regions and are predicted to alter CD24 expression in multiple tissues including the skin (Figure 2E, Table S9). The chromosome 11 variants lie within an intronic region of the RNA gene AP003066.1, and were predicted to modify expression of this RNA gene in regulatory T cells (Figure 2F, Table S9).

DISCUSSION

Our study interrogated >6 million common variants, considerably more than previous genome‐wide association studies of ASM‐induced SJS/TEN, which have typically included .5–1.2 million variants. , We also have one of the largest cohorts to interrogate the whole genome of ASM‐SJS/TEN cases and matched drug‐tolerant controls. The previous papers by McCormack et al. included 22 hypersensitivity syndrome subjects (12 with SJS/TEN) and 2691 healthy controls for their association. Unlike previous studies, we purposely recruited a large proportion of HLA‐B*15:02 carriers who tolerated treatment. This large number of tolerant individuals with the HLA‐B*15:02 allele enabled secondary analyses stratified by HLA‐B*15:02 carrier status. We were able to perform treatment‐specific analysis for carbamazepine and phenytoin only, as there were too few cases for other drugs, and could not stratify analysis based on HLA‐B*15:02 carrier status by drug. The analyses identified novel significant and suggestive variants located in non‐protein‐coding regions of the DNA, some of which have evidence to support functional effect on cell type‐specific gene expression.

Despite our relatively small sample size compared to other genome‐wide association studies, the primary analysis identified nine significant variants. This seemingly large number of variants is not unusual for detecting pharmacogenomic responses that are not as tightly controlled by natural selection, unlike disease‐related variants, and aligns with our power calculation. This has been observed in the previous genome‐wide association studies in this area with a large number of variants despite having a small sample size. , Based on the lambda value of 1.015, we concluded that the variants identified were not being driven by population stratification, as it was less than the commonly applied 1.05 threshold. However, there may be cryptic population structure that is unresolvable with the current methodology. Interestingly, in HLA‐B*15:02 carriers, two of these variants reduced the relative risk of developing SJS/TEN and might explain why some of these carriers do not go on to develop SJS/TEN following treatment with a high‐risk drug. These findings also suggest that carrying specific combinations of significant variants can modify predisposition to a severe drug reaction.

Our study specifically analyzed risk in a cohort of HLA‐B*15:02‐negative individuals. Although rarer, as the negative predictive value approaches 100% without HLA‐B*15:02, some of these individuals are nevertheless at risk of SJS/TEN. This analysis identified eight variants that can potentially explain this risk, including one, chr8:98263823, that also reached genome‐wide significance in the primary analysis. Interaction modeling showed that risk for HLA‐B*15:02‐negative individuals was increased, irrespective of zygosity, if they carried this variant, or one of six other genome‐wide significant variants identified in the primary analyses.

Evidence for functional analysis was found for only one of the 15 unique significant variants identified: the chr4:820728 variant that falls within an intronic region of CPLX1. This variant contributes to expression variation of multiple nearby genes, including SLC26A1 and IDUA, in several tissues including skin. It is feasible that a variant impacts on the function of distal genes rather than the gene it falls within; for instance, the modulation of IRX family gene expression by intronic variants within the "obesity" FTO gene has been well described.

Although these variants have relatively high OR for a whole genome sequencing study, they may only marginally alter the absolute risk of SJS/TEN due its overall rare occurrence in the population. Translating the OR of any association into predictive value needs to be interpreted in the context of the pretest probability of the condition in the test population. Although the negative predictive value of HLA‐B*15:02 is very high for carbamazepine‐induced SJS/TEN, this has not been demonstrated for SJS/TEN induced by other aromatic ASMs, and approximately one third of our cases did not carry the HLA‐B*15:02 allele. The positive predictive value of HLA‐B*15:02 for carbamazepine‐induced SJS/TEN is estimated to be <10%. Therefore, additional risk variants identified in this population may still be clinically meaningful in refining risk stratification, and this dataset highlights why segregating HLA‐B*15:02 carrier status is important to further understanding the complexities of SJS/TEN inheritance.

The quartile–quartile plot indicated enrichment of variants in subthreshold regions, which suggests a polygenic inheritance pattern. This pattern has been observed in other complex inherited traits such as schizophrenia and human height. In our study, the quartile–quartile plot variant enrichment indicated that a cluster of variants impacts the same genomic region, but need not be the same in all individuals. The polygenic architecture does not yield specific candidates for clinical screening; however, it provides insights into candidate DNA regions with biological importance for risk of ASM‐induced SJS/TEN.

The primary analysis identified an area of enrichment for suggestive variants on chromosome 12 between the CD9 and PLEKHG6 genes. eQTL and sQTL data show that the variants potentially impact the regulation of CD9 and PLEKHG6 expression in skin tissue. Notably, CD9 is an exosomal marker, and CD9‐expressing exosomes identified in the serum of SJS/TEN patients have recently been published. The authors attribute keratinocyte cell death to ferroptosis, which is a type of iron‐dependent regulated cell death, whereby the deleterious factors are delivered in exosomes internalized by keratinocytes. Our data suggest that CD9 may be an important gene for SJS/TEN caused by aromatic ASMs.

The suggestive variants identified in the HLA‐B*15:02‐positive subanalysis might affect splicing of the DOCK10 gene. The encoded protein is a member of the DOCK family of Rho GTPase activators and may have a role in activating T cells and the development of CD8+ T cells. DOCK10 is also differentially downregulated in contact dermatitis when presented with sensitizers. Our data suggests that genetic variation may modify risk of developing ASM‐induced SJS/TEN via perturbation of DOCK10 activity in HLA‐B*15:02‐positive individuals.

Suggestive variants identified in the carbamazepine subanalysis were predicted to impact the function of HLA molecules. In particular, five variants have the potential to modify the splicing ratio of HLA‐E, a nonclassical MHC class I molecule. Peptides can be presented on HLA‐E either via the CD94:NKG2 receptor or the T‐cell receptor. Keratinocytes from affected skin in patients with SJS/TEN and other exanthemas have been shown to express HLA‐E and be sensitized to killing by cytotoxic T lymphocytes expressing the CD94/NKG2C receptor. Although HLA‐B*15:02 remains a robust marker of carbamazepine‐induced SJS/TEN reactions, these findings suggest that genetic variants within regulatory DNA regions of chromosome 6 may indirectly contribute to risk of pathogenesis by modifying the function of additional HLA molecules.

The subanalysis of individuals exposed to phenytoin identified suggestive variants with potential to alter the expression of CD24 in skin tissue. Interestingly, CD24 has been found to be significantly differentially expressed with dithranol treatments in keratinocytes from psoriasis patients. CD24 can be expressed by many immune cells, including T cells, where it acts as a costimulatory molecule. B cells expressing CD24+ from patients with blistering autoimmune skin disease (bullous pemphigoid) fail to repress CD4+ T‐cell inflammatory proliferation. If mice CD4+ T cells expressed CD24, they were not protected from concanavalin A‐induced liver injury. Therefore, genetic variation that modifies the function of CD24 in either T cells, B cells, or keratinocytes may contribute to susceptibility to phenytoin‐induced SJS/TEN.

Information on the potential functional consequence of the non‐protein‐coding variants identified is limited. Here, we have made inferences about possible consequence based on the most reliable data available. Functional studies that could include CRISPR/Cas9, ChIP‐seq experiments, or functional genomic approaches (i.e., differential expression studies utilizing either transcriptomics or proteomics) are needed to validate our conjectures.

The functional consequences associated with suggestive variants support the conjecture that non‐protein‐coding variants modify susceptibility to ASM‐induced SJS/TEN. It should also be noted that current databases report on potential eQTL effects at baseline conditions only, and the functional consequences of variants that become apparent under certain dynamic conditions such as drug exposure, which have been termed “response eQTL,” remain elusive. The genes implicated in this study warrant further functional interrogation to better understand the molecular mechanisms underlying ASM‐induced SJS/TEN.

In conclusion, our analysis highlights that mixing HLA‐B*15:02 carrier status in previous studies might have masked variants and genes contributing to susceptibility. Our interaction analyses suggest that presence of additional variants may reduce risk conferred by HLA‐B*15:02, explaining why some carriers tolerate ASMs. These findings, coupled with variant enrichment observed in our quantile–quantile plots, suggest that inheritance of risk in SJS/TEN may involve multiple markers of small effect size, as has been observed in other complex inherited traits such as human height. In this scenario, the combination of risk variants may be heterogeneous across those susceptible, but may disrupt pathways common to the pathology.

CONFLICT OF INTEREST

Z.C. is supported by an Early Career Fellowship from the National Health and Medical Research Council of Australia (GNT1156444), and he/his institution has received consultancy fees and/or research grants from Arvelle Therapeutics and UCB Pharma. None of the other authors has any conflict of interest to disclose.

AUTHORS CONTRIBUTIONS

K.A.M. processed, analyzed, and interpreted the whole genome sequence data and was a major contributor in writing the manuscript. A.A. also contributed to whole genome sequence data processing, analysis, and interpretation, as well as contributing to the writing process. Z.C. contributed to the analysis and interpretation of the Stata results. P.K., Y.‐W.S., C.‐C.N., L.B., S.C., S.P., P.C.S., K.‐S.L., and W.‐P.L. aided with the acquisition of samples, project design, and draft revisions. P.K. conceived and supervised the study. J.‐H.D. contributed to the whole genome sequence data processing and initial analysis. All authors read and approved the final manuscript.

Fig S1

Click here for additional data file.

Fig S2

Click here for additional data file.

Fig S3

Click here for additional data file.

Fig S4

Click here for additional data file.

Fig S5

Click here for additional data file.

Supplementary Material

Click here for additional data file.

Supplementary Material

Click here for additional data file.