Assessment of Molecular Subtypes in Thyrotoxic Periodic Paralysis and Graves Disease Among Chinese Han Adults

Key Points Question Is thyrotoxic periodic paralysis (TPP) a molecular subtype of Graves disease? Findings In this case-control study in a Chinese Han population, 5 TPP susceptibility loci were identified, including 3 specific loci and 2 loci shared by Graves disease and TPP. The ratio of persistent thyrotropin receptor antibody positivity was higher in TPP than in Graves disease, and TPP could be predicted from Graves disease using TPP-specific loci. Meaning A complete genetic architecture will be helpful to understand the pathophysiology of TPP, and a useful prediction model could prevent the onset of TPP, suggesting TPP as a molecular subtype of Graves disease.


The characters of study participants
A total of 537 patients with TPP, 1,519 patients with GD having no history of TPP, and 3,249 healthy participants were recruited from Chinese Han population through collaboration with the hospitals in China (Table 1). All samples were obtained under local institutional review board approval and with documented informed consent.
All the cases were sporadic. We collected 5-ml blood samples from each participant for DNA preparation and biochemical measurements.
TPP was diagnosed in patients with episodic flaccid paralysis, serum potassium concentration less than 3.5 mmol/liter (normal concentration, 3.5-5.0 mmol/liter) during at least one of the attacks, and thyrotoxicosis.
Diagnosis with GD was based on documented clinical and biochemical evidence of hyperthyroidism, diffuse goiter, and the presence of at least one of the following clinical characters: positive thyroid-stimulating hormone receptor antibody (TRAb) tests, exophthalmos, or pretibial myxedema. 1-5 All individuals classified as having TPP or GD were interviewed and examined by the experienced clinicians.
Control participants were over 35 years old. Given that TPP, GD, and other autoimmune thyroid disease have preponderance in the young female population, this age criteria could reduce the number of controls who might develop GD later. For excluding clinical or subclinical autoimmune thyroid disease, the levels of sensitive TSH and TPOAb in control participants were measured using chemiluminescence immunoassay in our laboratory. 1,2

Evaluation of population structure and Quantile-quantile plots
To evaluate the population structure in samples, principal component analysis (PCA) and multidimensional scaling (MDS) analysis using a subset of pruned and unlinked markers were performed by SmartPCA (one of the program modules in the EIGENSOFT software package) 7 and PLINK 6 , respectively. For the PCA, our samples and the HapMap samples were plotted using the first two eigenvectors produced by smartPCA (eFigure 2A in the Supplement). For MDS analysis, dimensions were calculated based on the IBD pairwise distance among all individuals in both our cohort and the HapMap participants using PLINK. 6 The first two dimensions of the result were plotted (eFigure 2B in the Supplement). SmartPCA was also used to identify potential genetic outliers. The top ten principal components were computed using HapMap populations followed by projection of current cohort participants onto those principal components, and then run the outlier removal (σ threshold =6 with five iterations).
Finally, all participants for the association analysis were restricted to Chinese Han Populations ancestry.
The distribution of observed P values (on the -log10 scale) of given SNPs were plotted against the theoretical distribution of expected P values to construct quantile-quantile plots. To account for relatedness and stratification within our TPP cases and healthy control sample sets, we calculated the genomic control inflation factor (λ) by dividing median χ2 statistics by 0.4563 based on chip markers. The calculations and the plots were done using PLINK and R statistics packages (eFigure 3 in the Supplement).

Imputation and statistical analysis in the discovery stage
The genotype imputation for the GD GWAS section has been carried out as described previously. 1,2 And genotype imputation for the TPP GWAS section was also performed using IMPUTE2 software 8 in the 158 TPP cases and 803 controls. The 1000G phase 1 integrated variant set (Mar 2012) were also used as a reference. For the imputed SNPs, we first removed the indels and SNPs with a relatively low confidence (estimated probability ≤ 0.9). Then the total of 33,676,913 SNPs with genotype call rates below 98%, or with MAF less than 0.01, or with significant deviation from Hardy-Weinberg equilibrium in the controls (P < 10 -6 ) or in cases (P < 10 -10 ), were excluded. Finally, the association analysis for 2,752,055 SNPs among the 171 GD cases with TPP, 1,404 GD cases without TPP, and 2,160 healthy controls were carried out using the SNPTEST v2 software (Frequentist Association Tests with score method). 9 For the 2,692,307 autosomal SNPs in 171 GD cases with TPP and 2,160 controls, or 1,404 GD cases without TPP, and 2,160 healthy controls, we also carried out the association analysis using the Cochran-Armitage trend test in PLINK, 6 respectively. Whereas, we analyzed 59,748 SNPs located in X chromosome using the logistic regression in PLINK according to the following principles: for alleles A and a, males were coded A to 0 and a to 2, and females were coded AA to 0, Aa to 1, and aa to 2, and additionally sex (male =0, female =1) was also automatically included as a covariate. 6 Genome-wide P-value and regional plots were generated using R software (LocusZoom) or SNAP, version 2.2 software. 10,11 SNPs selection for the second stage study Among 2,752,055 SNPs in the discovery stage, there are 1,731 autosomal SNPs with discovery stage P < 0.0001 and the P value between patients with TPP and patients with GD having no history of TPP less than 0.05. Firstly, we defined the chromosomal regions based on HapMap recombination rates (the recombination rate more than 30 cM/Mb) and found 1,731 SNPs with discovery stage P < 0.0001 located in 163 chromosomal regions. Next, we performed the forward logistic regression analysis using SNPs with the discovery stage P < 0.01 to identify the potential causal SNPs using R package. Then, we performed the sample size analysis and further selected 104 SNPs which can meet the genome-wide association significance level (P ≤ 5 × 10 -8 ) within the samples size of TPP cases less than 1,000. Also, we exclude 3 SNPs without significant difference for the MAF between 171 GD cases with TPP and 1,404 GD cases without TPP. Given the association of 17q24.3 locus with TPP has been reported in the previous paper, a total of 12 tagSNPs at 17q24.3 (including rs312691) were also selected and genotyped in the second stage study (eTable 2 in the Supplement). Unfortunately, the probe of rs79449227 at donated by 25 individuals from the same group for gene expression assay in distinct subpopulations of PBMCs.
The CD4+, CD8+, CD14+ and CD19+ subsets of PBMCs were isolated using MACS Column kits (Miltenyi Biotec) according to the manufacturer's instructions. The purity of each cell subpopulation was determined by an LSR II Flow Cytometer (BD Biosciences), and the cell subpopulations with a purity of over 90% were used for real time  Real time quantitative reverse transcription PCRs (RT-PCRs) for a series of genes or transcript units were performed in duplicate using the SYBR Green and ABI 7900HT Fast Real-Time PCR System.
Expression levels in all samples were normalized to the relative expression level of GAPDH. Primers used for the amplifications were shown in Table S8. We performed statistical analysis of expression data using ANOVA and an unpaired Student's t-test (the two-tail P-value was indicated on the figures). The genotypes of the SNPs were determined using the ABI 7900HT Fast Real-Time PCR System.

Sequencing of two TPP susceptibility genes and gene level association analysis
To further evaluated the association of two TPP susceptibility genes (DCHS2 and KCNJ2) and a novel long intergenic noncoding RNA (lincRNA, CTD-2378E21.1) with TPP, we first amplified and sequenced the exonic regions of three TPP genes in 34 patients with TPP and 102 healthy controls using PCR and Sanger sequencing and the primers were provided in the online eTable 10 in the Supplement. Next, we compared the frequency difference of rare nonsynonymous variants (MAF ≤ 1%), based on 1000 Genomes East Asian (EAS), ESP6500, ExAc EAS, dbSNP 139 and our inhouse database et al, in the three genes using Fisher's method to combine burden test.

The prediction models
We constructed the prediction model using the weighted genetic risk score (wGRS) 16 and two different markers.
Firstly, three independent TPP specific SNPs (rs1352714, rs2186564, and rs312691) in three chromosomal loci were used as the markers to construct the wGRS model. Then, a total of 11 independent SNPs in eight chromosomal loci with Preplication < 0.05, with Pcombined < 0.0005 (Bonferroni-corrected significance in combined populations), and with PTPPvsGD < 0.05, were used as the markers to construct the wGRS model. As for the wGRS model, a weighted score was calculated by multiplying each participant's allele score (0, 1, 2) by the SNP's relative effect size (the natural log of OR in our study) obtained from our association study. Then, the value of the weighted score was rescaled by dividing all values by the sum of the effect sizes and then multiplying by the total number of SNPs, thus obtaining the final weighted GRS. 16 Next, logistic regression was used to calculate the OR and P values for each wGRS and gender as a covariant. Subsequently, we generated the receiver operator characteristic (ROC) curves and calculated the area under the curve (AUC), the sensitivity and specificity of each model to determine how well the models discriminates between TPP cases and Patients with GD without TPP.
Finally, we categorized the risk scores as four different groups according to the mean wGRS and S.D.: group one, wGRS less than the mean; group two, wGRS between the mean and the mean + 1 S.D.; group three, wGRS from the mean + 1 S.D. to the mean + 2 S.D.; group four, wGRS greater than the mean + 2 S.D. P values, OR and 95% confidence intervals were evaluated using group 1 as the reference (eTable 9 in the Supplement).

Web resources
Software and Web resources utilized in the study are provided below: