A GHRHR founder mutation causes isolated growth hormone deficiency type IV in a consanguineous Pakistani family

Background Isolated growth hormone deficiency (IGHD) is caused by a severe shortage or absence of growth hormone (GH), which results in aberrant growth and development. Patients with IGHD type IV (IGHD4) have a short stature, reduced serum GH levels, and delayed bone age. Objectives To identify the causative mutation of IGHD in a consanguineous family comprising four affected patients with IGHD4 (MIM#618157) and explore its functional impact in silico. Methods Clinical and radiological studies were performed to determine the phenotypic spectrum and hormonal profile of the disease, while whole-exome sequencing (WES) and Sanger sequencing were performed to identify the disease-causing mutation. In-silico studies involved protein structural modeling and docking, and molecular dynamic simulation analyses using computational tools. Finally, data from the Qatar Genome Program (QGP) were screened for the presence of the founder variant in the Qatari population. Results All affected individuals presented with a short stature without gross skeletal anomalies and significantly reduced serum GH levels. Genetic mapping revealed a homozygous nonsense mutation [NM_000823:c.G214T:p.(Glu72*)] in the third exon of the growth-hormone-releasing hormone receptor gene GHRHR (MIM#139191) that was segregated in all patients. The substituted amber codon (UAG) seems to truncate the protein by deleting the C-terminus GPCR domain, thus markedly disturbing the GHRHR receptor and its interaction with the growth hormone-releasing hormone. Conclusion These data support that a p.Glu72* founder mutation in GHRHR perturbs growth hormone signaling and causes IGHD type IV. In-silico and biochemical analyses support the pathogenic effect of this nonsense mutation, while our comprehensive phenotype and hormonal profiling has established the genotype–phenotype correlation. Based on the current study, early detection of GHRHR may help in better therapeutic intervention.

Both IGHD types IA and IB are autosomal recessive conditions characterized by short stature. In type IA, there is an absence of serum GH, and those affected produce anti-GH antibodies after GH treatment (5,6). In type IB, there are low (but detectable) serum GH levels and no evidence of antibody production after GH treatment (7). IGHD type II, however, is an autosomal dominant condition; as with type IB, low serum GH levels are detectable and no anti-GH antibodies are produced upon GH treatment (8). IGHD type III usually segregates in an X-linked manner and is often associated with agammaglobulinemia (9). IGHD type IV is a recessive condition characterized by early and severe growth failure. Those affected exhibit a reduced GH response to various provocation tests (e.g., tests for determining growth hormone level) and low insulin-like growth factor-I (IGF1) and IGF-binding protein-3 (IGFBP3) levels, but a good response to GH treatment. At the cellular level, we know that human GH binds to human growth receptor (GHR) molecules and induces signal transduction through receptor dimerization (10). When GHRH interacts with its corresponding transmembrane domains on somatotropic (GH-producing) cells, a G proteinmediated interaction with ion channels causes an increase in intracellular cAMP accumulation, which ultimately promotes GH release from secretory granules (10)(11)(12). Indeed, elevated cAMP causes protein kinase A to phosphorylate and activate CREB (13, 14), whose target genes include the pituitary-specific transcription factor Pit-1 (also known as GHF-1) (15)(16)(17). Pit-1 is a prototypic POU domain protein that is required for the proper regulation of GH gene activity in somatotropic cells, thereby providing a pathway by which a GHRH signal can lead to increased pituitary GH synthesis. Somatostatin, an inhibitory peptide, is thought to interact with this same signaling pathway via G protein-mediated suppression of the cAMP pathway (18,19). Indeed, the malfunctioning or underexpression of endogenous CREB protein in pituitary somatotropic cells causes somatotroph hypoplasia and dwarfism in mice (20). It is now clear that any disruption to this multistep GH signaling pathway can result in GH deficiency and ultimately lead to short stature and various other clinical problems.
Here, we analyzed a Pakistani family comprising four affected individuals with IGHD4. Whole-exome sequencing in this family revealed a nonsense mutation NM_000823:c.G214T:p.Glu72* in the third exon of GHRHR. The identified mutation presumably creates a premature terminator codon in the extracellular domain of the GHRHR and results in the synthesis of a truncated and nonfunctional receptor. In-silico findings and biochemical analysis support the pathogenic effect of the reported nonsense mutation.

Study design, declarations, and approvals
The base call files were converted to FASTQ files using bcl2fastq v2.20.0.422 (https://support.illumina.com/sequencing/sequencing_ software/bcl2fastq-conversion-software.html). The sequence reads in the FASTQ files were aligned to the human reference genome (GRCh37/hg19) using BWA-mem 0.7.17 (arXiv:1303.3997v2 [qbio.GN]) to generate BAM files. The BAM files were further processed for variant calling to generate VCF files using GATK v.3.8 (21). Finally, variant annotation and interpretation was performed using the EVIDENCE tool developed by "3billion" (South Korea) (22). To evaluate the authenticity of these predictions, variant annotation was carried out in parallel using wANNOVAR (https://wannovar.wglab.org/index.php) (23) and VariantStudio software (version 3.0) with the reference assembly "hg19." Variant prioritization was phenotype driven (i.e., based on the disease diagnosis), and previously reported genes were used as the training set (retrieved from the OMIM database, accessed on 20 August 2022). Primarily, the genes reported for IGHD were screened for the presence of possible candidate variants before the variant analysis was extended exome-wide to identify protein-disrupting mutations. For candidate variant prioritization, the analysis focused on filtering homozygous, missense and indel variants, and low-frequency and protein-encoding alleles. Finally, splice variants and untranslated region (UTR) variants were also investigated to exclude their possible involvement.

Preliminary computational screening
Before segregation analysis and 3D protein modeling, primary pathogenicity validation parameters were assessed. Initially, the biological significance of candidate variants was predicted using VarSome (24), which integrates multiple online pathogenicity prediction tools in a single platform. The gnomAD database (https:// gnomad.broadinstitute.org/) was then consulted to assess the population allele frequency.

Segregation analysis
To validate the segregation of the candidate variant(s) with the disease phenotype, Sanger DNA sequencing was performed on DNA samples provided by all available family members: III-1, IV-1, V-4, V-6, and V-10. The Primer3web tool (version 4.1.0) (26) was used to design suitable primers, and the mutation analysis was performed using BioEdit software (version 7.1) (27, 28).

Biomedical analysis
To assess receptor functionality and GH deficiency, a clinical test for serum GH level was performed for two affected participants (V-10 and V-11). Blood samples (5 ml each) were collected from both the affected individuals under fasting conditions and referred to the CAPcertified (College of American Pathologists) laboratory at Aga Khan University Hospital, Karachi, Pakistan. X-ray images of the hands and feet were taken for all affected individuals to identify any limb abnormalities, e.g., bone degeneration (rheumatoid arthritis) and abnormal bone fusion (syndactyly).

In-silico functional assays
Computational 3D structural modeling, protein-protein interaction analyses, and molecular dynamic simulation analyses were performed to predict the pathogenic effects of the identified variant(s) on the structure and function of the protein.

Protein structure prediction
Structural modeling of normal and mutant GHRHR was performed using the I-TASSER tool (29) and SWISS-MODEL (30). After structure prediction, 3D models of the wild-type and mutant GHRHR were superimposed to detect gross structural anomalies.

Molecular docking
ClusPro (31) was used to investigate the GHRH-GHRHR interaction. Namely, protein-protein docking of wild-type and mutant GHRHR was performed with its close functional interactor, GHRH, which was identified through the STRING matching database (32). Molecular visualization was performed using Chimera 1.13.1 (33) and LigPlot+ (version 2.1) (34).

Molecular dynamics simulation
AMBER package (version 16) (35) was used for molecular dynamic (MD) simulations of GHRHR models (wild type and mutant) and their docked complexes (35). All models were solvated in a rectangular box (8.0 Å) in the presence of counter ions to neutralize the system. Prior to MD simulations, each system was separately minimized, heated, and equilibrated. While keeping the protein model fixed, the studied systems were minimized. After minimization, the systems were heated and equilibrated at 300 K for 100 ps with a constant volume. Once the systems were equilibrated, all restraints were removed under constant pressure before running further equilibration steps at 300 K. MD simulation trajectory files were saved at 0.5 ps intervals. The PTRAJ/CPPTRAJ module (36) was used to analyze the saved MD simulation trajectory files, and Xmgrace (37) was used for visualization. The root mean square deviation (RMSD) and root mean square fluctuation (RMSF) were analyzed with the CPPTRAJ module. The RMSD values of both the wild-type and mutant structures were calculated using the starting structure as a reference frame and the deviation of the coordinates of a given set of atoms in a time interval. The RMSF was calculated to measure the fluctuations of each residue from their mean positions.

Qatar Genome Program data analysis
The Qatar Genome Program (QGP) is a population-focused study that aims to generate a whole genome sequence for all participants in the Qatar Biobank (QBB) (38,39). Qatar is a multicultural society where South Asian expats are considered to be one of the founder populations (40). The present study was conducted using information from~15,000 QGP genomes that have been comprehensively phenotyped based on biomedical information deposited in the QBB. The whole genomes and phenotype data were investigated to determine the frequency of rs121918117 in the local Qatari population as well as their associated clinical phenotypes (including hormonal profiles and liver and kidney functions).

Clinical findings of patients with IGHD type IV
A five-generation consanguineous Pakistani family affected by IGHD type IV was recruited through a local survey. The family consisted of five affected individuals, namely, four surviving patients [one female (V-4) and three male patients (V-10, V-11, V-12) aged 37, 26, 11, and 12 years ( Figure 1A), respectively] and one who died of an unrelated cause. The parents were first cousins with no previous history of IGHD and were asymptomatic with normal height, stature, and physique. The disease had emerged only in the most recent generation, and none of their predecessors exhibited this phenotype. All the patients had been born at full term with a normal delivery and no pregnancy complications. Genealogical analysis indicated a pattern of autosomal recessive disease inheritance ( Figure 1B).
We first collected anthropometric, hormone, and general clinical measurements from affected individuals with familial IGHD-associated dwarfism (Table 1). We saw that all the surviving patients were affected by severe growth retardation and short stature (i.e., dwarfism) ( Figure 1A). This finding supported the low or negligible serum GH levels, indicating abnormal GHRHR function. Specifically, all patients measured −6.9 to −8.0 SD below the mean height for the relevant age group in the healthy population (https://simulconsult.com/resources). By comparison, all unaffected individuals from general population had an average height of~5.5 feet. Digit curvature, namely, bilateral 4/5 clinodactyly, was observed in the hands of V-12; no other digit anomalies such as polydactyly or brachydactyly were observed. Moreover, no other gross anomalies were detected by X-ray analysis of the appendicular or axial skeleton.
The body mass index (BMI) of all patients was in the normal range, and no one exhibited abdominal fat deposits. The head-to-body ratio was normal and non-indicative of micro-or macrocephaly; the thyroid gland examination was also normal. All patients exhibited a delay in the eruption of permanent dentition, yet the skull shape/size and jaw structure were normal in all patients, with no symptoms of hypoplastic maxilla. The radiological findings of the affected individuals revealed normally developed carpal/metacarpal and tarsal/metatarsal bones. There was no polydactyly, syndactyly, or brachydactyly of the hands or feet, thus suggesting normal digit development ( Figure 1C).
The general clinical assessment of the patients did, however, indicate an underdevelopment of the body structure. All patients had a high-pitched voice. In medical examinations, no frontal bossing, depressed nasal bridge, or depressive behavior was observed. While fertility was not assessed in these patients, it was noted that all had reached puberty by 18 years of age (on average).
The findings of ophthalmologic, otorhinolaryngologic, neurologic, skeletal, and visceral organ examinations were normal. The patients had no current or previous history of diabetes, hypoglycemia, hypertension, cardiac problem, allergies (food or drug), dermal lesions, or autoimmune disorders.

Molecular analysis found founder mutation in GHRHR
Having comprehensively assessed the clinical phenotypes of all four affected patients, we next performed a WES analysis to determine the genetic basis of the disease. Here, we identified a recurrent homozygous nonsense mutation [NM_000823:c.G214T:p.(Glu72*); (rs121918117)] in the third exon of the GHRHR gene. Sanger sequencing confirmed the co-segregation of the identified mutation with the disease phenotype in all patients ( Figure 1D). This variant was not found in the unaffected family members or in our in-house database comprising an ethnically matched control population. The minor allele frequency of this allele in gnomAD was 0.0001755, with no known homozygote. According to a database survey, the variant was reported in South Asian regions only and not found in any European, American, Ashkenazi Jews, or East Asian communities, thus indicating that the variant is Asian-specific (allele count 34). Indeed, the (p.Glu72*) nonsense mutation has previously been mapped in three South Asian families (Indian, Sri Lankan, and Pakistani), which suggests a founder effect of this mutation in IGHD.
The ClinVar, VarSome, and InterVar databases reported rs121918117 as a pathogenic variant. Indeed, the premature stop codon presumably truncates the protein by deleting the C-terminus GPCR transmembrane (TM) domains and, in part, the GPCR-2 extracellular domains. According to TM domain integrity prediction, the truncated protein is unable to anchor within the membrane due to loss of all the downstream TM domains, causing loss of receptor function (Figure 2). We thus speculate that the premature stop codon and early truncation of the protein product caused by the identified nonsense mutation (p.Glu72*) might activate the nonsense-mediated mRNA decay pathway, resulting in low GH production in affected patients.

Protein structural analysis
Molecular modeling of wild-type and mutant GHRHR proteins and hormone-receptor docking We next aimed to determine the structure and interaction studies. We first found that superimposition of the 3D models of both wildtype and mutant GHRHR was not possible (Figures 3A-C), confirming that the c.G214T: p.Glu72* nonsense mutation causes structural distortion of GHRHR due to protein truncation.
We next monitored the interactions between the wild-type and mutant GHRHR receptor proteins with GHRH. Compared with the full-length wild-type protein, the truncated GHRHR protein showed a different interaction pattern with respect to residue numbers, residue identity, and bonding category ( Figures 3D, E). The wild-type receptor could interact with GHRH via 14 amino acid residues (Asp395, Lys130, Asp183, Ser124, Arg384, Asp190, Asp191, Lys329, Arg405, Lys417, Glu123, Gln381, Thr128, and Tyr125) and 27 bonding forces involving 23 H-bonds and four salt bridges. By contrast, the mutant receptor docked with GHRH via only six amino acids (Leu63, Cys64, Asp37, Arg4, Asp2, and Glu72) and eight hydrogen bonds. These findings suggest that the mutation protein binds more strongly with the interactor compared to the wild-type GHRHR.

Findings on molecular dynamic simulation analysis
Next, we evaluated the dynamic behavior of the wild-type and mutant GHRHR proteins and their interactions with GHRH by 100ns MD simulation. We generated models of the wild-type-GHRHR-GHRH and mutant-GHRHR-GHRH complexes ( Figures 3D, E,  respectively). We determined the atomic positions of the modeled GHRHR protein structure and compared them with those of the mutant-GHRHR protein to examine the convergence and durability of the MD simulations. The trajectories of the MD simulations were analyzed in a multistep process to decode the backbone stability and residual flexibility. The RMSD was first calculated as a measure of the average distance between the backbone carbon alpha atoms in the overlay frames. The average distance between the backbone carbon alpha atoms of the overlay frames was then calculated, and the stability of the simulated C atoms of the protein was determined by plotting the modeled wild-type-GHRHR, mutant-GHRHR protein, wild-type-GHRHR-GHRH protein, and mutant-GHRHR-GHRH proteins as a function of time. We found evidence for a slightly greater deviation of the wild-type-GHRHR (RMSD mean 5.93364 Å, max. 8.2271 Å at 85,361 frames) than the mutant-GHRHR (RMSD mean 5.02231 Å, max. 7.2316 Å at 8,858 frames) (Figure 4, upper panel). In both the wild-type-and mutant-GHRHR, a constant displacement was observed in the first 25 ns of the MD simulation, indicating that the atoms were out of alignment and the structure was unstable. While the mutant-GHRHR showed minor fluctuations for the remaining simulation period, the wild-type-GHRHR-GHRH showed considerable deviations (RMSD mean 8.79908 Å, max. 8.79908 Å at 42,247 frames), indicating that the docked complex was stable after 60 ns (Figure 4, lower panel). By contrast, the RMSD plot for the mutant-GHRHR-GHRH revealed an unstable system that stabilized during the MD simulation and remained stable until the end. The plot showed considerable deviations (RMSD mean 6.73827 A, max. 8.7608 Å at 97,250 frames), indicating that the atoms were out of synchrony and the system was unstable. We also performed an RMSF analysis to assess the fluctuation of residues in the wild-type and mutant docked complexes. Here, the mutant structure showed much greater stability than the wild-type structure, which instead underwent marked fluctuations. Membrane topology of GHRHR: images of the transmembrane topology of the GHRHR protein.

Estimation of different binding free energies using MMGB/PBSA
We assessed the binding free energies of the wild-type and mutant GHRHR proteins complexed with GHRH. We calculated these based on the simulated trajectories using the MMGBSA and MMPBSA approaches. The net total energy for the wild-type-GHRHR-GHRH structure was −215.747, while that for the mutant-GHRHR-GHRH structure was −99.2622. In comparison to the MMGBSA approach, the MMPBSA net binding energy showed that the interaction energies of both systems were highly favorable and stable ( Table 2).

QGP data outcomes
To determine the frequency of rs121918117 (GHRHR) in the population, we leveraged the QGP database and identified two heterozygous individuals (QGP1 and QGP2) with an alternate allele frequency of 0.00006817. Both QGP1 and QGP2 were clinically normal and asymptomatic although their seated height was slightly reduced compared to that of the average value of the QGP-cataloged individuals. In terms of standing height, QGP1 was in the normal range, while QGP2 had a short stature (−0.5 SD). All other available health indicators were in the good and excellent categories, with the exception that QGP2 was classified as obese and QGP1 had systolic/ diastolic BP and pulse rates tending to 89/62 mmHg and 68 beats per minute.

Discussion
IGHD is an inherited condition of inadequate secretion of GH by the pituitary gland (41) that results in severe short stature (9). To date, 82 mutations in GHRHR have been identified in IGHD patients of multiple ethnicities. Of these mutations, six are nonsense, 12 are splice disrupting, 44 are missense, 16 are small indels, and four are   MD simulation properties: molecular dynamic simulation graphs of wild-type, mutant, and their docked complexes with GHRH.     (Table 3). Here, we have added to the portfolio of GHRHR mutations underlying IGHD after performing a WES analysis of a large, consanguineous Pakistani family comprising four individuals suffering from IGHD4. Our WES and Sanger sequencing analysis revealed a nonsense GHRHR mutation [NM_000823:exon3:c.G214T:p.Glu72*] that segregated with the disease phenotype. This nonsense mutation results in the generation of a truncated protein devoid of the Cterminus membrane spanning TM domains. We speculate that this truncation triggers the nonsense-mediated decay of defective mRNA, which is the most likely reason for the loss of functional integrity of the GHRHR and hormonal signaling. Indeed, this phenomenon can be correlated with reduced or absent GH levels in patient serum and is consistent with the extremely reduced levels of serum GH in affected individuals observed in our hormone analysis. In expression studies of truncated and chimeric epitope-linked GHRH receptors (i.e., GHRHR), DeAlmeida and Mayo (79) first identified the regions vital for the interaction of the receptor with GHRH. They investigated the fate of two truncated receptors; GHRHdelta-N, without the N-terminal domain (the region between the signal sequence and the first TM domain) and GHRH-delta-C, with trimmed downstream TM domains. In both experiments, they found diminished receptor-ligand interactions, showing that neither the Nterminus extracellular domain nor the membrane-spanning Cterminus domains were independently capable of interacting with GHRH. Instead, it was determined that although the N-terminal extracellular domain is necessary for ligand binding, the transmembrane domains and associated extracellular loop regions of the GHRHR are vital for the specific interaction with GHRH (79). The current findings also support this notion, wherein the C-terminal deletion mutation p.Glu72* also exhibits diminished interaction of GHRHR-GHRH as presented by the reduced level of growth hormones.
The p.Glu72* mutation detected in this study was previously identified by Wajnrajch et al. in an Indian family with profound GH deficiency (55), with the patient's phenotypes found to be similar to the "little mouse" harboring Ghrhr gene mutation(s). Furthermore, the human counterparts of the mutations in the "little mouse" model (especially Asp60Gly mutation in mouse being closest to human Glu72*), were associated with comparable phenotypic features in terms of stature and hormone level. The phenotypes of patients enrolled in our study showed considerable similarities with those reported by Wajnrajch et al.: all patients exhibited phenotypes as seen in the mouse model although they did not exhibit obesity or frontal bossing. Genetic studies, including our own, have shown that carriers of GHRHR mutations are asymptomatic and of normal height. In a clinical investigation of a cohort of 76 adults (aged 25-75 years) carrying a heterozygous GHRHR mutation, Pereira et al. (80) found no association between the GHRHR mutation and short stature. We also found that all the family members (including both wild-type and heterozygous carriers) were under 5 feet in height. However, the slightly reduced sitting height of individuals from the QGP cohort remains to be investigated.
The (p.Glu72*) nonsense mutation we identified here has also been mapped in three different South Asian families from Bombay (India), Delft (an island located between Sri Lanka and India), and Sindh (Pakistan). The relative geographic connection of these families indicates the possibility of a common origin. Wajnrajch et al. performed a microsatellite-based haplotype analysis to evaluate this hypothesis among these three apparently unconnected families from the Indian subcontinent harboring the common GHRHR mutation (p.Glu72*) (81). Their evolutionary analysis revealed that these families originated from a common ancestor between 1,350 and 2,700 years ago when the three traditionally and linguistically distinct families separated from their common ancestry (47). Our identification of the same p.Glu72* mutation in a Pakistani family increases the probability of its founder effect. Furthermore, since the family in our study originated from the southernmost part of Pakistan near the Sindh Province, our study also supports the possibility of Indo-Dravidian (or Indo-Aryan) migration (81).

Conclusions
In the current study, we identified a nonsense mutation (p.Glu72*) of GHRHR in a consanguineous Pakistani family affected by IGHD4. Biochemical analysis revealed extremely low GH levels, suggesting functional loss of the GHRHR. We confirmed the pathogenicity of this variant by in-silico analysis. To the best of our knowledge, this is the fourth report of a GHRHR mutation NM_000823:c.G214T:p.(Glu72*), and together with previous reports, our data provide further support for the founder effect of p.Glu72* in South Asian families. Our data will help to establish the genotype-phenotype correlation in IGHD type IV and, thus, will support molecular diagnosis in clinical practice. Evaluating GHRHR gene integrity is an important aspect of genetic evaluation in those of short stature, as early genetic diagnosis might support precision medicine treatment options. Indeed, successful treatment, i.e., hormonal therapy based on GHRHR detection, may help maintain a normal growth pattern and would subsequently enable affected patients to live a confident life in society and avoid any secondary implications of their disease, such as delayed puberty, fertility, and physique.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/supplementary material.

Ethics statement
This medico-genetic study was approved by the Institutional Ethical Review Board of Gomal University D.I. Khan (IRB# 04/ ERB/GU), KPK, Pakistan. Informed consent was obtained from all study participants and their parents/guardians. The patients/ participants provided their written informed consent to participate in this study. Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.

Author contributions
SA, MA, IAhma, and ShA assisted with the patient recruitment, conducted the experiments, and performed the data analysis. SWA performed the bioinformatic analysis. IAhme, SN, and MZ performed the sequence data analysis. MK, AAA, and KF conceptualized and supervised the study, obtained funding to support the study, and drafted the manuscript. All authors contributed to and approved the final version of this manuscript.

Funding
This study was supported by Sidra Medicine, through funding number 20.6545-632157 (provided to AAA, Sidra Medicine, Doha, Qatar).