As of December 2022, among 2,764 patients undergoing treatment for congenital scoliosis at Peking Union Medical College Hospital, 40 Han Chinese patients were confirmed to have congenital kidney anomalies with associated urinary tract abnormalities. The phenotype evaluation for congenital scoliosis and congenital anomalies of the kidney and urinary tract was performed by different specialists in orthopedics and nephrology. Patients were categorized to different phenotypical subgroups. CS type I was characterized by simple vertebral dysplasia, CS type II by vertebral segmentation defects, and CS type III represented a combination of the previous two conditions. To note, severe infantile idiopathic scoliosis (IS) was also included. The diagnostic criteria have been reported in prior study (Li et al., 2022). Written informed consent was obtained from all participants at the time of recruitment. This study was approved by the Ethics Committee of Peking Union Medical College Hospital (I-22PJ976).

Next-generation sequencing (NGS), encompassing both exome sequencing (ES) and genome sequencing (GS), was carried out on DNA samples extracted from peripheral blood of the study participants. ES was conducted in paired-end mode on the Illumina HiSeq 2000 platform, maintaining an average depth of coverage at 70X. To create paired-end libraries, three different capture kits were utilized: the xGEN targeted capture kit v2 (IDT, Coralville, IA, United States), VCRome SeqCap EZ Choice HGSC 96 Reactions (Roche, Pleasanton, CA, United States), and SureSelect Human All Exon V6 + UTR r2 core design (Agilent, Santa Clara, CA, United States).

For GS, sequencing libraries were generated following the manufacturer’s recommended protocol using the KAPA Hyper Prep kit (KAPA Biosystems, Kusatsu, Japan). Subsequently, PCR-free sequencing was executed on either the Illumina HiSeq X-Ten or NovaSeq platform, achieving an average depth of coverage of 35X.

The study employed the Peking Union Medical College Hospital pipeline for variant calling and quality control (Chen et al., 2021; Zhao et al., 2021). We utilized the XHMM tool to analyze copy number variations in whole-exome sequencing (WES) data and employed the CNVpytor tool to analyze copy number variations in whole-genome sequencing (WGS) data. Raw sequencing reads were aligned to the GRCh37 genome, and variants were identified using Sentieon software (Supplementary Methods). Following rigorous quality control of both cariant and sample level, 12,041 variants from 40 cases and 2,764 controls were retained (Supplementary Methods). In brief, this process involved hard filtering, population-based filtering, and Variant Quality Score Recalibration (VQSR). The final dataset underwent principal component analysis (PCA) and identity-by-descent (IBD) analysis for further refinement.

The annotation of small variants (i.e., SNVs and indels) was carried out using the Variant Effect Predictor (VEP), which predicted the effects of the variants on transcript and protein sequences. The LofTee (https://github.com/konradjk/loftee) and dbNSFP plugins were used to generate bioinformatic predictions. Transcript-level information was annotated using the NCBI RefSeq database (https://www.ncbi.nlm.nih.gov/refseq/, accessed on 2022/05/25), with the RefSeq transcript labeled as “canonical” by ENSEMBL (https://www.ensembl.org/, accessed on 2022/05/25) being selected. CNVs were annotated using AnnotSV (v3.0), which compiled genomic and genetic information including genes and functional elements involved, gene-based annotations, genomic content around breakpoints, and overlapping records from population and disease databases.

The Genome Aggregation Database (gnomAD, https://gnomad.broadinstitute.org/) was utilized for annotating the population frequencies of variants. Rare variants with a gnomAD population-max allele frequency ≤ 0.1% and cohort allele frequency ≤ 0.1% were selected for gene-based genetic burden analysis. Subsequently, the retained variants were further annotated with transcript-level information from the NCBI RefSeq database (https://www.ncbi.nlm.nih.gov/refseq/). Each variant was categorized into different mask levels (from 1 to 6) and assigned weight values (ranging from 1 to 0) based on variant type and bioinformatic prediction results. The mutational burden of a given gene in each individual was designated as the maximum weight value among all ultra-rare variants carried by that individual. The ACAT package 50 was employed to conduct weighted mutational burden tests for each gene, with a minor allelic count (MAC) threshold set to 5 (weight ≠ 0) (See Supplementary Methods for details).

We obtained aggregated mutation burdens of gene sets related to biological processes, molecular pathways, or disease phenotypes from the following databases: BioCarta (http://www.biocarta.com/), Gene OntoManhattan plot depicting pathway-based rare variant logy (GO, http://geneontology.org/) for biological processes, Human Phenotype Ontology (HPO, https://hpo.jax.org/), Kyoto Encyclopedia of Genes and Genomes (KEGG, https://www.genome.jp/kegg/), Mouse Phenotype Ontology (MPO, https://www.mousephenotype.org/), Pathway Interaction Database (PID, http://pid.nci.nih.gov/), Reactome (https://reactome.org/), and WikiPathways (WP, https://www.wikipathways.org/). For a specific set of genes, the mutation value (equivalent to the gene-based test) for each gene was combined (for each person) into a total score. This score was then employed in a logistic regression model to estimate the likelihood of a case or control status. Enrichment in the analyzed gene sets is indicated by odds ratios and p-values.

For each detected CNV, the genomic position was matched against known gene boundaries to classify CNVs that are entirely within a gene (intragenic), span multiple genes, or lie between genes (intergenic). For each gene, CNVs from individual samples were aggregated to produce a gene-level CNV score. This was done to evaluate the collective impact of CNVs on individual genes and to facilitate gene-based association testing. Gene-based CNV scores were statistically assessed for their association with the coexistence of CS and CAKUT using the chi-square test.

In our study cohort, we examined a total of 40 patients with CS-CAKUT. The gender distribution in this cohort roughly follows a male-to-female ratio of 2:1. Notably, among the patients categorized under the scoliosis subtype, CS Type I exhibited the highest prevalence, constituting 40% of the entire patient population. These spinal curvatures primarily manifested within the thoracic region, specifically between vertebrae T1–T12. As for the CAKUT phenotype, unilateral renal agenesis displayed the most frequent characteristic in the cohort, affecting 32.5% of the patient cohort. These findings are summarized in Table 1.

Additionally, the characterization of patient traits highlighted that those with CS Type I in conjunction with renal abnormal (KA and EK-HK) displayed a syndromic deformity frequency of 33.3%, accounting for 4 out of the 12 patients in this subgroup (Figure 1). Consequently, it is common to observe variations in other systemic phenotypes in these individuals. Moreover, among the subset of 5 patients diagnosed with OU CAKUT, 3 patients experienced syndromic deformities, constituting 60% of this subgroup. Particularly noteworthy is that 2 of these patients exhibited respiratory system phenotypes, specifically pulmonary hypoplasia or restrictive ventilatory defects. It is crucial to emphasize that both of these patients were diagnosed with CS Type III in conjunction with OU. Due to the more severe symptoms associated with CS Type III, our findings indicate that it often presents with comorbidities such as Rib Fusion and Hypoplasia of the ribs. Another notable discovery was the occurrence of dextrocardia in two KA patients. Specifically, one had CS Type I combined with a horseshoe kidney, and the other had CS Type III combined with unilateral renal agenesis. Both individuals also present with severe scoliosis and comorbidities beyond skeletal and urogenital issues (Figure 1).

In summary, our study has revealed the intricate and diverse phenotypic variations among patients with CS-CAKUT, encompassing aspects of the spine, urinary system, and other bodily systems.

Building on these clinical observations, we employed a weighted gene-based burden test to compare the genome-wide frequencies of rare coding variants between CS—CAKUT cases and controls.

In comparison to the control group, SPRTN exhibited the strongest correlation (p = 1.41E-26) (Table 2; Figure 2A). The protein encoded by SPRTN potentially participates in DNA repair during replication in the presence of damaged DNA. Mice lacking a similar protein display chromosomal instability and a premature aging phenotype. Mutations in this gene are associated with Ruijs-Aalfs syndrome (RJALS), characterized by intellectual developmental delay, atypical facial features, short stature, other skeletal abnormalities, and neurological issues. Previous studies have also linked SPRTN to premature aging and hepatocellular carcinoma (Hiom, 2014; Lessel et al., 2014). Particularly intriguing was our discovery concerning protein tyrosine phosphatase non-receptor 11 (PTPN11) (p = 1.95E-21), an intracellular protein tyrosine phosphatase (PTP) (Kamiya et al., 2014) (Table 2; Figure 2A). Germline point mutations in the PTPN11 gene lead to Noonan syndrome (NS; OMIM 163950) and LEOPARD syndrome (LS; OMIM 151100), both sharing several common features, notably skeletal abnormalities, consistent with our study cases. Therefore, PTPN11 emerges as a compelling candidate gene for CS-CAKUT.

All identified deleterious variants, p. Asn308Asp, p. Arg351Leu, and p. Thr468Met, in CS-CAKUT patients are located within the Phosphatase domain (Figure 2B). This domain plays a vital role in promoting extracellular matrix mineralization and is linked to hypophosphatasia (Yadav et al., 2011). The three patients carrying loss-of-function variants exhibit severe scoliosis and KA, with additional phenotypes mainly affecting the cardiovascular and pulmonary domains. Notably, one of these patients has dextrocardia. Previous reports have closely associated PTPN11 p. Asn308Asp with Noonan syndrome, characterized by typical facial features, cardiac issues, growth problems (Beneteau et al., 2009), and varying cognitive impairments, aligning with our observed phenotypes. Furthermore, PTPN11 p. Thr468Met has been found to cause not only Noonan syndrome but also exhibit a close connection with LEOPARD syndrome (Conboy et al., 2016; Zhang et al., 2016). This highlights the significant role of PTPN11 in skeletal growth and remodeling.

In our effort to thoroughly investigate functional pathways enriched with genetic variations, we conducted a comprehensive analysis using the pathway-based gene burden approach. Figure 3 meticulously delineates the top 15 pathways that emerged as significant contributors in our study.

Among these pathways, three stood out due to their significant enrichment in cell signaling mechanisms: “Netrin Mediated Repulsion Signals (OR = 25.3, p = 0.0003),” “SPRY Regulation of FGF Signaling (OR = 16.3, p = 0.0004),” and “Raf Independent MAPK1_3 Activation (OR = 11.8, p = 0.0003) (Supplementary Table S2).” These findings hold substantial importance as they are closely tied to cell signaling cascades crucial for the activation of the Mitogen-Activated Protein Kinase (MAPK) pathway. The centrality of the MAPK/ERK pathway cannot be overstated, as it plays a pivotal role in orchestrating the intricate process of ureteric bud (UB) branching morphogenesis—a fundamental aspect determining the kidney’s size, structure, and the composition of its renal units (Kurtzeborn et al., 2019). Additionally, the tyrosine phosphatase SHP2, encoded by the PTPN11 gene, is essential for MAPK cascade activation (Roberts et al., 2007). A prior study by Du and colleagues demonstrated that Robo2 and Gen1 jointly regulate UB initiation, significantly increasing the incidence of CAKUT. Importantly, the MAPK/ERK pathway serves as a central mechanism for these two genes (Du et al., 2023). In line with Du’s findings, our study suggests that PTPN11 may potentially influence UB development leading to CAKUT through the aforementioned three putative cell signaling pathways, in addition to modulating the MAPK/ERK pathway.

This discovery illuminates the intricate interplay between these pathways and the development of both the spine and the kidneys, providing valuable insights into the complex molecular mechanisms governing the morphogenesis of these vital structures. For a more intuitive understanding of these pivotal pathways, please refer to Figures 3A, B, which offer compelling insights into the potential keys for unraveling the complexities of spine and kidney development and function.

To obtain a wide genomic understanding of the genetic pathogenisis of CS-CAKUT, we performed association analysis of CNVs that are most likely to disrupt gene function.

Our findings revealed significant copy number loss signals in NTRK3 (p = 1.29E-05) and KMT2B (p = 3.26E-06) (Supplementary Table S3). NTRK3 encodes the neurotrophic tyrosine receptor kinase (NTRK), and KMT2B is responsible for encoding histone H3 methyltransferase. Both genes have the potential to influence intracellular phosphorylation signaling pathways, regulating crucial protein phosphorylation and dephosphorylation events. These pathways are vital across various biological processes and have implications for conditions such as childhood muscle tone disorders and congenital mesoblastic nephroma (CMN) (Knezevich et al., 1998; Meyer et al., 2017).

Additionally, copy number loss of IGFLR1 (p = 3.26E-06) emerges as a noteworthy result (Supplementary Table S3). IGFLR1, also known as Transmembrane Protein 149 (TMEM149), is a protein-encoding gene located on chromosome 19. It is widely expressed in lymph nodes, spleen, and kidneys and has been implicated in promoting the progression of renal cell carcinoma (Fagerberg et al., 2014; Song et al., 2020). Despite its relatively limited attention in previous research, this uncharacterized gene possesses structural characteristics akin to the tumor necrosis factor receptor family. Tumor necrosis factor, in turn, plays a significant role in spinal deformities (Kudo et al., 2018). Our analysis has suggested that haploinsufficiency of IGFLR1 might hold a close association with the pathogenesis of CS-CAKUT, shedding new light on a previously underexplored genetic component within this complex syndrome.