Rare variant enrichment analysis supports GREB1L as a contributory driver gene in the etiology of Mayer-Rokitansky-Küster-Hauser syndrome

Abstract

Mayer-Rokitansky-Küster-Hauser (MRKH) syndrome is characterized by aplasia of the female reproductive tract; the syndrome can include renal anomalies, absence or dysgenesis, and skeletal anomalies. While functional models have elucidated several candidate genes, only WNT4 (MIM: 603490) variants have been definitively associated with a subtype of MRKH with hyperandrogenism (MIM: 158330). DNA from 148 clinically diagnosed MRKH probands across 144 unrelated families and available family members from North America, Europe, and South America were exome sequenced (ES) and by family-based genomics analyzed for rare likely deleterious variants. A replication cohort consisting of 442 Han Chinese individuals with MRKH was used to further reproduce GREB1L findings in diverse genetic backgrounds. Proband and OMIM phenotypes annotated using the Human Phenotype Ontology were analyzed to quantitatively delineate the phenotypic spectrum associated with GREB1L variant alleles found in our MRKH cohort and those previously published. This study reports 18 novel GREB1L variant alleles, 16 within a multiethnic MRKH cohort and two within a congenital scoliosis cohort. Cohort-wide analyses for a burden of rare variants within a single gene identified likely damaging variants in GREB1L (MIM: 617782), a known disease gene for renal hypoplasia and uterine abnormalities (MIM: 617805), in 16 of 590 MRKH probands. GREB1L variant alleles, including a CNV null allele, were found in 8 MRKH type 1 probands and 8 MRKH type II probands. This study used quantitative phenotypic analyses in a worldwide multiethnic cohort to identify and strengthen the association of GREB1L to isolated uterine agenesis (MRKH type I) and syndromic MRKH type II.

Keywords: Alu-Alu mediated rearrangement; developmental genomics; exonic deletion CNV; genitourinary development; human phenotype ontology; infertility; rare variant gene enrichment; sex-limited trait; transcriptional regulation.

Figure 1

Genes with significant rare variant mutational enrichment in the Americas and European MRKH cohort The y axis represents the -log of the adjusted p value and x axis the ratio of the frequency of rare variants passing parsing steps in the MRKH cohort to the control cohort. p value testing was done using a Fisher’s exact test and p values were adjusted using the Holm-Bonferroni method. A total of 3,003 genes, which were identified to have rare variants present in the MRKH cohort using the parsing criteria described (materials and methods), were tested. Of these, 3,001 genes had an adjusted p value of 1.

Figure 2

Pedigrees representing families with predicted deleterious GREB1L variant alleles observed in the MRKH cohort Identifiers for probands are located to the left of each respective proband in the pedigree. (A) Sanger dideoxy traces for each variant allele from individuals with DNA available for testing are represented below the respective pedigree symbol. (B) Variant data are represented in a table at the bottom of the figure. All variants refer to the GRCh37/Hg19 genome assembly build and to RefSeq transcript ID NM_001142966.3. AF, allele frequency; BHCMG, Baylor-Hopkins Center for Mendelian Genomics Database; CADD, combined annotation dependent depletion score; CHR, chromosome; POS, position; tR, total reads; vR, variant reads; Zyg, zygosity.

Figure 3

GREB1L probands cluster within a phenotypic group including uterine and renal malformations Heatmap represents symmetric Lin score for proband HPO terms queried against each other. Clusters are demarcated by black outline. Top left cluster: Includes probands (cluster #1) that share unilateral renal agenesis as a feature. Bottom right cluster: Includes probands that share scoliosis (top left group of bottom right cluster; i.e., cluster #2) or ectopic kidney (bottom right group of bottom right cluster; cluster #3). GREB1L probands are represented by red font.

Figure 4

GREB1L deletion CNV in proband BH9101-1 (A) Droplet digital PCR results for BH9101-1 and NA15501 (an unaffected control sample) for exon 18 of GREB1L and a control biallelic gene, RPP30, are shown. Event number is displayed on the y axis and BH9101-1 and NA15501 results for GREB1L and RPP30 are displayed on the x axis. Exon 18, denoted by a red star throughout Figure 4, ddPCR of BH9101-1 resulted in significantly fewer positive droplets (events) than controls. (B) Visualization of the read ratio of BH9101-1 and BH9105-1 (a control from the MRKH cohort without a deletion called) raw reads to the average raw reads across 100 samples with the most similar read profile across the exome to BH9101-1 or BH9105-1, and raw read count of BH9101-1 and BH9105-1 across GREB1L and ESCO1 exons called by XHMM to be deleted and the gene closest to the 3′ end of GREB1l, ROCK1. The exon structure of ROCK1, GREB1L, and ESCO1 are shown at the bottom. Raw read count for each exon is shown in red for BH9101-1 and green for BH9105-1. (C) HMZDelFinder visualization of RPKM across the coding exons of GREB1L and flanking genes is shown in red for BH9101-1, green for BH9105-1, and black for the 100 samples with the most similar read profile to BH9101-1. RPKM across exons 3–33 of GREB1L and 4–12 of ESCO1—representing all coding exons of both genes—is significantly reduced in BH9101-1. Log₁₀(RPKM +1) is shown on the y axis and exon structure and number of exons shown across the x axis. Arrowheads show the direction of transcription. The region with a called CNV is represented by a red bar along the x axis. Below the HMZDelFinder visualization, the UCSC Genome Browser of the region corresponding to deletion is shown. UCSC and RefSeq gene transcripts are shown above and repetitive elements are depicted at bottom. The scale at top represents 100kb. (D) Top: Genome-wide 1M Probe aCGH results for BH9101-1 are shown for the region of chromosome 18 with deletion. Minimum and maximum deletion intervals are shown based on the last adjacent probes with a deletion call and normal copy number state on 3′ and 5′ ends of the deletion. Log₂ ratio is shown on the y axis and chromosome position on the x axis. Genes within the region are shown at bottom and the deletion region highlighted by a green bar at top. Bottom: Visualization of Illumina short-read whole genome sequencing read depth data in BH9101-1 was performed using VizCNV (https://github.com/BCM-Lupskilab/VizCNV). A red line depicts the average log₂ ratio of read depth at a locus to the median read depth of chromosome 18 (y axis) within the genomic region chr18:20,999,991-21,996,196 (GRCh38; x axis). Black dots represent the log₂ ratio at each locus with a diploid copy number state. Green dots represent the log₂ ratio at each locus within the deletion CNV. (E) Breakpoint junction Sanger sequencing chromatogram for BH9101-1 is shown at left. The 18-bp microhomology sequence is denoted by a black rectangle. On the right, a simplified genomic view of the deletion CNV encompassing GREB1L coding exons (3–33, shown in purple) and ESCO1 coding exons (4–12, shown in orange) is depicted. The Alu elements AluY and AluSx are shown in red and green, respectively. A representation of the AAMR generated deletion and recombinant chimeric Alu (red/green hybrid) is shown below the simplified genomic view, with the individual Alu involved in the recombination shown in red and green, respectively. The deletion spans 302,966 bp (chr18:18,856,932-19,159,898). At bottom right, the sequences of the proximal (red) and distal breakpoints (green) are shown above and below the AAMR product sequence (red and green). Microhomology sequence is shown in blue and denoted by a black rectangle. Genomic coordinates for the proximal breakpoint are shown above the proximal breakpoint sequence, and those for the distal breakpoint are shown below the distal breakpoint sequence.

Figure 5

Rare variants of GREB1L are associated with female reproductive tract, urinary tract, skeletal/digital, and ear/CN VII malformations (A) Venn diagram depicting the number and overlap of variants associated with each phenotypic trait. Purple denotes female reproductive tract malformations, orange urinary tract malformations, blue skeletal/digital malformations, and yellow ear/CN VIII malformations. (B) A representation of disease-associated variant type, with LoF/PTV variants colored red, missense blue, and nonframeshift green, within coding exons 3–33 of GREB1L. Phenotypes associated with variants are depicted in a heatmap below variant type diagram, and orange shading within a row marks association of a variant with a phenotype trait. Exons 6, 16, 19, 21, 23, 26, 27, and 31 are separated from adjacent exons in the heatmap to better delineate phenotype overlap as these exons encompass four or more disease-associated variants. CNV deletion GREB1L null allele not shown.

Figure 6

Phenotypic grid of GREB1L variant-associated HPO terms Left: GREB1L pathogenic variant alleles are labeled at left, and HPO phenotype terms labeled at the bottom of the phenotype grid. Red denotes the presence of a term and gray absence. Terms are ordered by frequency, which is displayed by black bars on a 0–1 scale at the top of the heatmap. Asterisk (right of grid) represents the specific row for the GREB1L deletion CNV null allele. Right: The phenotypic grid for GREB1L probands has been collapsed into broader phenotypic categories labeled at the bottom of the heatmap. Frequency is displayed at top by black bars on a 0–1 scale. Asterisk (right of grid) represents the specific row for the GREB1L deletion CNV null allele.