Height-Related Polygenic Variants Are Associated with Metabolic Syndrome Risk and Interact with Energy Intake and a Rice-Main Diet to Influence Height in KoGES

Adult height is inversely related to metabolic syndrome (MetS) risk, but its genetic impacts have not been revealed. The present study aimed to examine the hypothesis that adult height-related genetic variants interact with lifestyle to influence adult height and are associated with MetS risk in adults aged >40 in Korea during 2010–2014. Participants were divided into short stature (SS; control) and tall stature (TS; case) by the 85th percentile of adult height. The genetic variants linked to adult height were screened from a genome-wide association study in a city hospital-based cohort (n = 58,701) and confirmed in Ansan/Ansung plus rural cohorts (n = 13,783) among the Korean Genome and Epidemiology Study. Genetic variants that interacted with each other were identified using the generalized multifactor dimensionality reduction (GMDR) analysis. The interaction between the polygenic risk score (PRS) of the selected genetic variants and lifestyles was examined. Adult height was inversely associated with MetS, cardiovascular diseases, and liver function. The PRS, including zinc finger and BTB domain containing 38 (ZBTB38)_rs6762722, polyadenylate-binding protein-interacting protein-2B (PAIP2B)_rs13034890, carboxypeptidase Z (CPZ)_rs3756173, and latent-transforming growth factor beta-binding protein-1 (LTBP1)_rs4630744, was positively associated with height by 1.29 times and inversely with MetS by 0.894 times after adjusting for covariates. In expression quantitative trait loci, the gene expression of growth/differentiation factor-5 (GDF5)_rs224331, non-SMC condensin I complex subunit G (NCAPG)_rs2074974, ligand-dependent nuclear receptor corepressor like (LCORL)_rs7700107, and insulin-like growth factor-1 receptor (IGF1R)_rs2871865 was inversely linked to their risk allele in the tibial nerve and brain. The gene expression of PAIP2B_rs13034890 and a disintegrin and metalloproteinase with thrombospondin motifs-like-3 (ADAMTSL3)_rs13034890 was positively related to it. The PRS was inversely associated with MetS, hyperglycemia, HbA1c, and white blood cell counts. The wild type of GDF5_rs224331 (Ala276) lowered binding energy with rugosin A, D, and E (one of the hydrolyzable tannins) but not the mutated one (276Ser) in the in-silico analysis. The PRS interacted with energy intake and rice-main diet; PRS impact was higher in the high energy intake and the low rice-main diet. In conclusion, the PRS for adult height interacted with energy intake and diet patterns to modulate height and was linked to height and MetS by modulating their expression in the tibial nerve and brain.


Introduction
Adult height is achieved by growth during childhood and adolescence, with full height attained at age 18-20. The adult height of a child or adolescent between the age of 4 to 17.5 years can be predicted based on the current height, body mass, chronological age, and parents' height. The genetic impact on height increases can be estimated from adult height. However, after 30 years of age, people, especially women, gradually lose height. A loss of height of about 2.54 cm for men and 5.08 cm for women occurs between ages 30 to

Adult Height Criteria
The short stature (SS; n = 51,165; control) and tall stature (TS; n = 7536; case) groups of the participants were divided by the cutoff of ≥175 cm for men and ≥163 cm for women as tall stature, and below that for short stature. These heights were the 85th percentiles of adult height for people aged >40 between 2010-2014 in Korea.

Survey Questionnaires and Anthropometric and Biochemical Measurements
On their first visit, the participants completed the demographic questionnaires. Height and weight were measured with an automatic digital machine wearing light gowns and bare feet. The appendicular skeletal muscle mass in the city hospital cohort was predicted using a prediction model generated by the XGBoost algorithm in the Ansan/Ansung cohort [14]. Blood pressure was determined with a sphygmomanometer in a sitting position three times, and average systolic (SBP) and diastolic blood pressure (DBP) were recorded. The participants fasted for over 12 h. The serum and plasma were collected from the blood drawn without and with ethylenediaminetetraacetic acid (EDTA) and heparin, respectively. Biochemical parameters in the serum were measured with previously described standard methods [15]. They included glucose, triglyceride, total cholesterol, HDL, alanine aminotransferase (ALT), aspartate aminotransferase (AST), and creatinine concentrations in the serum, and blood hemoglobin A1c (HbA1c) and white blood cell (WBC) contents. The estimated glomerular filtration rate (eGFR) was calculated using the equation generated by the Modification of Diet in Renal Disease study [16]. The insulin resistance was calculated using the prediction model generated by the homeostatic model assessment of insulin resistance (HOMA-IR) in the Ansan/Ansung cohort from a previous study [17].
Details of lifestyle factors, including food intake, were collected using questionnaires during a health interview by a trained technician. Daily alcohol and coffee intakes were calculated by multiplying their intake amounts with frequencies. Smoking status was classified into non-, past, and current smokers. Smokers were specified as smoking >100 cigarettes in their lifetime, and current smokers who had not smoked for 6 months were considered past smokers. Exercise status was classified into regular exercise and no exercise based on the cutoff of 30 min of moderate-intensity activity such as water pushing a lawn mower, brisk walking, hiking, aerobics, riding a bike, dancing, doubles tennis, or rollerblading for ≥3 days per week.

Usual Food Intake Measurement
The usual food and nutrient consumption during the last 12 months was determined based on the response to a semi-quantitative food frequency questionnaire (SQFFQ) generated for Koreans. It included 106 foods Koreans commonly consumed. The SQFFQ was confirmed with three-day food records for the four seasons in a previous study [17]. The food frequencies were divided into 8 categories, such as seldom, once per month, two to three times per month, once or twice per week, three or four times per week, five or six times per week, daily, twice a day, and ≥3 times a day. One portion size of each food was classified into either more than, equal to, or half of the portion size shown in a photograph. Participants selected the frequencies and one portion size of each food item. Their intake Nutrients 2023, 15, 1764 4 of 23 was calculated by multiplying the median of the frequencies by one portion size of each food. Daily nutrient intake was estimated from the calculated daily food amount using the computer-aided nutritional analysis program CAN-Pro 2.0 (Korean Nutrition Society, Seoul, Republic of Korea). Daily nutrient intake was calculated by summing the individual nutrient of each food.

Dietary Patterns and Dietary Inflammatory Index (DII)
The dietary patterns were clustered with principal component analysis (PCA) of 30 predefined groups categorized from the 106 food items. Four dietary patterns were designated based on eigenvalues >1.5, orthogonal rotation procedure (varimax), and ≥0.40 factorloading values. The four dietary patterns were named as Korean-balanced diet (KBD), Western-style diet (WSD), plant-based diet (PBD), and rice-main diet (RMD), according to foods. The dietary inflammation index (DII) was calculated with a prediction equation generated with the food and nutrient intake and divided by 100 [18].

Genotyping, Its Quality Control, and Genotype-Tissue Expression (GTEx)
Genomic DNA was isolated from the volunteers' blood and genotyping was performed using a Korean Chip (Affymetrix, Santa Clara, CA, USA). It was made to evaluate the genetic impact of single nucleotide polymorphisms (SNPs) on metabolic diseases in Koreans at the Center for Genome Science at the Korea National Institute of Health [19]. The genotypes were imputed based on the Korean Haplotype Map (HapMap) data [20]. The genotyping accuracy was assessed with a Bayesian learning algorithm for Robust General Linear Models (RGLMs) [21]. Genotyping accuracy was checked with the exclusion criteria as follows: <98% genotyping accuracy, ≥30% heterozygosity, ≥4% missing genotype call rate, ≤1% minor allele frequency (MAF), and p ≤ 0.05 Hardy-Weinberg equilibrium (HWE) [21]. The significance distribution of the genetic variants was visualized with a Manhattan plot by the R program. A Q-Q probability plot indicates the goodness of fit of the actual data distribution to the theoretical data distribution. When the lambda value in the Q-Q plot was close to 1, the genotypes were ideal. Among the selected genotypes, the pathway associated with height was selected with satisfying the P value for Bonferroni correction using the multi-marker analysis of genomic annotation (MAGMA) gene-set analysis in the SNP2GENE function of the functional mapping and annotation (FUMA) web application (https://github.com/Kyoko-wtnb/FUMA-webapp/, accessed on 20 April 2022).
Genotype-Tissue Expression (GTEx) was used to show the gene expression with different genetic variant alleles in various tissues using the GTEx expression quantitative trait loci (eQTL) calculator (https://gtexportal.org/home/testyourown, accessed on 3 May 2022). Figure 1 provides the optimal genetic model selection process. First, genetic variants (n = 1299) influencing height were selected using the SS and TS groups in the city hospital-based cohort at a significance level of p < 5 × 10 −7 , as defined above. Of these, 956 genetic variants satisfied MAF (>1%) and HWE (p > 0.05), and the 79 unduplicated gene names were identified using g:Profiler (https://biit.cs.ut.ee/gprofiler/snpense, accessed on 7 March 2022). The genetic variants with high linkage disequilibrium (LD; D ≥ 0.2) were excluded due to giving the same genetic impact on the height. After exclusions, 37 genetic variants remained to meet the criteria (D < 0.2) using Haploview 4.2 in the PLINK toolset, and 15 with the same or unidentified gene names were eliminated. Nutrients 2023, 15, x FOR PEER REVIEW variants remained to meet the criteria (D' < 0.2) using Haploview 4.2 in the PLINK and 15 with the same or unidentified gene names were eliminated. Figure 1. Flow chart to generate polygenic risk score (PRS) associated with stature by SN interaction and its interaction with lifestyle factors. The short stature (SS; n = 51,165; cont tall stature (TS; n = 7536; case) groups of the participants were divided by the cutoff of ≥17 men and ≥163 cm for women, which were the 15th percentiles of adult height in people a years in 2010-2014 in Korea.

Selection of Genetic Variants to Influence Adult Height and Their Optimal Model with the SNP-SNP Interaction
Ten SNPs were included in the optimal genetic model by selecting them ba their interaction with each other from 22 genetic variants in the generalized mul dimensionality reduction (GMDR) software version 2.0 using the exhaustive sear after adjusting for covariates. Covariate set 1 contained age, gender, residence are cation, and income, and covariate set 2 included covariates in set 1 plus energy alcohol intake, regular exercise, and smoking status for model 2. The selection were a p < 0.001 for the sign test of the testing balanced accuracy (TEBA) and ten validation consistency (CVC) in the ten-fold cross-validation [22]. Among the o models to satisfy the significance level of the sign test and CVC, the model with the number of genetic variants and highest odds ratios (ORs) was selected as the o model.
The PRS for the optimal genetic model was calculated as a sum of the numbe risk alleles from each selected SNP in the optimal SNP-SNP interaction model [23-the SNP risk allele was G, the genetic scores of "AA", "AG", and "GG" were 0, 1 respectively. Each model was divided into Low-, Middle-, and High-PRS. The PR four-SNP model was classified into Low-, Middle-, and High-PRS based on scores Figure 1. Flow chart to generate polygenic risk score (PRS) associated with stature by SNP-SNP interaction and its interaction with lifestyle factors. The short stature (SS; n = 51,165; control) and tall stature (TS; n = 7536; case) groups of the participants were divided by the cutoff of ≥175 cm for men and ≥163 cm for women, which were the 15th percentiles of adult height in people aged >40 years in 2010-2014 in Korea.
Ten SNPs were included in the optimal genetic model by selecting them based on their interaction with each other from 22 genetic variants in the generalized multifactor dimensionality reduction (GMDR) software version 2.0 using the exhaustive search type after adjusting for covariates. Covariate set 1 contained age, gender, residence area, education, and income, and covariate set 2 included covariates in set 1 plus energy intake, alcohol intake, regular exercise, and smoking status for model 2. The selection criteria were a p < 0.001 for the sign test of the testing balanced accuracy (TEBA) and ten cross-validation consistency (CVC) in the ten-fold cross-validation [22]. Among the optimal models to satisfy the significance level of the sign test and CVC, the model with the lowest number of genetic variants and highest odds ratios (ORs) was selected as the optimal model. The PRS for the optimal genetic model was calculated as a sum of the number of the risk alleles from each selected SNP in the optimal SNP-SNP interaction model [23][24][25]. As the SNP risk allele was G, the genetic scores of "AA", "AG", and "GG" were 0, 1, and 2, respectively. Each model was divided into Low-, Middle-, and High-PRS. The PRS in the four-SNP model was classified into Low-, Middle-, and High-PRS based on scores of 0-3, 4-5, and ≥6, respectively; and that in the seven-SNP model was categorized into 2-7, 8-10, and ≥11, respectively. The PRS of the model was further used for analyzing the interaction with the lifestyle parameters.

Molecular Docking of Wild and Mutated GDF5 with Food Compounds
The wild and mutated chemical structures of the proteins GDF5_rs224331 (Ala276Ser) were made in the Protein Data Bank (PDB) format from the Iterative Threading Assembly Refinement (I-TASSER) website (https://zhanggroup.org/I-TASSER/ (accessed on 20 June 2022). Food compounds (n = 20,000) were downloaded from the fooDB website (https://foodb.ca/, accessed on 9 June 2022). Water molecules attached to the ligands were eliminated using the pleomorphic analysis methodology (PyMOL) software version 2.0 (DeLano Scientific LLC, South San Francisco, CA, USA) [26]. The PDB structure of the protein and ligand was changed into a protein data bank, partial charge (Q), and atom type (T) (PDBQT) files lattice format using AutoDock Tools 1.5.6 (Molecular Graphics Laboratory, Scripps Research Institute, La Jolla, CA, USA) [26]. The active sites of the proteins GDF5 were identified using the ProteinsPlus website (https://proteins.plus/, accessed on 8 July 2022), and they and the mutated sites (rs224331) were included for molecular docking.
Food components with a <−10 kcal/mol binding free energy were selected [27]. The binding free energy at the active site indicated the binding affinity of GDF5 to the food compound. The lower the binding free energy, the tighter the binding and affinity.

Molecular Dynamics Simulation (MDS)
The conformational changes in the GDF5 structures were examined using MDS to detect the changes in their activity. After the top docking poses with the selected food compounds were inserted, simulations were conducted on the GDF5 active site and docked complexes. Additionally, the Chemistry at Harvard Macromolecular Mechanics (CHARMM) force field was added to each molecular structure generated by "Simulation", and the protein was solvated by "Solvation". The "Standard Dynamics Cascade" was used to set the molecular dynamics simulation parameters for the protein added to the solvent system. The ramp-up time, equilibration time, simulation sampling time, and simulation step size were set to 40 ps, 400 ps, 10,000 ps, and 2 fs, respectively, and other parameters were set to default values. The hydrogen bond values, root mean square deviation (RMSD), and root mean square fluctuations (RMSF) were analyzed after the 10 ns simulation.

Statistical Analysis
The sample size (n = 58,701) was sufficient to show significance at α = 0.05 and β = 0.99 when assigned an odds ratio (OR) of 1.08 in the logistic regression analysis using a G-power calculator. Descriptive statistical analysis was conducted using SAS (version 9.3; SAS Institute, Cary, NC, USA). Frequency distributions were used for the categorical variables, and their statistical differences between the SS and TS groups were calculated with the Chi-square test. The adjusted means and standard deviations for the continuous variables were assessed after adjusting the covariates, including age, gender, education, income, survey year, BMI, osteoporosis, energy intake, exercise, alcohol intake, and smoking status. The statistical differences between the genders and adult height groups were analyzed using a two-way analysis of covariance (ANCOVA) with adjustment for the covariates [28]. When ANCOVA was significant, multiple comparisons according to the genders and height groups were performed using Tukey's test.
The association of height with anthropometric and biochemical parameters was analyzed using logistic regression analysis after adjustment for covariates. The results are provided with each metabolic parameter's adjusted odds ratios (ORs) and 95% confidence intervals (CI). Two different covariate sets were included according to the covariates. Covariate set 1 included age, place of residence, survey year, osteoporosis, BMI, education, and income as covariates. Covariate set 2 was calculated with covariates of model 1 plus energy intake, physical activity, smoking status, and alcohol intake.
The two-way ANCOVA was conducted with the main effect terms of PRS and lifestyles, their interaction term, and covariates. If the interaction terms showed statistical significance, each lifestyle parameter was categorized into high or low groups according to the dietary reference intake [29] or their 30th percentile. The ORs and 95% CI of the height with PRS were also assessed with logistic regression analysis in each of the high or low groups in the lifestyle-related parameters. Any significant difference in the height was determined according to the PRS groups using the χ 2 test or ANCOVA in the lifestyle-related parameters' low-or high-groups.

Demographic Characteristics and Lifestyles According to Gender and Adult Height
The age was higher in the SS group than in the TS group in both genders ( Table 1). The proportion of tall persons was higher than that of short persons in men, and it was the opposite in women. The proportion of both men and women in the TS who lacked higher education (≤middle school) was much lower than in the SS group (Table 1). However, the proportion of men with a high income was lower in the TS than in the SS group, but it was the reverse in women. Smoking status was not significantly different between the SS and TS groups in both genders ( Table 1). The proportion of the participants with regular physical exercise was small but significantly higher in the TS group than in the SS group, but only in women. Only in men, alcohol intake was higher in the TS group than in the SS group (Table 1).

Nutrient Intake According to Gender and Adult Height
Energy intake was lower in the SS than in the TS group in both genders. Regarding the intake of macronutrients, fat intake was lower in the SS group than in the TS group only in women, and the trend of carbohydrate intake was the opposite of that of the fat intake (Table 1). However, protein intake was not significantly different between the SS and TS groups for both genders. Dietary fiber intake was lower in the SS group than in the TS group only in men, and calcium and vitamin C and D intakes were lower in the SS group than in the TS group for both genders ( Table 1). Intakes of DII and flavonoids did not differ between the TS and SS groups ( Table 1). The proportion of KBD was lower in the SS group than in the TS group only in men, but that of PBD and RMD was lower in the SS group than in the TS group only in women ( Table 1). The proportion of WSD was much lower in the SS group than in the TS group in both genders. Coffee intake, but not tea intake, was lower in the SS group than in the TS group in women (Table 1).

Prevalence of MetS and Its Related Parameters According to Gender and Adult Height
As expected, height was much greater in the TS group than in the SS group in both genders. When the increase in height ceased at age 18, the weight was not significantly different between the TS and SS groups ( Table 2). In women, the BMI during the study was higher in the SS group than the TS group and was inversely associated with height (adjusted ORs = 0.908). Waist circumferences were lower in the TS group than in the SS group for both genders. The inverse association of waist circumference with height was much more significant than the association with BMI (adjusted ORs = 0.327) ( Table 2). The skeletal muscle index (SMI) and fat mass were significantly lower in the TS group than the SS group in both genders and was inversely associated with height. White blood corpuscles (WBC), which reflect the immune status, showed a trend similar to the SMI and fat mass ( Table 2). However, there was no significant difference in serum high-sensitive C-reactive protein (hs-CRP) concentrations, an inflammatory marker, between the TS and SS groups ( Table 2).

Genetic Variants Linked to Adult Height
The genetic variants affecting adult height are shown as a Manhattan plot in Supplementary Figure S1A, and those satisfying the selection criteria were located in chromosomes 2, 3, 5, 6, 12, 15, and 20. The Q-Q plot in Supplementary Figure S1B indicated that the genetic variants did not deviate from the expected values (lambda = 1.032).

SNP-SNP Interaction by GMDR
The genetic models, including 4, 7, 8, 9, and 10 genetic variants, met the selection criteria for the optimal model. The four-SNP model included ZBTB38_rs6762722, PAIP2B_rs13034890, CPZ_rs3756173, and LTBP1_rs4630744 whereas the seven SNP model contained ZBTB38_rs6762722, PAIP2B_rs13034890, GDF5_rs143384, LCORL_rs7700107, DIS3L2_rs1249260, LTBP1_rs4630744, and NCAPG_rs2074974 (Supplementary Table S1). Height was 1.29 times and 1.22 times positively associated with PRS for models 4 and 7, respectively, after adjusting for covariate groups (group 1 included age, sex, weight, education, income, and place of residence; group 2 contained covariate group 1 plus energy intake, alcohol intake, physical activity, and smoking, as shown in Figure 2). Therefore, the adjusted OR of the PRS was higher than that of the single SNP model. However, the ORs of the 4-SNP and 7-SNP models were similar, and the PRS of the 4-SNP model was appropriate. The four-SNP model was, therefore, more appropriate to explain the genetic impact on height. Nutrients 2023, 15, x FOR PEER REVIEW 12 o Figure 2. Adjusted odds ratio (ORs) and 95% confidence intervals (CI) of 4-SNP PRS and 7-SNP P models for tall stature PRS was generated with the sum of the number of risk alleles in each S and it was classified as Low-PRS, Middle-PRS, and High-PRS according to the range 0-3, 4-5, ≥6 in the four-SNP model and 0-5, 6-7, and ≥8 in the six-SNP model, respectively. Covariate group 1 included age, gender, weight, residence area, education, and income, and those of grou contained those in group 1 plus energy intake, exercise, alcohol drinking, smoking, and osteopor incidence.

Expression of Quantitative Trait Loci (eQTL) of the Selected Genes According to the Allel
The genetic expression of the variant alleles of GDF5, NCAPG, LCORL, PAIP ADAMTSL3, and IGF1R was different across different tissues. The tibial nerve showed maximum differences in the gene expression of the alleles. The gene expression NCAPG_rs2074974 decreased with the risk allele (A) compared to the non-risk allele the tibial nerve (slope = 0.14, p = 0.003), and IGFR1 also showed a similar trend in the tib nerve (slope = 0.11, p = 0.011) ( Figure 3). However, ADAMTSL3_rs1600640 showed an creased expression with the risk allele compared to the non-risk allele in the tibial ne (slope = 0.16, p = 0.0015) and the arterial appendage of the heart (slope = 0.16, p = 2.1 × 10 The risk allele of PAIP2B_rs13034890 had a lower expression than the non-risk allele the pituitary (slope = −0.29, p = 0.000083). Interestingly, the expression of GFD5_rs224 showed a significant decrement in the risk allele compared to the non-risk allele in vari tissues, especially the cortex, hippocampus, and amygdala of the brain, tibial nerve, t roid, and heart (slope = −0.12~−0.41, p = 0.0095-1.1 × 10 −8 ) (Figure 3). These results indica that the gene expressions in the tibial nerve could act as central modulators of height a that the selected genetic variants in the tibial nerve play a critical role in height.

Figure 2.
Adjusted odds ratio (ORs) and 95% confidence intervals (CI) of 4-SNP PRS and 7-SNP PRS models for tall stature PRS was generated with the sum of the number of risk alleles in each SNP, and it was classified as Low-PRS, Middle-PRS, and High-PRS according to the range 0-3, 4-5, and ≥6 in the four-SNP model and 0-5, 6-7, and ≥8 in the six-SNP model, respectively. Covariates of group 1 included age, gender, weight, residence area, education, and income, and those of group 2 contained those in group 1 plus energy intake, exercise, alcohol drinking, smoking, and osteoporosis incidence.

Expression of Quantitative Trait Loci (eQTL) of the Selected Genes According to the Alleles
The genetic expression of the variant alleles of GDF5, NCAPG, LCORL, PAIP2B, ADAMTSL3, and IGF1R was different across different tissues. The tibial nerve showed the maximum differences in the gene expression of the alleles. The gene expression of NCAPG_rs2074974 decreased with the risk allele (A) compared to the non-risk allele in the tibial nerve (slope = 0.14, p = 0.003), and IGFR1 also showed a similar trend in the tibial nerve (slope = 0.11, p = 0.011) ( Figure 3). However, ADAMTSL3_rs1600640 showed an increased expression with the risk allele compared to the non-risk allele in the tibial nerve (slope = 0.16, p = 0.0015) and the arterial appendage of the heart (slope = 0.16, p = 2.1 × 10 −7 ). The risk allele of PAIP2B_rs13034890 had a lower expression than the non-risk allele in the pituitary (slope = −0.29, p = 0.000083). Interestingly, the expression of GFD5_rs224331 showed a significant decrement in the risk allele compared to the non-risk allele in various tissues, especially the cortex, hippocampus, and amygdala of the brain, tibial nerve, thyroid, and heart (slope = −0.12~−0.41, p = 0.0095-1.1 × 10 −8 ) (Figure 3). These results indicated that the gene expressions in the tibial nerve could act as central modulators of height and that the selected genetic variants in the tibial nerve play a critical role in height.

Binding Affinity of Hydrolyzable Tannins to GDF5_rs224331
The wild and mutated GDF5_rs224331 exhibited somewhat different binding free energy to food agents (Table 4). Among them, rugosin A is shown in Figure 4 as an example. The wild type of GDF5_rs224331 had a decreased binding free energy (<−10.7 kcal/mol) to hydrolyzable tannins such as stachyurin, lambertianin B, sanguiin H6, lambertianin A, mongolicain A, casuariin, punicacortein D, rugosin A, rugosin E, valolaginic acid, rugosin D, cinnamtannin II, eugenigrandin A, rugosin A, Chinese tannin, and gemin D (Table 4). However, some hydrolyzable tannins, including rugosin A, rugosin D, rugosin E, and valolaginic acid, increased binding free energy to the mutated type, indicating that they had a lower binding affinity to mutated GDF5. The binding energy between the wild-type GDF5 protein and rugosin A with hydrogen bond in Figure 4A and the pink and green parts indicated hydrogen donor and acceptor, respectively. Their binding position and intermolecular force are shown in the two-dimensional picture in Figure 4B. The binding and interactions of binding affinity between rugosin A and the mutated type of GDF5 are

Binding Affinity of Hydrolyzable Tannins to GDF5_rs224331
The wild and mutated GDF5_rs224331 exhibited somewhat different binding free energy to food agents (Table 4). Among them, rugosin A is shown in Figure 4 as an example. The wild type of GDF5_rs224331 had a decreased binding free energy (<−10.7 kcal/mol) to hydrolyzable tannins such as stachyurin, lambertianin B, sanguiin H6, lambertianin A, mongolicain A, casuariin, punicacortein D, rugosin A, rugosin E, valolaginic acid, rugosin D, cinnamtannin II, eugenigrandin A, rugosin A, Chinese tannin, and gemin D (Table 4). However, some hydrolyzable tannins, including rugosin A, rugosin D, rugosin E, and valolaginic acid, increased binding free energy to the mutated type, indicating that they had a lower binding affinity to mutated GDF5. The binding energy between the wild-type GDF5 protein and rugosin A with hydrogen bond in Figure 4A and the pink and green parts indicated hydrogen donor and acceptor, respectively. Their binding position and intermolecular force are shown in the two-dimensional picture in Figure 4B. The binding and interactions of binding affinity between rugosin A and the mutated type of GDF5 are also shown in Figure 4C,D.  Figure 4E,F showed the root mean square deviation (RMSD) and root mean square fluctuation (RMSF) for GDF5 wild and mutated types binding to rugosin A. RMSD for GDF5 wild type binding with rugosin A was sustained close to 3 Å during 100 nanoseconds. RMSF for GDF5 wild-type binding with rugosin A also did not exceed 3 nm, except at the 580 residue index in the RMSF graph. These results suggest that rugosin A stably bound to the GDF5 wild type.      Figure 4E,F showed the root mean square deviation (RMSD) and root mean square fluctuation (RMSF) for GDF5 wild and mutated types binding to rugosin A. RMSD for GDF5 wild type binding with rugosin A was sustained close to 3 Å during 100 nanoseconds. RMSF for GDF5 wild-type binding with rugosin A also did not exceed 3 nm, except at the 580 residue index in the RMSF graph. These results suggest that rugosin A stably bound to the GDF5 wild type.

Association of PRS with the Risk of Metabolic Syndrome and Its Components
Height was greater with PRS in the ascending order. The PRS for the four-SNP model was not related to body composition, including BMI, waist circumferences (Table 5), SMI, and fat mass, except height. Osteoporosis and arthritis showed no significant change with changes in PRS. However, PRS was inversely related to WBC, an immunity index ( Table  5). The PRS was inversely linked to MetS, blood HbA1c, and serum glucose concentrations, but was not significantly associated with insulin resistance. However, lipid profiles, blood pressure, and serum ALT and AST concentrations were not associated with PRS (Table 5). Therefore, the association of height-related PRS with glucose metabolism is likely related to insulin secretion. However, serum HDL, LDL, and triglyceride concentrations did not differ among the PRS groups. SBP and DBP showed a similar trend (Table  5). However, serum ALT concentrations but not AST were lower in the High-PRS compared to the Low-PRS groups and were inversely associated with PRS (Table 5). Table 5. Adjusted means and odds ratio (ORs) of metabolic syndrome (metS)-related parameters according to polygenic risk score (PRS) of the four-SNP model generated from adult height.

Association of PRS with the Risk of Metabolic Syndrome and Its Components
Height was greater with PRS in the ascending order. The PRS for the four-SNP model was not related to body composition, including BMI, waist circumferences (Table 5), SMI, and fat mass, except height. Osteoporosis and arthritis showed no significant change with changes in PRS. However, PRS was inversely related to WBC, an immunity index ( Table 5). The PRS was inversely linked to MetS, blood HbA1c, and serum glucose concentrations, but was not significantly associated with insulin resistance. However, lipid profiles, blood pressure, and serum ALT and AST concentrations were not associated with PRS (Table 5). Therefore, the association of height-related PRS with glucose metabolism is likely related to insulin secretion. However, serum HDL, LDL, and triglyceride concentrations did not differ among the PRS groups. SBP and DBP showed a similar trend (Table 5). However, serum ALT concentrations but not AST were lower in the High-PRS compared to the Low-PRS groups and were inversely associated with PRS (Table 5). Table 5. Adjusted means and odds ratio (ORs) of metabolic syndrome (metS)-related parameters according to polygenic risk score (PRS) of the four-SNP model generated from adult height.

Interaction between PRS and Lifestyle Factors for Adult Height
Among the lifestyle factors, energy intake and rice-main diet interacted with PRS to influence height (p = 0.0078 and p = 0.0095) ( Table 6). In people with high energy intake, PRS was positively related to height by 1.414 times but not in those with low energy intake. PRS interacted with the rice-main diet to affect height among the dietary patterns (Table 6). In adults with a low rice-main diet, PRS was positively linked to height by 1.318 times but not in those with a high rice-main diet. The results suggested that the genetic impact for adult height was not shown to be related to low energy intake, mainly containing rice.

Discussion
The adult height peaked in the 20s and gradually reduced after 30 years of age. Although the loss of height varies from individual to individual, various factors are involved in the decrement. The impact of genetics explains about 10% of the height variation among people of East Asian, Hispanic, African, and South Asian ancestry [30]. An average adult loses about two centimeters in height every decade after age 30. Women tend to lose more height than men due to an accelerated rate of bone loss during menopause [31]. The loss of height can also be accelerated by certain health conditions such as osteoporosis, poor nutrition, lack of exercise, and some medical therapies such as chemotherapy and radiation therapy [32]. Therefore, the present study explored the genetic variants for tall stature using adult height and adjusting for covariates related to metabolic syndrome. Our results indicated that the selected genetic variants could explain the genetic impact on tall stature in Asians.
Short height in adults is linked to an increased risk of MetS. Previous studies have demonstrated that being short is associated with higher body fat, higher triglycerides, lower serum HDL cholesterol concentrations, and an increased risk of hypertension and type 2 diabetes, all of which are components of MetS [33,34]. However, tall stature is associated with obesity during childhood [35]. The association between stature and MetS and the impact of genetics remains unclear. The present study exhibited that tall stature was inversely linked to MetS, its components, CVD, fat and skeletal muscle mass, serum ALT and AST concentrations, and arthritis, but not osteoporosis, after adjusting height-related parameters. It was related to increased insulin resistance, which might be a primary factor for the inverse relationship between tall stature and the risk of MetS and CVD. An earlier German study demonstrated that an additional height of 10 cm is associated with a 41% and 33% lower risk of type 2 diabetes among men and women, respectively. Moreover, the association of height with a lower risk of type 2 diabetes was also seen in overweight or obese men and women (36% in men and 30% in women), although the decrease in risk was lower than that seen in normal-weight adults [33]. Therefore, increased insulin resistance in short stature elevated MetS risk, which was consistent with the results of the current study.
Common genetic variants significantly influence height. Genetic variants affect the trait by the cumulative effects of many alleles at multiple loci rather than a single genetic variant. In European ancestries, one of the most commonly identified variants is the rs1042725 variant in the HMGA2 gene, which is associated with a 0.4 cm increase in height [36,37]. However, the present study did not include HMGA2_rs1042725 as a heightrelated genetic variant in Asians, indicating that Asians have different genetic variants for height traits. In the genetic investigation of the anthropometric traits (GIANT) consortium, 12,111 common SNPs are associated with the height trait, accounting for 10-40% of all height variations depending on the person's ancestry [30]. Previous studies have shown that genetic variants for the height trait are somewhat different between Asians and Caucasians. The genetic variants in the loci of LCORL, CABLES1, CDK10, ZBTB38, ZNF638, and TSEN15 are linked to stature in Han Chinese from the Beijing study [9]. The polygenic loci of LCORL, DIS3L2, EFEMP1, ZBTB38, high mobility group AT-hook 1 (HMGA1), citrate synthase (CS), and GDF5 are also shown to be linked to stature in Taiwanese [38] and Japanese [39]. Some genetic variants related to height in Asians (ADAMTSL3, ZBTB38, LCORL, DIS3L2, and GDF5) overlapped with those in the present study. Therefore, polygenic variants linked to height vary according to ethnicity.
Height is inversely associated with insulin resistance, and the waist-to-height ratio is used to predict insulin resistance [40]. The present study also demonstrated that height and waist-to-height ratio were inversely associated with insulin resistance. Some polygenic variants related to height, such as IGF1R, GDF5, DIS3L2, and ADAMTSL3, have been reported to modulate serum glucose and insulin concentrations, potentially [41,42]. Interestingly, the PRS of the common polygenic variants related to height was associated with MetS and its components, fasting serum glucose and blood HbA1c concentrations, but not insulin resistance in the present study. The PRS was also associated with MetS, a consequence of insulin resistance. However, PRS was not associated with insulin resistance. PRS could be associated with insufficient insulin secretion, commonly seen in Asians with type 2 diabetes, and it suggests that height-related genetic variants could be associated with insulin secretion.
Genetic variants affect gene function by modifying gene expression and/or resulting in catalytic activity by altering protein conformation, especially in missense mutations. Gene expression having risk alleles is different in different tissue types. In the present study, the polygenic variants selected for stature were expressed mainly in the tibial nerves. In the tibial nerve, the risk alleles of NCAPG, LCORL, IGF1, and GDF5 had a lower expression than the non-risk allele, while the risk allele of ADAMTSL3 showed a higher expression than the non-risk allele. NCAPG is responsible for condensing and stabilizing chromosomes during mitosis and meiosis. It is also involved in the carcinogenesis and progression of tumors [43]. LCORL is a transcription factor involved in spermatogenesis related to skeletal frame size and adult height [44,45]. The NCAPG-LCORL locus is associated with body growth and feed intake in cattle [46]. The IGF1R → insulin/phosphatidylinositol-3 kinase (PI3K) → protein kinase B (Akt) signaling pathway is a critical pathway for the growth of long bones [47]. The mutation in c.926C > T of IGF1R is involved in severe short stature in Chinese [48]. However, it did not show a significant relationship to adult stature in the present study. On the other hand, the rs2871865 variant in the intron of the IGF1R was significantly associated with stature. Its expression of the risk allele was lower than that of the non-risk allele in the tibial nerve. These results suggest that the genetic variations for height were mainly expressed in the tibial nerve, which receives the message from the brain for movement of the legs, feet, and toes. The tibial nerve conduction velocity is reported to be inversely associated with height [49].
Among the selected genetic variants for height, the expression of GDF5 was the most prominent and has been revealed to be involved in height and osteoarthritis. GDF5 variants are expressed in tissues such as the tibial nerve, brain, pituitary, thyroid, and adipose tissue. Its risk allele lowered the GDF5 protein expression more than the non-risk allele [50][51][52]. In addition, GDF5_rs224331 could affect the catalytic activity since the SNP site is a missense mutation. Hydrolyzable tannins are reported to increase the expression of collagen I, a primary component of the extracellular matrix found in the bones, to improve bone growth and ameliorate osteoarthritis [53]. They also stimulate osteoblast proliferation and differentiation and suppress osteoclast activity [54,55]. Hydrolyzable tannins have been demonstrated to have anti-inflammatory and antioxidant properties, which may help in bone growth [56]. It suggests hydrolyzable tannins may affect GDF5 to modulate bone growth and osteoarthritis. The present study demonstrated that some bioactive compounds had a decreased binding free energy with the GDF5_rs224331 wild type. Hydrolyzable tannins such as stachyurin, lambertianin A and B, sanguiin H6, mongolicain A, casuariin, punicacortein D, cinnamtannin II, eugenigrandin A, Chinese tannin, and gemin D had a lower binding free energy. However, the mutated GDF5 had a decreased binding affinity (increased binding free energy) with rugosin E, rugosin A, rugosin D, and valolaginic acid. It indicated that GDF5_rs224331 might achieve its activity by modifying binding free energy. Bioactive compound effects can be altered with the GDF5_rs224331 mutation. Further experimental studies are needed to confirm this.
The study is novel in showing that polygenic variants are involved in height through SNP-SNP interactions in Asians and were also associated with immunity and glucose metabolism. Interestingly, their expressions were mainly linked to the tibial nerve. GDF5_ rs224331 is a missense mutation; its binding affinity to some hydrolyzable tannins, such as rugosin E, rugosin A, rugosin D, and valolaginic acid, was lower in the wild type than the mutated one. Although further studies are needed, it can be suggested that bone growth is related to modulating the expression and binding affinity of GDF5_rs224331. However, the limitations of the present study were as follows: First, adult height was used for estimating genetic impact, which could provide good results since it was adjusted for the covariates affecting reduction in height after age 30. Second, it was conducted as a cross-sectional study; although the sample size was large (n = 58,701), the specialists gathered the samples uniformly from the volunteers. Third, data on usual food intake were gathered using the SQFFQ, in which usual food intake could be underestimated or overestimated, although it was designed for Koreans and validated with 3-day food records for four seasons.

Conclusions
The genetic impact on tall stature was found to be 1.29 times with the four-SNP model, including ZBTB38_rs6762722, PAIP2B_rs13034890, CPZ_rs3756173, and LTBP1_rs4630744. Although the GDF_rs224331 and IGF-1R_rs2871865 were significantly associated with tall stature, they did not interact with other SNPs to lower the adjusted OR in model 7. Furthermore, the PRS was inversely associated with MetS, hyperglycemia, and WBC risk. The SNPs in the model were explicitly expressed in the tibial nerve, which is associated with increased height. Hydrolyzable tannins lowered the binding free energy with the wild type of GDF5_rs143384. However, some hydrolyzable tannins (rugosin E, rugosin A, rugosin D, and valolaginic acid) did not decrease binding free energy with the mutated gene. Therefore, the genetic variants for tall stature may modulate not only height growth but also MetS, glucose metabolism, and immunity by altering the gene expression and/or their activity. Adults with PRS for short-stature need to be more cautious of the MetS risk. These results can be used in precision nutrition after further clinical study has been conducted.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/nu15071764/s1. Table S1: Generalized multifactor dimensionality reduction (GMDR) results of multi-locus interaction with genes related to adult height; Figure S1: Distribution of genetic variants for tall stature by a genome-wide association study. Informed Consent Statement: All participants signed a written informed consent form.

Data Availability Statement:
The raw data involved in this study will be available by the corresponding author to any qualified researcher.