Whole Exome Sequencing Uncovered the Genetic Architecture of Growth Hormone Deficiency Patients

Purpose Congenital growth hormone deficiency (GHD) is a rare and etiologically heterogeneous disease. We aim to screen disease-causing mutations of GHD in a relatively sizable cohort and discover underlying mechanisms via a candidate gene-based mutational burden analysis. Methods We retrospectively analyzed 109 short stature patients associated with hormone deficiency. All patients were classified into two groups: Group I (n=45) with definitive GHD and Group II (n=64) with possible GHD. We analyzed correlation consistency between clinical criteria and molecular findings by whole exome sequencing (WES) in two groups. The patients without a molecular diagnosis (n=90) were compared with 942 in-house controls for the mutational burden of rare mutations in 259 genes biologically related with the GH axis. Results In 19 patients with molecular diagnosis, we found 5 possible GHD patients received known molecular diagnosis associated with GHD (NF1 [c.2329T>A, c.7131C>G], GHRHR [c.731G>A], STAT5B [c.1102delC], HRAS [c.187_207dup]). By mutational burden analysis of predicted deleterious variants in 90 patients without molecular diagnosis, we found that POLR3A (p = 0.005), SUFU (p = 0.006), LHX3 (p = 0.021) and CREB3L4 (p = 0.040) represented top genes enriched in GHD patients. Conclusion Our study revealed the discrepancies between the laboratory testing and molecular diagnosis of GHD. These differences should be considered when for an accurate diagnosis of GHD. We also identified four candidate genes that might be associated with GHD.


INTRODUCTION
Congenital growth hormone deficiency (GHD) is a rare disease characterized by decreased growth hormone (GH) secretion of the anterior pituitary, which leads to growth impairment and metabolic dysfunction in children (1, 2). The causes of GHD include pituitary dysplasia and pathogenic mutations in GHinsulin like growth factor 1 (IGF1) axis-related genes, such as GH1, GHRHR (3,4). However, recent cohort studies of the genetic causes of GHD patients in South and East Asian populations reported that the proportion of patients with molecular diagnosis ranges from 4% to 43% (5,6). The undiagnosed rate indicates that there are still a lot of unknowns about the genetic predispositions of GHD. It is challenging to diagnose GHD only based on clinical symptoms. Clinical presentation of symptoms including neonatal hypoglycemia, midfacial defects such as cleft lip or palate, history of external head injuries, and vacuolated Sella turcica found in pituitary magnetic resonance imaging (MRI) or pituitary dysplasia can assist in GHD diagnosis (7)(8)(9). According to the guidelines set by the GH Research Society (GRS), the clinical diagnosis of GHD needs to be based on growth and development data, and data from various laboratory tests (i.e. decreased peak of GH provocation test with two different agents, serological detection of IGF1 and IGF1BP3, combined gonadal hormone when necessarily, and genetic test), and imaging data (craniocerebral MRI) (10). However, there is still a high falsepositive rate in the existing serological diagnosis methods (11)(12)(13)(14).
With the application of whole exome sequencing (WES) in the exploring the etiology of GHD, an increasing number of GHD-associated genes have been discovered (9). It has been recently reported that a variety of congenital diseases associated with GHD are caused by certain genetic variants. For example, choreoathetosis, hypothyroidism, and neonatal respiratory distress caused by NKX2-1 (MIM: 610978), neurofibromatosis, type 1 was caused by NF1 (MIM: 162200), growth plate disorders caused by FLNB (15)(16)(17). To further explore the contribution of genetic variations to GHD and comprehensively decipher the molecular basis of GHD at an exome level, we herein investigated the difference between molecular findings and clinical diagnosis criteria by whole exome sequencing (WES) data in a cohort of 109 possible GHD patients. To elucidate the genetic architecture of GHD patients, we further analyzed the mutational spectrum in patients without a molecular diagnosis by a mutational burden analysis.

Cohort Collection
From July 2014 to August 2018, we screened 561 patients with short stature from three centers in China [Maternal and Child Health Hospital of Guangxi Zhuang Autonomous Region, the Second Affiliated Hospital of Guangxi Medical University and Beijing Children's Hospital, all three centers as parts of the DISCO study (http://www.discostudy.org/)]. Patients with pituitary tumors, chronic diseases, iatrogenic shortness, and previous GH treatment were excluded from this cohort. Demographic information, physical examination results and clinical symptoms were taken into account. Details of cohort were reported in the previous study (18,19).
Radiographic evaluation included conventional radiographs of the anterior-posterior view of the left hand and pituitary MRI. Radiographs of the left hand and wrist for bone age were assessed independently by two pediatric clinicians. GH stimulation tests with L-dopa and insulin were performed. L-dopa was administered orally, and insulin was subcutaneously injected in the overnight fasting state. Blood samples were collected at 0, 30, 60, 90, and 120 min, respectively. GH concentration of each time was measured using a chemiluminescence method with a sensitivity of 0.01 ng/ml. Serum IGF1 was measured by the chemiluminescence immunometric method. Two GH stimulation tests were provided for all enrolled patients. This study was obtained approval from the ethics committee at the Maternal and Child Health Hospital of Guangxi Zhuang Autonomous Region (G-1-1), the Second Affiliated Hospital of Guangxi Medical University (2020-KY[0112]) and Beijing Children's Hospital (Y-028-A-01).
Based on GHD diagnostic criteria of other countries and reference data of serology in the East Asian populations (20,21), we recruited 109 unrelated patients with potential GHD using following eligibility criteria: 1) age of male patient < 12 yrs., female patient age < 10 yrs., 2) height standard deviation < -2 standard deviation (SDs), 3) the peak GH concentration < 7ng/ml in one stimulation test, 4) the standard deviation of IGF1 serum level < 0 SDs, 5) the anterior-posterior plain radiograph of left wrist and hand showed delayed bone age, 6) growth velocity < 5cm/yr, 7) availability of MRI examination of the pituitary gland. According to the latest diagnostic criteria of GRS in 2019 (10), all enrolled patients were grouped according to the following items: (1) the peak GH concentration < 7ng/ml in both provocation tests, (2) the standard deviation of IGF1 serum level < -2 SDs, (3) the delayed bone age was assessed independently by two different doctors. Recruited patients were further divided into two groups: Group I (n=45) included patients considered to have definitive GHD by currently stringent diagnosis; Group II (n=64) included patients with possible GHD (Figure 1).

Whole-Exome Sequencing
DNA samples from three centers were processed and prepared under the same protocol from patient blood. In total, 109 patients underwent proband-only WES, while 9 underwent

Variant Calling and Interpretation
The sequencing data were analyzed and annotated using an inhouse developed analytical pipeline, the Peking Union Medical College hospital Pipeline (PUMP) (

Variant Validation and Sanger Sequencing
All pathogenic and likely pathogenic variants were manually reviewed using Integrative Genomics Viewer (IGV) (26). Sanger sequencing was performed independently on available subjects and parental samples to validate variant interpretation by an orthogonal sequencing method and to investigate segregation according to Mendelian expectations for the identified variant allele(s).

Mutational Burden Analysis
To extend the spectrum of known mutations in genes contributing to GH synthesis and secretion, we performed mutational burden analysis to search for rare variants of uncertain significance (VUS) in 259 candidate genes in the WES data. We first selected 259 candidate genes related to the synthesis and secretion of GH according to the Human Phenotype Ontology (HPO, https://hpo. jax.org/), the OMIM (http://omim.org/) and the Kyoto Encyclopedia of Genes and Genomes (KEGG,https://www. genome.jp/kegg/) databases ( Table S1). The variants filtering criteria was based on dominant model and recessive model according to the allele mutation site. The filtering criteria of the dominant model for predicted pathogenic mutations include: 1) public database frequency was required to be < 10 -4 ; 2) The missense mutations with Combined Annotation Dependent Depletion (CADD) score > 20 were considered as high credibility; 3) truncating mutations; 4) intra-frame mutations require a range of influence of more than one amino acid; 5) splicing mutations include canonical splicing mutations or mutations with a spliceAI predictions score more than 0.5 (27). For the recessive model, we required that the public database frequency be less than 10 -3 , and the rest of the criteria are the same as the dominant model.

Statistical Analysis
Mean comparison of relevant features was conducted using Student's t-test. One-tailed p-values of < 0.05 were considered statistically significant. Fisher's exact test and chi-square test were used for genetic burden analysis.

Cohort Information and Clinical Characteristics
A total of 109 unrelated patients were enrolled from previous cohort, including 9 cases with un-affected parents (trios) and 100 singletons (total 127 subjects). Of the 109 patients, 87 were male and 22 were female, with a mean age of 7.6 ± 3.0 years. According to the diagnostic criteria, all patients were divided into two groups. Group I consisted of 44 patients (8 female and 36 male, average age 8.6 ± 2.7) and group II of 65 patients (14 female and 51 male, average age 7.0 ± 3.0) ( Table 1). There were no significant differences in sex-distribution and delayed bone age between Group I and II except for age, Z-score of height, peak GH concentration, IGF1 and growth velocity ( Table 1).

Spectrum of Causal Gene and Variants
After WES and variant interpretation based on ACMG guidelines (25), causal variants in 14 genes were identified in 19/109 (17.4%) of the patients. Of 14 genes, seven genes were associated with GH secretion and synthesis, and 10 variants in these genes were identified in patients from both groups ( Table 2). Interestingly, DISCO-S693 of group I and DISCO-S255 of group II harbored distinct variants in GHRHR and had different IGF1 serum levels (-1.61 SDs for DISCO-S255, -2 SDs for DISCO-S693), with only DISCO-S693 meeting stringent diagnostic criteria ( Table 2). Therefore, we speculated that the  adoption of strict clinical diagnostic criteria, particularly IGF1 being less than -2 SDs, could lead to a certain false-negative rate in the diagnosis of GHD. Interestingly, seven genes that were not previously associated with GHD were identified in 9 patients from both groups. Same variants (c.253G>A[p.Ala85Thr]) in IDS were identified in 2 patients, patient DISCO-S189 from Group I and patient DISCO-S032 from Group II ( Table 2). Seven patients carried pathogenic or likely pathogenic variants in 7 genes which had no clear evidence supporting their biological function in GH synthesis or secretion (Table 2, Figure S1). Our results showed that patients with GHD might have defects in gene expression relating to diverse biological functions.

Genetic Burden Analysis Revealed Novel Genes Biological Related to GH Axis
Due to the discrepancy between the clinical diagnosis and molecular diagnosis in GHD patients, the adoption of less stringent inclusion criteria might help explore the molecular mechanism of GHD more comprehensively. To this end, we compared the mutational burden of rare coding variants of each candidate gene between the patients without a molecular diagnosis and the control samples. Among all the genes related to GH synthesis and secretion, we found that rare VUS in 40 genes enriched in our GHD cohort (Table S2). Among 40 genes enriched in our cohort, four genes were found with trend toward significance: POLR3A (p = 0.005), SUFU (p = 0.006), LHX3 (p = 0.021), CREB3L4 (p = 0.040) (Tables 3 and S2).
POLR3A encodes the largest subunit of RNA transcriptase III. Bi-allelic mutation in this gene is known to cause two hereditary neurological diseases, leukodystrophy, hypomyelinating 7 (4H Leukodystrophy, MIM: #607694) and Wiedemann-Rautenstrauch Syndrome (MIM: #264090). Through mutation burden analysis, four rare heterozygous and in silico predicted deleterious variants in POLR3A were significantly enriched in GHD patients without molecular diagnosis, and these variants were located in the conserved sites based on GERP++ scores (28) ( Table 4). In addition to the peak concentration of GH in provocation test < 7ng/ml and serum IGF1 lower than 0 SD, three of the four patients showed abnormal pituitary MRIs. DISCO-S122 harbored c.200G>A(p.Arg67His) mutation with 29 CADD scores, presenting with a small pituitary and no neurohypophysis, and accompanied by small penis. Pituitary MRI found Sella turcica dysplasia accompanied by a cerebrospinal fluid hernia, small pituitary morphology and lip cleft in DISCO-S227 patient with c.1721G>T(p.Gly574Val) mutation. The CADD score for this mutation was 33. Pituitary MRI of DISCO-S596 patient with c.2672G>A(p. Arg891Gln) mutation, of which CADD score was 36, showed pituitary stalk interruption syndrome, as well as combined congenital hypothyroidism and polydipsia polyuria. The patient DISCO-S059 had a c.1676T>G(p. Phe559Cys) mutation with 24.2 CADD score, and showed a GH provocation test peak value of below 7ng/ml and decreased serum IGF1 level ( Table 4).
SUFU encodes a negative regulatory factor in the sonic hedgehog pathway. Germline and somatic mutation of this gene can result in developmental abnormalities and tumor predisposition syndrome. In 90 GHD patients without molecular diagnosis, two patients carried rare heterozygous variants, namely, c.373A>C(p. Lys125Gln) in patient DISCO-S150 and c.868A>G(p. Ser290Gly) in patient DISCO-S041 ( Table 4). The CADD scores for these two variants were slightly higher than 21, and both variants replaced conserved residues of the respective protein.
LHX3 encodes a LIM domain transcription factor, which is involved in the early steps of pituitary ontogenesis. It has been reported that mutation of this gene causes combined pituitary hormone deficiency (MIM: #221750) with impaired production of GH and one or more of the other five anterior pituitary hormones through a recessive mode (29). We found two predicted deleterious LHX3 missense variants in two patients. The c.1159C>A(p.Pro387Thr) in LHX3 was identified in DISCOS-085, with 25 CADD score ( Table 4). This male patient was 2 years and 10 months old, presenting with severe growth retardation. The serological examination showed a peak value of GH at 5.72 ng/ml, and IGF1 was -0.39 SDs, while imaging examination found that his bone age was delayed by 10 months and pituitary MRI had no obvious abnormality. The other patient, DISCO-S311 was an 11-year old male with short stature, and carried a heterozygous missense variants (c.611G>T [p.Arg204Leu]) in LHX3 with high CADD score of 28.8. Serological examination showed that the peak concentration of GH were 6.22 ng/ml and 6.54 ng/ml, IGF1 was -1.11 SDs lower than normal, imaging examination showed that the bone age was delayed by 4 years and the pituitary MRI had no obvious abnormality (Table 4).
CREB3L4, located in the region of chromosome 1q21.3, codes a member of cAMP response element-binding (CREB) proteins family. The encoded protein contains 1 bZIP domain and  Table 4). Taken together, these four genes represent potential novel genes for GHD patients.

DISCUSSION
GHD is a rare condition characterized by the limitation of growth and development in children caused by defects in GH synthesis and secretion. However, the methods for diagnosing GHD are based mainly on clinical examinations rather than etiologically findings. In this study, we recruited 109 possible GHD patients to evaluate the consistency between clinical diagnosis and WES-based molecular testing, and how molecular testing can compensate for limitations in the clinical diagnosis of GHD. We also deciphered the molecular architecture of GHD by analyzing variants in genes biologically related with the GH axis. Due to the wide variety of clinical symptoms of GHD and the pulsatile mode of GH secretion, there is no standard for the diagnosis of GHD (3). As a result, a certain false positive rate exists in a variety of diagnostic methods, even if combined with the measurement of serum IGF1 levels and the GH provocation test (11,13,30). A previous study of molecular screening of GHD in 80 Morocco patients reported that 10% patients had monogenic defect, and serologic results of 2 patients with molecular finding in this study did not meet current diagnosed criteria (31). In our cohort study, five patients with pathogenic mutations in GH-IGF1 axis-related genes did not show a phenotype that fully met the stringent clinical diagnostic criteria of GHD. This indicates that genetic testing in GHD patients is a useful tool for clinical diagnosis and genetic counseling, which can assist clinical diagnosis.
At the same time, our results indicate that results from serologic laboratory testing still have a certain possibility of false negatives. Especially patients with mutations in the same gene cannot be fully diagnosed GHD, such as IDS and GHRHR. Recent genotype-phenotype correlation study in 24 patients with GHRHR mutation found that the different mutation types could not fully explain the phenotypic discrepancy (32), which is consistent with our findings. In addition, the inconsistent severity of phenotypes was also observed in two IDS patients with the same mutation. This might be caused by variable penetrance, where the same variant in a certain gene might give rise to different manifestations and severity of the disease. The underlying mechanism of the variable penetrance might be the genetic background and modifiers, which need to be further studied using larger sample size of phenotype-genotype correlation analysis. However, these results indicates that the stringent clinical diagnostic methods have a considerable false negative rate. Therefore, we propose that precise genetic testing may be helpful to further understanding the molecular etiology and mechanism of GHD development. Although the underlying pathogenesis of GHD developmental defects is associated with pituitary dysplasia, Sella dysplasia and GH-IGF1 axis gene mutations, most genetic causes of GHD have not been identified, which may lead a certain false-negative rate in diagnostic approach. Recently, it has been reported that some patients with Potocki-Lupski syndrome or Prader-Willi syndrome have similar phenotypes as GHD patients (33,34), and recessive mutations in RNPC3 could also lead to the occurrence and development of GHD (35). In our study, we found seven genes (PTPN11, ROR2, ACAN, COL11A1, TRPS1, IDS, MMP13) might associate with dysfunction of GH-IGF axis, indicating GHD may exist in a variety of Mendelian syndromes. Pathogenic mutations in these seven genes also might be incidental finding, which need more evidence for validation.
To further elucidate the mutational spectrum of GHD patients, we performed gene-based mutation burden analysis and identified predicted deleterious variants in 4 genes (POLR3A, SUFU, LHX3 and CREB3L4) that were enriched in GHD patients, implicating that these four genes are likely novel effectors for GHD. However, the small sample size limited the power of detecting GHD-associated genes. We choose the variants in the top four genes with p<0.05 as the statistical trend. Recently, a multicenter retrospective study suggested that patients with 4H leukodystrophy (hypomyelination, hypodontia, and hypogonadotropic hypogonadism) may have more obvious GH secretion reduction and deficiency of various anterior pituitary hormones, such as FSH and LH (36). Our results suggest that rare heterozygous variants in POLR3A may lead to different degrees of recessive phenotype in patients through dose-effect models. Due to the 'gate keeper' role of SUFU in sonic hedgehog pathway, mutations in the SUFU gene were mostly studied in the context of susceptible tumor diseases (37). As far as we know, this is the first time that this gene is reported to be related to GHD. Previous reports found that the KO mouse model of this gene complete loss of white fat (38), which is partially consistent with the phenotype we observed in our patient. Also, a cohort study found that mutations in genes related to the hedgehog pathway were significantly enriched in patients with pituitary stalk interruption syndrome, a common pituitary change of GHD (39). Combined with our results, this suggests that the hedgehog pathway may be critical to the development of the pituitary gland. However, the specific mechanism of GHD caused by this gene requires further verification through in vivo and in vitro functional tests.
Biallelic mutations in the LHX3 have been associated with combined pituitary hormone deficiency (MIM: #221750), characterized by impaired production of GH and other anterior pituitary hormones (29). It has been previously reported that there is variable phenotypic expressivity associated with different mutations in this gene (40,41). Recently, Jullien et al. reported dominant pattern of LHX3 might lead to a mild phenotype of combined pituitary hormone deficiency (41). Our findings support this new pathogenic pattern and the effect of heterozygous LHX3 variants on the synthesis and secretion of GH. The CREB3L4 gene is mainly expressed in prostate tissue. However, it is rarely expressed in brain tissue (42). Previous animal model studies found that male mice with Creb3l4 knockout had spermatogenesis disorder phenotype (43). In our study, the two GHD patients with CREB3L4 mutation were both females, and it is the first report that rare mutations of the gene are associated with human diseases. Therefore, we speculate that this gene may have sex-dependent expression characteristics, which may lead to the occurrence and development of GHD.

CONCLUSION
Our study was conducted based on a relatively large cohort of GHD patients with clinical characteristics and WES data. We found the discrepancy between serologic laboratory testing and molecular diagnostic methods. In addition to known causal genes, rare variants in 4 genes related to GH synthesis and secretion were identified by mutational burden analysis. These results suggest that predicted deleterious variants are enriched in GHD patients. Therefore, we recommend providing molecular testing for all possible GHD patients.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.