|Home | About | Journals | Submit | Contact Us | Français|
Recent success of genome-wide association studies (GWASs) on human height variation emphasized the effects of individual loci or genes. In this study, we used a developed pathway-based approach to further test biological pathways for potential association with stature, by examining ~370 000 single-nucleotide polymorphisms (SNPs) across the human genome in 618 unrelated elder Han Chinese. A total of 626 biological pathways annotated by any of the three major public pathway databases (KEGG, BioCarta and Ambion GeneAssist Pathway Atlas) were tested. The regulation-of-autophagy (ROA) (nominal P=0.012) pathway was marginally significantly associated with human stature after our family wise error rate multiple-testing correction. We also used 1000 random recruited US whites for further replication. Interestingly, the ROA pathway presented the strongest signals in whites for height variation (nominal P=0.002). The results correspond to biological roles of the ROA pathway in human long bone development and growth. Our findings also implied that multiple-genetic factors may work jointly as a functional unit (pathway), and the traditional GWASs could have missed important genetic information imbedded in those less significant markers.
Human stature has a long research history in the community, as it is the most visible traits with remarkable high heritability and can be easily and accurately measured.1 Specifically in Han Chinese, heritability of stature is > 0.6.2 The major driver of those studies on stature is to provide clues for human growth, bone development and onset of some diseases such as cancer, type 2 diabetes, H syndrome and osteoporosis.1,3,4
Recent success in genome-wide association studies (GWASs) has revealed that some pathways can influence the human height.5–10 Strong statistical and experimental evidence showed that pathways underlying bone and cartilage development have critical roles in the formation of adult height,1 such as pathways of Hedgehog signaling,6 bone morphogenic9,11 and extracellular matrix.1 Meanwhile, some pathways may have indirect impact on human height and highlight more fundamental biological processes, including the pathways of chromatin structure regulation,6 cell-cycle regulation1 and the C-type natriuretic peptide signaling.12
Our earlier GWAS performed in Han Chinese13 pinpointed some single-nucleotide polymorphisms (SNPs) that was influencing Chinese height variation. Nevertheless, the most significant SNPs can only explain limited variance (~0.3%) in stature.13 The pathway-based association test integrates current knowledge about functional interactions among genes and represents a promising approach to look into genome-wide data to find the potential joint effects or functional interactions of genes, which were imbedded in the remaining unreported SNPs that did not achieve the genome-wide significant level.14
In this study, the genome-wide pathway-based analysis was conducted to detect the potential association between 626 pathways and height variation in Han Chinese. We aimed to identify pathways that can be influencing Chinese stature and/or to provide novel hypotheses in the mechanisms of human growth. This study may provide new insights into height variation, especially for the Chinese ethnic, which was less studied and genetically differs from the whites.1,9,13
The study was approved by the Research Administration of Xi’an Jiaotong University to assure that the general ethnic rules were addressed. All study volunteers have signed the informed-consent forms before they entered the study. A total of 618 elder Han adults were randomly recruited from the Xi’an City and its surrounding areas in northern China. Basic characteristics of the study sample are given in Table 1a. Stature for each subject was measured using a stadiometer wearing no shoes. The daily calibrated stadiometer was in well maintenance and the long-term precision was monitored. The coefficient of variation of stature measurements was 5.8%. Strict data quality control was conducted by excluding the outliers and database inputting errors. The whole data set was double crosschecked by at least two different researchers.
All the individuals with chronic diseases and conditions that might potentially affect hormone level, bone structure or metabolism were excluded from this study, which included but not limited to chronic disorders involving vital organs, diabetes, hyperthyroidism, hypo-/hyper-parathyroidism and other serious skeletal diseases (Paget’s disease, osteogenesis imperfecta, rheumatoid arthritis, kyphosis, serious vertebral fracture), chronic use of drugs affecting bone metabolism (hormone replacement therapy, corticosteroid therapy, anticonvulsant drugs, anti-bone-resorptive or bone anabolic agents) and malnutrition conditions (such as anorexia nervosa).
Genomic DNA was extracted from whole human blood using a commercial isolation kit (Gentra systems, Minneapolis, MN, USA) following the standard protocols. For each sample, genotyping with GeneChip Human Mapping 500 K set containing 250 K Nsp array and 250 K Sty array (Affymetrix, Santa Clara, CA, USA) was performed using the standard protocols recommended by the manufacturer.
From the initial full-set of 500 568 SNPs, we discarded 18 067 SNPs with sample call rate <90%, 13 041 SNPs with allele frequencies extremely deviating from Hardy–Weinberg equilibrium (P<10−7) and 98 202 SNPs with minor allele frequencies <1%. The final SNP number for the analyses is 371 258. This SNP set covered 14 642 genes and yielded an average marker spacing of ~7.9 kb throughout the human genome.
We used STRUCTURE 2.2 program (http://pritch.bsd.uchicago.edu/software.html) to detect the potential population stratification of our sample. The program used a Markov chain Monte Carlo algorithm to cluster individuals into cryptic subpopulations on the basis of multi-loci genotype data.15 We performed the analysis assuming the number of population strata k=2 and using a set of 1000 unlinked SNPs randomly selected genome wide. EIGENSTRAT16 was used to crosscheck the stratification in our analyses.
To identify the pathways potentially contributing to Chinese stature variation, we comprehensively retrieved the pathways collected by three major pathway databases, KEGG (http://www.genome.ad.jp/kegg/pathway.html), BioCarta (http://www.biocarta.com/genes/index.asp) and Ambion GeneAssist Pathway Atlas (http://www.ambion.com/tools/pathway). The number of pathways finally retrieved from BioCarta, KEGG and Ambion were 364, 262 and 201, respectively. To relieve the multiple-testing issue and avoid testing over narrow or broad functional pathways, we only included those pathways with 10–200 genes and at least 85% of which were covered by at least one SNP of our qualified genotypic data. The pathways with the following characteristics were filtered: (1) the pathway contains < 10 or > 200 genes in the above three database; (2) after discarding those unqualified SNPs (bad sample call rate, no Hardy–Weinberg equilibrium, small minor allele frequencies), the pathway contains < 10 or > 200 genes in our genotypic data; (3) after filtering those SNPs that were 500 kb away from the involved gene, the pathway contains < 10 or > 200 genes in our cleaned data and (4) the pathway contains genes within 10–200 gene, but <85% of them were included in our analysis due to discarded SNPs or SNPs with >500kb distance. Overall, a total of 626 pathways/functional networks were finally analyzed.
Gender and age are significant fixed effectors on stature13 and they were used to adjust the raw stature data. The GWAS analysis for individual SNPs was achieved by the Wald t-test implemented in Plink (version 1.03).17 The used pathway-based GWAS approach was proposed and developed first by Wang et al.14
SNPs that are 500 kb away from any gene are not considered. The SNPs of our identified pathway were verified without overlapping or misrepresentation problem (Figure 2). Among all the SNPs of a given gene, the SNP achieving the smallest P-value (the highest Wald t-test statistic values) in the single-marker association tests was assigned to represent the magnitude of association evidence of the gene. If in very rare cases where one SNP is located within shared regions of two overlapping genes, we map the SNP to both genes. They will all get the same score. We ranked all genes by sorting their statistic values from the largest to smallest. Then, a weighted Kolmogorov–Smirnov-like running-sum enrichment score (ES) was calculated.14 Briefly, ES score was used to reflect the enrichment of a given gene set at the top of the entire genome-wide ranked gene list. It measures the deviation from zero to maximum level of the highest statistic value encountered in the random walk. High ES score means this set of genes (genes of pathway) concentrated at the top of ranked gene list.
The nominal significance level of ES for each pathway was estimated by shuffling phenotypic data. The procedures of ES calculation were repeated for 1000 times. Finally, an empirical distribution of ES (ESnull) was generated for each pathway. The nominal significance of an observed ES for a pathway (nominal P-value) was estimated as the percentage of permutations whose ESnull>observed ES.
Finally, two measures were used to adjust for multiple-hypothesis testing for more reliable results. We first normalized the ES to correct the influence of gene size, so that different gene sets are directly comparable to each other. The false-discovery rate (FDR) Q value (denoted as QFDR) and family wise error rate (FWER) P-value (denoted as PFWER) were generated for correcting the multiple testing.14 FDR is a procedure that is usually used to control the fraction of expected false-positive findings to stay below a certain threshold. FWER is another approach that is considered to be highly conservative seeking to ensure that the list of reported results does not include even a single false-positive gene set.
To identify the common pathway(s) that may influence height in distinct ethnics, and also as a replication effort, we compared the results of this study with those of the pathway-based GWAS analysis for height variation in our American white sample study, which was approved by the University of Missouri-Kansas City Institutional Review Board. A total of 1000 whites were randomly recruited from Midwestern of United States (Table 1b). The same genotyping and statistical approaches were applied across the Chinese and the American studies.5
A total of 618 individuals were included in the analysis. Basic characteristics of the study sample are shown in Table 1a.
The STRUCTURE program results showed that all the subjects were clustered together as a single group by using 1000 unlinked SNPs, which suggested that no significant population stratification was present in our sample. The analyses by EIGENSTRAT confirmed that only one principal component was significant (P<0.001), which suggested that our study sample originated from one common ancestry population. Consequently, the inflation factor λ calculated by the Genomic Control approach18 implemented in EIGENSTRAT is 1.02, close to 1 standing for no stratification.
Out of the 626 pathways, the regulation-of-autophagy (ROA) pathway achieved the strongest evidence of association with height variation (nominal P=0.012, QFDR=0.30 and PFWER=0.11). The ROA pathway was involved in 35 genes (data from KEGG database), with 28 genes covered by our SNP set (Table 2). Among the 28 genes, 5 genes (PIK3C3, ULK2, ATG5, GABARAPL1 and IFNG) got nominal P<0.01, and 14 genes were ranked among the top 5000 (leading edge) genes in the single-SNP GWAS. The ROA pathway achieved a high ES score of 0.47 with nominal P-value=0.012 (Figure 1). The individual gene with the most significant P-value among all the 28 genes is PIK3C3 (nominal P=0.0009). Other ROA genes adding to the ES score (that is the genes that ranked before or at the point of maximize ES achieved, also denoted as leading-edge genes) are listed in Table 2 and Figure 1. Information of the representative SNPs for the 17 leading-edge genes that contributed positively to the ES is given in Table 3.
In the American whites, the ROA pathway was the only one that consistently presented the strongest association signals with height variation (nominal P=0.002, ES score=0.55, QFDR=0.042, PFWER=0.016). The individual gene with the most significant nominal P-value is IFNA4 (nominal P=0.0026). Seven genes of the ROA pathway (IFNA4, IFNA21, PIK3C3, ULK2, IFNA14, IFNA13 and IFNG) closely associated with height in our GWAS with nominal P-value <0.01. The genes PIK3C3, ULK2 and IFNG presented significant nominal P-value <0.01 in both Han Chinese and American whites. The results from United States did not indicate other novel pathways or replicate the prior known ones regulating human height (all other tested pathways had QFDR>0.2 and PFWER>0.13).
We also compared the P-values of the same representative SNPs in Chinese with in US whites, and none of them were significant (except for gene PRKAA2) (Table 4). Thus, different sample may have different representative SNPs even for the same gene, because of the different minor allele frequency and the distinct ethnic background.
The pathway-based approach can work as a complementary strategy to traditional SNP-based analyses and could be helpful for the detection of novel mechanisms underlying complex diseases/phenotypes. In this study, this approach was used to analyze our genomic genotyping data for height variation in an elder Han Chinese population.13 The ROA pathway was suggested at the genome-wide level pathway-based approach after the multiple-testing correction. The result implied that autophagy-related biological processes may have important roles in human stature variation, and thus may participate in the development and growth disorders of bones.
The ROA pathway (Pathway ID hsa04140 in KEGG database) is the unique pathway related to autophagy among all the pathways tested. Autophagy is a conservative house-keeping catabolic process that involves degradation of a cell’s own components by formation and separation of a membrane around the targeted region of the cell. The resultant vesicle then fuses with a lysosome and then the inner contents degrade. The autophagy process was generally induced and enhanced by nutrient starvation, defense of infection or cell programmed death.
Involvement of the ROA pathway in human bone and cartilage growth has been demonstrated in earlier studies. Long bone growth and consequent increase in height take place through the activities of chondrocytes embedded in the epiphyseal growth plate. The embedded chondrocytes of long bones undergo autophagy and express autophagic proteins before the induction of osteogenesis.19 This is because autophagy enhances the survival of chondrocytes in the hypoxic microenvironments (the cartilage is avascular, thus contains little nutrition and oxygen) and induces the synthesis of the calcified extracellular matrix.20–22 In addition, the crosstalk of bone growth factors such as interleukin-3 with the ROA pathway regulates the utilization of exogenous nutrients to maintain cell viability.23 As a result, autophagy is an important intermediate stage during the terminal differentiation of chondrocytes and long bone development. Notably, the interferon-α family genes in the ROA pathway have some other effects on bone remodeling (Table 2). They can broadly suppress the proliferation of osteoprogenitor cells and inhibit the differentiation of monocytes to osteoclasts.24,25
Our study for the first time raises the importance of the ROA pathway to Chinese stature variation. The pathway was also highly suggested by our pathway-based genome-wide analysis that could influence the Caucasian height variation. In addition to highlighting its potential contribution to human height variation, the ROA pathway was demonstrated that associated with long bone mineral density at the ultra distal skeletal site (our unpublished data). Our findings in the bone mineral density study indicated that the ROA pathway may influence the height variation through different rate of long bone development and remodeling. Further functional researches are warranted.
A major strength of this study may be embodied by raised importance to height of the ROA pathway, which are largely neglected by previous GWASs.1 This may partially attribute to the pathway-based approach applied here that allows to spontaneously examine the genetic information imbedded in multiple SNPs and genes. This nature of the pathway-based method alleviates the punishment for the multiple-testing problem, and allows examining the combined effects other than individual genes whose effects are too modest to be detected through single-SNP tests. Given the complexity of biological systems, such circumstance should not be rare. Therefore, we argue that the pathway-based tests can serve as a useful complement of traditional GWAS. In addition, we emphasized that as the pathway analysis aims to jointly consider multiple variants and measures the enrichment of gene set at the top of the entire ranked gene list, only a few dramatically significant SNPs in GWAS study (for example SNPs in IHH gene of Hedgehog Signaling) may not necessarily pull the entire pathway to the top. This phenomenon cannot be defined as the ‘non-replication’ for earlier studies.
However, several concerns need to be clarified in our study. First, the average age is over 70 years for Chinese and >50 years for Whites. To use stature data of aged group cannot represent a typical adult height variation. As the bone remodeling is active in the elder people, we may actually detect a pathway that enhanced the bone absorbing process. Some famous validated pathways (for example Hedgehog Signaling) in whites were not found by our study, one of the potential confounding factors could be the sample age bias. The ethnic difference was also another factor for non-replication because those pathways were initially identified at the Caucasian ethnic context.1
Second, our sample size is small. Large sample size can dramatically improve the statistical power for detecting the potential pathways. The association signal of ROA in our study (618 individuals) was weak (nominal P=0.012) compared with our US replication study (nominal P=0.002, 1000 individuals). Therefore, we reported the ROA pathway based on the statistical evidence with a consistently increased trend.
Third, although pathway-based analysis is a promising strategy to data mining of GWAS data sets, it is currently unclear what the best strategy is. Several methods are developed for pathway hunting recently.26–30 All of these methods used relatively complex model or simulation to improve the accuracy and statistical power. In fact, how to properly account for the gene size, gene number and linkage disequilibrium (LD) and how to adopt the optimal SNP weighting scheme (for example the representative SNPs located in the coding sequence could be more important than others) are worth further investigating. We used this method because it was generally robust to gene size, LD, gene number and, thus, accepted by the community.14,31,32 In addition, it was the real ‘hypothesis-free’ approach for pathway investigation, which means it can comprehensively test the association across different up-to-date pathway databases without prior hypotheses. Some studies first identified promising genes in one hypothesized pathway, then assessed the association of those genes with specific traits such as Parkinson disease.33,34 However, similar approaches were limited for quantitative complex traits (height), because the unknown fundamental biological processes may influence the phenotypes, such as the EphrinA–EphR pathway for bone geometry,31 the C-type natriuretic peptide signaling for human height.12
Finally, there are eight genes (eight representative SNPs) of ROA pathway cluster in a small region as shown in Figure 2 and Table 4. In addition, we constructed the LD block map for them (Figure 3), and six of them were located in one strong LD block. Although the permutation approaches may relieve this issue,14,32 the cluster of SNPs could still influence the result, and we lack a reliable method to accurately evaluate the extent of the influence. In our recent study, we performed the pathway-based study through shuffling the genotypic data (not phenotypic as the original).35 This procedure can completely eliminate the influence of strong LD, because it breaks the structure of original LD. However, the simulated results showed that the false-positive rate was not reduced; rather, it was slightly larger than the one by shuffling phenotypic data. Therefore, the influence of strong LD actually is far more less than one can expect. To permutate, the phenotypic data were a more conservative pathway-based association approach. To report this result, we mainly based on the statistical evidence and the independent replication study. Therefore, our study may only be served as a tentative investigation using this kind of test strategy. The future studies with larger sample size are necessary to validate our result.
In summary, the ROA pathway was suggestively associated with the height variation in Han Chinese in our pathway-based GWAS analyses. The genetic findings were supported by their biological implications in long bone development and replication study in the whites. Further functional researches are still warranted to fill the biological mechanism details.
We thank two anonymous reviewers for their time and effort in helping us improve and clarify this manuscript. HWD was partially supported by grants from NIH (P50AR055081, R01AG026564, R01AR050496, RC2DE020756, R01AR057049 and R03TW008221) and Franklin D Dickson/Missouri Endowment. The study also benefited from grants from National Science Foundation of China, Huo Ying Dong Education Foundation, HuNan Province, Xi’an Jiaotong University and the Ministry of Education of China. This study was supported by NIH grants from R01 AR050496, R21 AG027110, R01 AG026564, P50 AR055081 and R21 AA015973.
CONFLICT OF INTEREST
The authors declare no conflict of interest.