Our index case, individual II.1, was extremely tall in adolescence and treated with high dose testosterone esters. At 56 years of age, we measured his height at 208.6 cm, which would translate to 211 cm (4.0 SDS) at 21 years (due to shrinking with age) (9). Yet, this may still be an underestimation in view of the treatment with testosterone esters in adolescence. Sitting height was 104 cm (sitting height/height -1.0 SDS) (10) and head circumference 60 cm (1.3 SDS) (11). He was married to a tall woman (II.2; 3.2 SDS) and had a tall son (III.1; 2.8 SDS). His sister (II.4; 3.2 SDS) married a man with normal height (II.3; 1.7 SDS) and had a daughter of normal height (III.2; -0.7 SDS). Height data were adjusted for secular trend and shrinking (9). Subjects II.1, II.4 and III.1 were examined in childhood and adolescence by experienced paediatric endocrinologists, who documented a stable height SDS and an absence of any signs of growth overproduction or dysmorphic syndrome. Further information on subjects and growth measurements, and on the calculation of adjusted height and shrinking, is given in the Supplementary Information .

Candidate genes identified in this study were also checked in all tall family members that we previously described in Weiss et al, 2021 (12).

Exome sequencing was conducted on genomic DNA from peripheral blood leucocytes. In brief, coding genomic regions were enriched using the SureSelect XT Human All Exon Kit V7 (Agilent Technologies, Santa Clara, CA, USA) for subsequent sequencing as 2x125 or 2x100 bp paired-end reads on a HiSeq2500 or NovaSeq6000 system (Illumina, San Diego, CA, USA). Generated sequences were analysed using the megSAP pipeline (https://github.com/imgag/megSAP). Clinical variant prioritization included different filtering steps and was conducted either A) in a single-mode (MAF [minor allele frequency] in 1000g, ExAC, gnomAD or in an in-house database set to <0.1% [autosomal-dominant] or <1% [autosomal-recessive]) or B) in a multi-mode (prioritizing only those variants shared by family members with tall stature), with three different sub-types: multi-sample dominant, multi-sample compound-heterozygous [both with a MAF cut-off of <0.1%] and multi-sample recessive [MAF cut-off: <1%]).

Variants of interest were confirmed by Sanger sequencing. Polymerase chain reaction (PCR) with primers indicated in Supplementary Table 1 was performed with HotStar Taq Polymerase and products analysed on agarose gels and subsequently sequenced using an ABI machine (Azenta Life Sciences).

We also checked for duplications including SHOX and IGF1R duplications and XYY in all the exome data sets and can exclude complete duplications. Partial duplications of these genes are also unlikely, but we cannot totally exclude them due to methodological limitations. We also checked for IGF-1R variants in the exome data set and re-checked the three known NPR2 mutations from the literature that are associated with tall stature and found no such variants in this family.

Variants were also checked in the “normal Non-Finnish control cohort” from gnomAD, which does not specify the height of the control individuals. Therefore, it is possible that single gene variants can be present in controls at low allele frequencies.

QIAGEN Ingenuity Pathway Analysis (IPA; https://digitalinsights.qiagen.com/) was applied to predict functional connections. IPA integrates selected omics data sets with mining techniques. Interpretation was in the context of protein networks that comprise protein-protein interactions and related biological functions and canonical signalling pathways.

TGP (https://browser.1000genomes.org); GnomAD (https://gnomad.broadinstitute.org/); CADD score (https://cadd.gs.washington.edu); dbSNP (https://www.ncbi.nlm.nih.gov); GTEx database (www.gtexportal.org); PROVEAN/SIFT (http://provean.jcvi.org/index.php); Polyphen2 (http://genetics.bwh.harvard.edu/pph2/); GWAS (www.ebi.ac.uk); Disease knowledge portal (https://cvd.hugeamp.org); IPA Ingenuity Systems (https://digitalinsights.qiagen.com/).

Gene expression analysis in mouse growth plate was carried out to analyse (1) growth plate specificity: comparing mRNA expression in growth plate of 1-week-old mice with expression in three different soft tissues (lung, kidney and heart) using microarray (13) (2); spatial regulation: comparing mRNA expression in resting, proliferative and hypertrophic zones in the 1-week-old mouse growth plate by RNA-Seq; and (3) temporal regulation: comparing mRNA expression at 1 and 4 weeks of age in proliferative and hypertrophic zones of mouse growth plate by RNA-Seq (14). For spatial comparison, chondrocytes from different zones of 1-week-old mouse proximal tibia growth plate were isolated by laser capture microdissection as previously described (15). The spatial analysis of gene expression in the growth plate has previously been validated by demonstrating that it correctly identifies the spatial expression pattern of known growth plate zone markers (14). Sequencing libraries were prepared using a TruSeq Stranded mRNA Prep Kit (Illumina). Sequencing was performed via a paired end 75 cycle on Illumina HiSeq 2500. RNA-Seq reads were trimmed with cutadapt (-AAGATCGGAAGAGCACACGTCTGAACTCCAGTCAAAGATCGGAAGAGCGTCGTGTAGGGAAAGA GTGT -overlap 6 -q 20 -minimum-length 25) and aligned using STAR (2 pass alignment) to mouse mm10 reference genome sequences. Log2 fold changes and p-values (adjusted for multiple comparison using the Benjamini Hochberg correction) were generated by DESeq2.

Exome sequencing was carried out in the index case (individual II.1) and three additional tall members of a three-generation family (individuals I.1, II.4, III.1; Figures 1A, B ). All individuals were healthy and without organomegaly nor dysmorphic features. Filtered variants predicted as damaging (based on Provean, SIFT and PPH2 scores) and not rated as polymorphisms (by gnomAD) and with high CADD scores, as well as shared among the related tall individuals, were first validated by Sanger sequencing and then checked for absence in a family member with normal stature (individual III.2, Figure 1A ). Missense variants in the following six genes were classified as likely pathogenic: CNGB1, GGTLC2, HSD3B2, IFT140, NAV2 and SCAF11 ( Table 1 ). Comparison to the normal Non-Finnish controls from gnomAD revealed allele frequencies of the individual variants between 10^-3 to 10^-5. The likelihood that a control individual carries several of these variants at the same time is minimal.

CNGB1 (cyclic nucleotide gated channel subunit beta 1) encodes an ion channel protein; GGTLC2 (gamma-glutamyltransferase light chain 2) is involved in RNA binding and spliceosomal complex assembly; HSD3B2 (3-beta-hydroxysteroid dehydrogenase) is essential for synthesis of multiple steroid hormones; IFT140 (Intraflagellar transport 140 homolog) is indispensable for intraflagellar transport in primary cilia; NAV2 (Neuron navigator 2) is an all-trans retinoic acid responsive gene; and SCAF11 (SR-related CTD associated factor 11) encodes an RNA binding protein and is a regulator of spliceosome assembly.

From all affected amino acids of the six genes, species conservation was highest for NAV2 (conserved between human and Drosophila) and lowest for HSD3B2 ( Table 2 ). GTEx expression analysis indicated that CNGB1 is expressed in a brain-specific way, GGTL2 in a testis-specific and HSDB2 in an adrenal gland-specific manner, whereas expression of IFT140, NAV2 and SCAF11 is considered ubiquitous ( Supplementary Table 2 ).

We next determined whether the identified six genes were associated with height in GWAS studies using the Disease knowledge portal (https://cvd.hugeamp.org). We found IFT140, NAV2 and SCAF11 to be associated with height with genome-wide significance (P-values between 10^-11 and 10^-14), while no associations were reported for CNGB1, GGTLC2 and HSD3B2 ( Supplementary Table 3 ). Because it is unknown whether common genetic variants (SNPs) associated with height are concentrated around the same genes where rare variants are identified, we decided that further evidence was needed.

To further address the biological relevance of the identified gene variants, network analysis by Ingenuity Pathway Analysis (IPA) was carried out. IPA was used on the six novel candidate genes for tall stature together with 86 genes previously associated with tall stature by GWAS and literature data ( Supplementary Table 4 ) (12, 20, 21). The newly identified six candidate genes were analysed to assess whether functional connections exist with known protein networks ( Figure 2 ; Supplementary Figure 1 ). Our network analysis indicated direct protein-protein interactions between GGTLC2 and TLX3 (22), between IFT140 and CRK (23), and between SCAF11 and YWHAE (24) ( Figure 2 ). Protein-protein interactions also exist between SCAF11 and the estrogen receptor ESR1 and the transcription factor TCF7 (25), and between HSD3B2 and ESR1 (26) ( Figure 2 ). The all-trans retinoic acid-responsive NAV2 protein interacts with YWHAE, KRAS, HRAS and the multiprotein complex Vitamin D receptor-interacting protein mediator/TRAP (for thyroid hormone receptor-associated protein), which is highly associated with height (3.79E-28) (27–29). TRAP interacts with various transcription factors and acts as a transcriptional co-activator. CNGB1 did not interact with any of the known tall stature- associated genes.

Together, these data provide experimentally-based evidence that five of the six candidate genes are involved in protein-protein interactions with other known genes previously associated with height.

To gain insight into the functional role of these candidate genes in bone growth, we assessed mRNA expression in murine growth plates. Growth plates are the cartilaginous structures residing near the ends of long bones, responsible for bone elongation and therefore linear growth. We evaluated microarray and RNA-Seq data (13, 14) in mouse growth plates (and other tissues) focusing on three comparisons. First, we compared 1-week old whole growth plate tissues with three soft tissues (lung, kidney, and heart) to compare expression of each gene of interest in the growth plate to expression in other tissues (13). Second, we compared expression in three different zones: the resting, proliferative, and hypertrophic zone in 1-week old mice by RNA-seq. We asked whether a gene showed a spatial gene expression pattern suggesting functional importance for stem cell maintenance (resting zone), cell division (proliferative zone), or terminal differentiation (hypertrophic zone). Third, we compared 1- versus 4-week old proliferative and hypertrophic zones using RNAseq (14), asking whether a gene showed a temporal decline or increase in gene expression that would suggest regulation of proliferation or differentiation, as growth slows with age.

Our analysis showed that Ift140 and Nav2 were well expressed in the growth plate, showing much higher levels compared to soft tissues analysed ( Table 3A ). Ift140 and Nav2, and to a lesser extent Scaf11, also showed significant downregulation of gene expression from proliferative to hypertrophic zone at 1 week. Ift140 with fold change of -0.86 and a p-value of 2.98E-06 as well as in particular, Nav2 with a fold change of -2.78 and a p-value of 1.91E-62 showed high changes between proliferative compared to hypertrophic zone, suggesting that they might be involved in chondrocyte proliferation and/or hypertrophic differentiation ( Table 3B ). Interestingly, Nav2 expression decreases in the proliferative zone from 1- to 4-week, further suggesting that Nav2 may act as an important regulator of chondrocyte proliferation ( Table 3C ).

In summary, our expression data are consistent with a functional role of both Ift140 and Nav2 in the growth plate, with Nav2 being the strongest of the two.

To functionally analyse NAV2 in an animal model, we used a CRISPR/Cas9-based knockout approach in the African clawed frog Xenopus laevis. Due to restrictions by animal protection, we could only raise tadpoles until stage 45, i.e. about 5 days after fertilization. Thus, we could test for a developmental function of the Xenopus nav2 gene during early body growth. To target the highly conserved exon harbouring the human arginine to cysteine missense mutation (R2142C), a single-guide RNA (sgRNA1) was designed. Injection of this sgRNA1 should generate a nonsense mutation or frameshift, leading to a truncated protein lacking the last six exons of the frog Nav2 protein, or alternatively, result in nonsense-mediated decay of the mRNA. We performed three injection experiments, injecting one-cell stage embryos and raised them until tadpole stage to quantify body sizes. Successful gene editing of injected specimens was evaluated by target sequence amplification and Sanger sequencing, demonstrating mutagenesis in about 35% of gene copies (for details, see Method section). Total tadpole body length of injected embryos showed increased body sizes in all three experiments (Exp1: +0.97%, Exp2: +3.36%, Exp3: +4.23%) ( Figure 3 ). As the expected size differences were small and wild type body sizes differ naturally from batch to batch due to minor age differences, we could only compare between mutant and the corresponding wildtype. To be able to combine the data of the three different clutches, we normalized body lengths of each experiment to its controls. Statistical analysis of the combination of relative body sizes showed a significant increase after nav2 gene editing.

We also designed a second sgRNA2 targeting a large 5’ exon, to cause a complete nav2 gene knockout (see Methods, for details). Sequencing confirmed a mutagenesis rate of 75-93%. Again, no early phenotype was observed in two different experiments, but when raised to tadpole stage 45, about 50% of injected specimens started to develop organ defects and edemas from stage 42 onwards ( Supplementary Figures 2A–C ). The massive body defects resulted in general body size reduction, potentially overruling any increase in length ( Supplementary Figure 2D ). Yet, when we compared the other 50%, crispants without gross defects, we also found mildly but not statistically significant increase in body sizes, as previously in the knockout approach with sgRNA1 (Exp4: +0.47%, Exp5: +1.49%) ( Supplementary Figure 2E ).

Taken together, our data on nav2 mutants provide independent functional verification of mildly increased body size of 1-4% in five different experiments already at early stages of tadpole development.

After we identified the candidate genes in this family, we re-examined exome sequencing data from a previously reported Dutch family with tall stature (12) and identified a rare NAV2 variant in two family members, but not in the other candidate genes of the current study. This NAV2 variant [c.G6112A, p.(E2038K)] was predicted as probably damaging by PPH2, damaging by Provean, deleterious by Sift and disease-causing by Mutation taster, with a CADD score of 26.8. The affected variant was conserved between human, mouse, chicken and Xenopus.

The variant [c.G6112A, p.(E2038K)] was identified in the two tallest family members, the index case (height 3.5 SDS) and his father (height 3.2. SDS). The variant was not carried by the sister of the index case (height 2.6 SDS) nor in the paternal uncle (height 2.9 SDS). All tall family members, however, shared damaging heterozygous variants in five different genes. This provides further evidence that damaging NAV2 variants, in conjunction with other genes, can be associated with tall stature phenotypes.