Context

Human height is an inheritable, polygenic trait under complex and multi-locus genetic regulation. Familial short stature (FSS; also called genetic short stature) is the most common type of short stature and is insufficiently known.

Objective

To investigate the FSS genetic profile and develop a polygenic risk predisposition score for FSS risk prediction.

Design and Setting

The FSS case group of Han Chinese ancestry was diagnosed by pediatric endocrinologists in Taiwan.

Patients and Interventions

The genetic profile of 1,163 FSS cases was identified by using a bootstrapping sub-sampling and genome-wide association studies (GWAS) method.

Main Outcome Measures

Genetic profile, polygenic risk predisposition score for risk prediction.

Results

Ten novel genetic SNPs and 9 reported GWAS human height-related SNPs were identified for FSS risk. These 10 novel SNPs served as a polygenic risk predisposition score for FSS risk prediction (area under curve (AUC): 0.940 in the testing group). This FSS polygenic risk predisposition score was also associated with the height reduction regression tendency in the general population.

Conclusion

A polygenic risk predisposition score composed of 10 genetic SNPs is useful for FSS risk prediction and the height reduction tendency. Thus, it might contribute to FSS risk in the Han Chinese population from Taiwan.

Ethics and consent

This study is a cross-sectional study on the clinical, biochemical, and genetic findings collected from cases of familial short stature (FSS) sequentially identified from the Children’s Hospital, China Medical University, Taichung, Taiwan, from August 1999 to September 2018. This study was approved by the institutional review board and the ethics committee of Human Studies Committee of China Medical University Hospital. Written informed consent was obtained from the participants, their parents, or legal guardians according to institutional requirements and Declaration of Helsinki principles (https://www.wma.net/policies-post/wma-declaration-of-helsinki-ethical-principles-for-medical-research-involving-human-subjects/).

Participants

The case group used in this study was diagnosed by pediatric endocrinologists in Taiwan (Fig. 1 and Supplementary Fig. S1). The case group included FSS cohort (N = 1,163). All recruited cases were of Han Chinese ancestry and diagnosed with FSS (1,20,22). The selection criteria for FSS was (1) height less than the 3rd percentile (Supplementary Fig. 7), (2) Their fathers’ and/or mothers’ height less than the 3^rd percentile, (3) bone age appropriate for chronologic age, (4) normal onset of puberty, (5) normal annual growth rate, and (6) normal results of clinical biochemistry examination. Excluded individuals were those with all other abnormal morphology and karyotyping results, abnormal bone age or puberty stage, or who had abnormal serum or plasma levels of clinical biochemistry examinations for GH-IGF1 axis, thyroid function or precocious puberty (Supplementary Fig. S1). Finally, 1,163 FSS subjects were served as the cases in the study.

In this study, the controls consisted of 4,168 individuals from the Taiwan Biobank (TW-Biobank; http://www.twbiobank.org.tw/new_web/index.php) and our type 2 diabetes cohorts (Fig. 1 and Supplementary Fig. 6). Furthermore, the selection criteria for controls with their top 75^th percentile of human height in our study (N = 1,071) was (1) no history of FSS diagnosis, (2) height exceeding that of the top 75^th percentile of the general population in Taiwan, and (3) age <61 years. All case and control groups in this study were of Han Chinese origin based on principal component analysis of genome-wide data (Supplementary Fig. S2).

Genotyping and quality control

Genomic DNA was extracted from the blood samples of participants according to standard protocols using the Qiagen genomic DNA isolation kit (Qiagen, Taichung, Taiwan). Each FSS case (N = 1,163) was genotyped at the National Genotyping Centre at Academia Sinica (Taipei, Taiwan) using the Axiom genome-wide CHB array plate, according to the manufacturer’s procedure. For the control group from Taiwan Biobank, the GWAS data of each sample was genotyped using the Axiom genome-wide TWB array plate. For the control group from our type 2 diabetes cohorts, the GWAS data of each type 2 diabetes patient was genotyped using the Axiom genome-wide TWB array plate, the Affymetrix genome-wide human SNP array 6.0, and the Illumina HumanHap550-Duo BeadChip according to the manufacturer’s procedure.

Because GWAS data were from different genotyping platforms, genotype imputations were performed in both FSS cases and controls according to a two-step genotype imputation approach. We used SHAPEIT2 to pre phase the study genotypes into full haplotypes(25). We then performed imputation using IMPUTE2 and Phase I 1000 Genomes Project reference panel (June 2011 interim release) consisting of 1,094 phased individuals from multiple ancestry groups [The 1000 Genomes Project Consortium, 2010](26). Finally, we used the GTOOL software (http://www.well.ox.ac.uk/~cfreeman/software/gwas/ gtool.html) to homogenize strand annotation by merging the imputed results obtained from each set of genotyped data.

Genotype and imputed genotype data were quality controlled and genetic variants were excluded for further analysis if (1) only one allele appeared in cases and/or controls; (2) the total call rate was less than 95% for both cases and controls; (3) the minor allele frequency was less than 0.5% in the controls in the Han Chinese population; (4) genetic variants significantly departed from Hardy-Weinberg equilibrium proportions (p < 0.01).

Genetic predisposition score calculation

The genetic predisposition score (also known as polygenic risk score or genetic risk score) is calculated by multiplying each beta-coefficient (log OR) value by the number of the corresponding risk allele under the additive model for each individual and then summing the products for the risk alleles identified from the multiple susceptible genetic variants (27). In this study, the genetic predisposition score was calculated based on the 10 genetic variants (SNPs). Each genetic variant was given a weightage based on the average effect size (beta-coefficient) for the FSS obtained from our study (Table 1). The genetic predisposition score was calculated by multiplying each beta-coefficient by the number of corresponding risk alleles (risk allele homozygote (the risk genotype is coded as “2”), risk allele heterozygote (the risk genotype is coded as “1”), and non-risk allele homozygote (the non-risk genotype is coded as “0”) according to the additive inheritance model) and then summing the products from these 10 genetic variants weighted by their estimated effect sizes (log OR).

Statistical analysis

All the genotyped and imputed GWAS results of FSS cases and their controls were used for association studies using a regression framework implemented in PLINK under the additive inherited genetic model (28). The difference in allelic frequency in the additive model between the cases and controls were measured by odds ratios (ORs) with 95% confidence intervals (CIs) using logistical regression models (Tables 1-2 and Supplementary Tables 1-6 and 8-9). All data management and statistical analyses were performed using PLINK and SAS software (version 9.4; SAS Institute, Cary, NC, USA).

For haplotype block analysis, the Lewontin D′ and R² values were used to evaluate the intermarker coefficient of linkage disequilibrium (LD) in both FSS cases and controls (29). The confidence interval for LD was estimated using a resampling procedure and was used to construct the haplotype blocks (Supplementary Fig. S3) (30,31).

The risk prediction model predicts the health outcome by using several predictor variables based on the observed patient’s characteristics (32). Risk prediction was evaluated by the area under the receiver operating characteristic (ROC) curves (AUCs). The AUC ranged from 0.5 (total lack of discrimination) to 1.0 (perfect discrimination). AUCs were calculated for the predicted risks of 10 novel, 9 reported, and combined SNPs, respectively (Figs. 2a and 2b).

For the linear human height curve model, participants with the genetic predisposition score calculated from the 10 novel SNPs (Fig. 3a) and 9 human height-related SNPs (Fig. 3b) were used as continuous variables in a linear regression with human height (cm) as the dependent variable, respectively.