Diverse ancestry whole-genome sequencing association study identifies TBX5 and PTK7 as susceptibility genes for posterior urethral valves

Abstract

Posterior urethral valves (PUV) are the commonest cause of end-stage renal disease in children, but the genetic architecture of this rare disorder remains unknown. We performed a sequencing-based genome-wide association study (seqGWAS) in 132 unrelated male PUV cases and 23,727 controls of diverse ancestry, identifying statistically significant associations with common variants at 12q24.21 (p=7.8 × 10^-12; OR 0.4) and rare variants at 6p21.1 (p=2.0 × 10^-8; OR 7.2), that were replicated in an independent European cohort of 395 cases and 4151 controls. Fine mapping and functional genomic data mapped these loci to the transcription factor TBX5 and planar cell polarity gene PTK7, respectively, the encoded proteins of which were detected in the developing urinary tract of human embryos. We also observed enrichment of rare structural variation intersecting with candidate cis-regulatory elements, particularly inversions predicted to affect chromatin looping (p=3.1 × 10^-5). These findings represent the first robust genetic associations of PUV, providing novel insights into the underlying biology of this poorly understood disorder and demonstrate how a diverse ancestry seqGWAS can be used for disease locus discovery in a rare disease.

Keywords: CAKUT; PUV; developmental biology; genetics; genome-wide association study; genomics; human; posterior urethral valves; seqGWAS; whole-genome sequencing.

Figure 1.. Study workflow.

The flowchart shows the number of samples included at each stage of filtering, the analytical strategies employed, and the main findings (blue boxes). PUV, posterior urethral valves; MAF, minor allele frequency; GWAS, genome-wide association study; EUR, European; cCRE, candidate cis-regulatory element.

Figure 1—figure supplement 1.. Principal component analysis (PCA) showing the first eight principal components for matched cases (blue) and controls (black) and unmatched cases (orange) and controls (grey).

Two cases and 2579 controls were excluded from downstream analyses.

Figure 2.. Manhattan plot for mixed-ancestry sequencing-based genome-wide association study (seqGWAS).

A genome-wide single-variant association study was carried out in 132 unrelated posterior urethral valves (PUV) cases and 23,727 controls for 19,651,224 variants with minor allele frequency (MAF) >0.1%. Chromosomal position (GRCh38) is denoted along the x axis and strength of association using a –log₁₀(p) scale on the y axis. Each dot represents a variant. The red line indicates the Bonferroni adjusted threshold for genome-wide significance (p<5 × 10^–8). The gene in closest proximity to the lead variant at significant loci is listed.

Figure 2—figure supplement 1.. Quantile-quantile (Q–Q) plot for the mixed-ancestry genome-wide association study (GWAS) displaying the observed vs. the expected –log₁₀(p) for each variant tested.

The grey shaded area represents the 95% confidence interval of the null distribution.

Figure 2—figure supplement 2.. Power calculations for the mixed-ancestry genome-wide association study (GWAS) were performed at various minor allele frequencies (MAF) using 132 cases and 23,727 controls under an additive genetic model to achieve genome-wide significance of p<5 × 10^–8.

Figure 3.. 12q24.21.

Regional association plot with chromosomal position (GRCh38) denoted along the x axis and strength of association using a –log₁₀(p) scale on the y axis. The lead variant (rs10774740) is represented by a purple diamond. Variants are coloured based on their linkage disequilibrium (LD) with the lead variant using 1000 Genomes data from all population groups. Functional annotation of the lead prioritized variant rs10774740 is shown, intersecting with CADD score (version 1.6), PhastCons conserved elements from 100 vertebrates, and ENCODE H3K27ac ChIP-seq, H3K4me3 ChIP-seq, and DNase-seq from mesendoderm cells. ENCODE cCREs active in mesendoderm are represented by shaded boxes; low DNase (grey), DNase-only (green). Genome-wide association study (GWAS) variants with p<0.05 are shown. Note that rs10774740 has a relatively high CADD score for a non-coding variant and intersects with a highly conserved region. PP, posterior probability derived using PAINTOR; cCRE, candidate cis-regulatory element.

Figure 3—figure supplement 1.. Heatmap of Hi-C interactions from H1 BMP4-derived mesendoderm cells demonstrating that rs10774740 is located within the same topologically associating domain (TAD) as TBX5.

TADs are represented by blue triangles. Protein-coding genes are denoted in blue, non-coding genes in green.

Figure 3—figure supplement 2.. Circos plot illustrating significant chromatin interactions between 12q24.21 and the promoter of TBX5.

The outer layer represents a Manhattan plot with variants plotted against strength of association. Only variants with p<0.05 are displayed. Genomic risk loci are highlighted in blue in the second layer. Significant chromatin loops detected in H1 BMP4-derived mesendoderm cultured cells are represented in orange.

Figure 4.. 6p21.1.

Regional association plot with chromosomal position (GRCh38) along the x axis and strength of association using a –log₁₀(p) scale on the y axis. The lead variant (rs144171242) is represented by a purple diamond. Variants are coloured based on their linkage disequilibrium (LD) with the lead variant using 1000 Genomes data from all population groups. Functional annotation of the lead prioritized variant rs144171242 is shown intersecting with ENCODE H3K27ac ChIP-seq, H3K4me3 ChIP-seq, and DNase-seq from mesendoderm cells. ENCODE cCREs active in mesendoderm are represented by shaded boxes; low DNase (grey), DNase-only (green), and distal enhancer-like (orange). ChromHMM illustrates predicted chromatin states using Roadmap Epigenomics imputed 25-state model for mesendoderm cells; active enhancer (orange), weak enhancer (yellow), strong transcription (green), transcribed and weak enhancer (lime green). Predicted transcription factor-binding sites (TFBS) from the JASPAR 2020 CORE collection (Fornes et al., 2020) are indicated by dark grey shaded boxes. Genome-wide association study (GWAS) variants with p<0.05 are shown. Note that rs144171242 intersects with both a predicted regulatory region and TFBS. PP, posterior probability derived using PAINTOR; cCREs, candidate cis-regulatory elements.

Figure 4—figure supplement 1.. Sequence logos representing the DNA-binding motifs of transcription factors FERD3L and ZNF317.

The black boxes indicate where the risk allele [G] may disrupt binding.

Figure 5.. Principal component analysis for the replication cohort.

Principal component analysis showing the first two principal components for a subset of cases (red) for whom genome-wide genotyping data was available (n=204), and the control (grey) cohort from 100,000 Genomes Project (100KGP) (n=4151) projected onto samples from the 1000 Genomes Project (Phase 3). Both cases and controls had confirmed European ancestry.

Figure 6.. Manhattan plot for European sequencing-based genome-wide association study (seqGWAS).

A genome-wide single-variant association study was carried out in 88 cases and 17,993 controls for 16,938,500 variants with MAF ≥0.1%. All cases and controls had genetically determined European ancestry. Chromosomal position (GRCh38) is denoted along the x axis and strength of association using a –log₁₀(p) scale on the y axis. Each dot represents a variant. The red line indicates the Bonferroni adjusted threshold for genome-wide significance (p<5 × 10^–8). The gene in closest proximity to the lead variant at significant loci are listed.

Figure 6—figure supplement 1.. Quantile-quantile (Q-Q) plot displaying the observed vs. the expected –log₁₀(p) for each variant tested.

The grey shaded area represents the 95% confidence interval of the null distribution.

Figure 6—figure supplement 2.. Comparison of (A) −log₁₀(p) and (B) BETA from the diverse ancestry and European-only genome-wide association study (GWAS).

All variants with p<10^–5 in both cohorts are shown. The shaded grey area represents the 95% confidence interval.

Figure 6—figure supplement 3.. Ancestry-specific minor allele frequencies for (A) rs10774740 (T) at 12q24.21 and (B) rs144171242 (G) at 6p21.1.

Error bars represent 95% confidence intervals. The lead variant in (B) was not identified in individuals with African ancestry in this study. AFR, African ancestry (11 cases; 483 controls); EUR, European ancestry (89 cases; 17,993 controls); SAS, South Asian ancestry (18 cases; 2948 controls).

Figure 6—figure supplement 4.. Forest plots demonstrating ancestry-specific odds ratios for (A) rs10774740 (T) and (B) rs144171242 (G).

Error bars represent 95% confidence intervals. The lead variant in (B) was not identified in individuals with African ancestry in this study. AFR, African ancestry (11 cases; 483 controls); EUR, European ancestry (89 cases; 17,993 controls); SAS, South Asian ancestry (18 cases; 2948 controls); ALL, mixed-ancestry cohort (132 cases; 23,727 controls).

Figure 6—figure supplement 5.. Linkage disequilibrium (LD) plots for 503 European (EUR), 489 South Asian (SAS), and 661 African (AFR) ancestry individuals from the 1000 Genomes Project (Phase 3).

Haploview (version 4.2) was used to compute pairwise LD statistics (r²) between variants for each population. The darker the shading, the higher the LD between variants. Black outlined triangles indicate haploblocks. (A) LD plot for chr12:114,663,967–114,667,916 (GRCh37) with the position of the lead variant rs10774740 represented by a green arrow; (B) LD plot for chr6:43,084,099–43,092,650 (GRCh37) with the lead variant rs144171242 represented by a green arrow. rs144171242 was not seen in the AFR population group.

Figure 7.. Manhattan plot of PheWAS for rs10774740 (T) at 12q24.21.

Plot downloaded from https://pheweb.org/UKB-SAIGE. The PheWAS was performed using imputed data from ~400,000 White British participants in the UK Biobank using SAIGE. The triangles indicate the direction of effect. The dashed grey line indicates a Bonferroni adjusted significance level of p< 3.6 x 10^-5 (1403 phenotype codes).

Figure 8.. Immunohistochemistry in human embryogenesis.

(A) Overview of transverse section of a normal human embryo 7 weeks after fertilization. The section has been stained with haematoxylin (blue nuclei). Boxes around the urogenital sinus and the mesonephric duct mark similar areas depicted under high power in (B–E). In (B–D), sections were reacted with primary antibodies, as indicated; in (E) the primary antibody was omitted. (B–E) were counterstained with haematoxylin. In (B–E), the left-hand frame shows the region around the mesonephric duct, while the right-hand frame shows one lateral horn of the urogenital sinus. (B) Uroplakin 1b immunostaining revealed positive signal (brown) in the apical aspect of epithelia lining the urogenital sinus (arrows, right frame), the precursor of the urinary bladder and proximal urethra. Uroplakin 1b was also detected in the flat monolayer of mesothelial cells (left frame) that line the body cavity above the mesonephric duct. (C) There were strong PTK7 signals (brown cytoplasmic staining) in stromal-like cells around the mesonephric duct (left frame), whereas the epithelia of the duct itself were negative. PTK7 was also detected in a reticular pattern in epithelia lining the urogenital sinus (right frame) and in stromal cells near the sinus. (D) A subset of epithelial cells lining the urogenital sinus (right frame) immunostained for TBX5 (brown nuclei; some are arrowed). The mesothelial cells near the mesonephric duct (left frame) were also positive for TBX5. (E) This negative control section had the primary antibody omitted; no specific (brown) signal was noted. Bar is 400 μm in (A), and bars are 100 μm in (B–E). ugs, urogenital sinus; md, mesonephric duct; hg, hindgut; u, ureter.

Figure 8—figure supplement 1.. Immunohistochemistry of a second 7-week human embryo counterstained with haematoxylin.

(A) View of the mesonephric duct (the epithelial tube with md in its lumen). Note the prominent signal (brown) for PTK7 in the stromal cells surrounding the duct. (B) View of the urogenital sinus (ugs) with a subset of nuclei (three shown by arrows) in its monolayer epithelium that stain (light brown) for the transcription factor TBX5. The hindgut (hg) is nearby. Bars are 100 μm.

Figure 9.. Manhattan plot of exome-wide rare coding variant analysis.

Chromosomal position (GRCh38) is shown on the x axis and strength of association using a –log₁₀(p) scale on the y axis. Each dot represents a gene. The red line indicates the Bonferroni adjusted threshold for exome-wide significance (p=2.58 × 10^–6). Genes with p<10^–4 are labelled.

Figure 9—figure supplement 1.. Quantile-quantile (Q-Q) plot displaying the observed vs. the expected –log₁₀(p) for each gene tested.

The grey shaded area represents the 95% confidence interval of the null distribution.

Figure 10.. Rare structural variant burden analysis.

The proportion of individuals with ≥1 rare autosomal structural variant intersecting with an ENCODE candidate cis-regulatory element (cCRE) in cases and controls was enumerated using a two-sided Fisher’s exact test. Note that inversions affecting cCRE are enriched in PUV. Vertical black bars indicate 95% confidence intervals. Unadjusted p-values shown are significant after correction for multiple testing (p<2.5 × 10^–3). CNV, copy number variant; DEL, deletion; DUP, duplication; INV, inversion; PUV, posterior urethral valves; dELS, distal enhancer-like signature; pELS, proximal enhancer-like signature; PLS, promoter-like signature; cCRE, candidate cis-regulatory element.