This study was approved by the Ethical Review Board of Gomal University, Dera Ismail Khan, Pakistan. Informed consent to perform genetic, molecular, and clinical analyses and to publish patient data and images was obtained from all study participants. The family was identified in Tehsil Paroa of District Dera Ismail Khan, Pakistan, through a local street-to-street survey. The genealogy was ascertained to assess the mode of disease inheritance and determine the level of consanguinity between the parents. Then, the blood samples were collected, and DNA was isolated using a GeneJET Genomic DNA purification kit (Thermo Fisher Scientific, USA, Cat# K0721), according to the manufacturer’s instructions.

All the patients (V-4, V-10, V-11, V-12) were siblings and exhibited the same phenotype ( Table 1 ). Therefore, due to the high probability of harboring common genetic variants, WES was performed on two randomly selected patients (V-4 and V-10). A sequencing library was constructed using an xGen Exome Research Panel v2.0 Kit (Integrated DNA Technologies, Coralville, IA, USA) and sequenced on a NovaSeq 6000 (Illumina, San Diego, CA, USA). 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 [q-bio.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.

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. To validate the clinical significance of identified candidate variants [as per ACMG guidelines (https://www.acmg.net/)], InterVar (25), VarSome (https://varsome.com/), and ClinVar databases (https://www.ncbi.nlm.nih.gov/clinvar/) were explored. Being a receptor protein, the integrity of transmembrane domains was assessed and visualized using TMHMM (https://services.healthtech.dtu.dk/service.php?TMHMM-2.0), Protter (https://wlab.ethz.ch/protter/start/), and Phobius (https://phobius.sbc.su.se/) prediction tools.

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).

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 CAP-certified (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).

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.

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).

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.

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.

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 ).

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.