Detection of de novo genetic variants in Mayer-Rokitansky-Küster-Hauser syndrome by whole genome sequencing.

OBJECTIVE: The aim of this study was to use whole genome sequencing (WGS) help detect de novo mutations or pathogenic genes of Mayer-Rokitansky-Küster-Hauser syndrome type 1(MRKH syndrome type 1). STUDY
DESIGN: This was a case-parent trios study. Nine unrelated probands, with MRKH syndrome type 1 and their parents were enrolled. The enrollment, sequencing process, establishment of the de novo mutations detecting procedure and experiment part were performed over a 2-year period.
RESULTS: we detected 632 de novo single nucleotide variants (SNVs), 267 de novo small insertions/deletions (indels), 39 de novo structural variations (SVs) and 28 de novo copy number alterations (CNAs). Three novel damaging coding de novo SNVs with three damaging coding de novo genes (PIK3CD, SLC4A10 and TNK2) were revealed. Two SNVs were annotated of the promoter region of gene NBPF10 and 3'UTR of NOTCH2NL, potentially contributing to the pathogenesis of MRKH.
CONCLUSION: We identified five de novo mutations in BAZ2B, KLHL18, PIK3CD, SLC4A10 and TNK2 by performing WGS, the functional involvement of all deleterious mutations in MRKH candidate genes of the trios warrant further study. WGS may complement conventional array to capture the complete landscape of the genome in MRKH.

Introduction

Mayer–Rokitansky–Küster–Hauser (MRKH) syndrome (OMIM # 277000) is characterized by congenital hypoplasia of the uterus, cervix and upper two-thirds of the vagina during fetal development [1]. MRKH syndrome is a rare birth defect that affects approximately 1 in 4,500–5,000 newborn girls [2]. Typically, individuals with MRKH syndrome have normal ovarian function and a normal female karyotype (46,XX) [3]. MRKH syndrome can be classified into type 1 (isolated; OMIM 277000), and type 2 (associated with renal, skeletal, cardiovascular and other malformations). MRKH Type 2 is also known as MURCS (Mullerian Duct Aplasia, Unilateral Renal Agenesis, And Cervicothoracic Somite Anomalies (OMIM 601076).

MRKH has been reported to be a sporadic malformation with characteristic familial clustering and a high risk of associated malformations in trios of affected women. The etiology of MRKH syndrome remains unknown. Several genes have been linked to MRKH, such as HOXA7 [4], HOXA13, PBX1 [5], AMH [6], RBM8A andTBX6 [6]. These studies have failed to provide sufficient insight into MRKH syndrome and the primary developmental anomalies. De novo mutations play an important role in human disease [7]. Furthermore there are few genetic studies exploring the possible genetic etiology of MRKH syndrome [8]. Whole exome and whole genome studies are two approaches to identify de novo mutations in candidate genes for particular developmental conditions. This approach has been used extensively in uncovering potential new genes involved in autism spectrum disorder [9]. A recent study using whole-exome sequencing (WES) combined with SNP arrays identified five novel mutated genes among seven patients with MRKH type 1 [10]. These mutations were recurrent in more than two individuals, but the study does not include analysis of the corresponding parent samples.

In this study, we hypothesize that the MRKH type 1 may be due to the occurrence of de novo mutations. We used the first whole genome sequencing (WGS) in full parent-child MRKH type 1 trio families to unequivocally identify these events.

Material and methods

Samples from MRKH syndrome-affected families

Nine unrelated Chinese patients with MRKH Type 1 and their parents were recruited from March 2014 to November 2014 in the Department of Obstetrics and Gynecology of Luohu Hospital. Blood samples and phenotypic data were collected from probands and their parents (Table S1). All participants provided written consent. Approval for the study was obtained from the Ethics Committees of Luohu Hospital.

Sequencing, alignment and variant discovery

Genomic DNA of probands and their parents (trio families) were extracted from blood and sequenced using HiSeq X Ten sequencing system (500 bps library, 150 bps paired reads). The genome mean depth was at least 50× each sample (probands: ˜57×, maternal: ˜56× and paternal: ˜53×), and the coverage of each sample was > 99% (Fig. S2 and Table S2). Reading were aligned with the reference genome (hg19, GRCh37) by the Burrows-Wheeler Aligner, BWA, version 0.59 [11]. Sequence Alignment/Map package (SAM tools [12]: Picard tools (version 1.52)) was used to remove duplications. The single nucleotide variants (SNVs) and small insertions/deletions (Indels) were identified using the Genome Analysis Toolkit, GATK version 3.4 [13]. All SNV/Indels mutation events were presented in Table S3. The effects (e.g., missense, nonsense, or frameshift mutations) and classifications (e.g., in exonic, intronic, or intergenic regions) of variants across the genome were annotated by ANNODB (unpublished software) and Oncotator [14]. The flowchart of the analysis protocol is shown in Fig. S1.

De novo SNV detection

A potential de novo SNV mutation (DNM) was identified following both of the two conditions: (1) one site is heterozygous (novel allele frequency >30%) or homozygous with non-reference SNV in the proband. (2) This site is homozygous (novel allele frequency <5%) with reference nucleotide in both parents. In addition, other filters: (1) the read depth (DP) at the variant position of each member of the trio families should be more than 10×. (2) a phred-scaled quality score(PSQ) more than 30 [15]. (3) a QualByDepth (QD: variant confidence from the QUAL field/unfiltered depth) of more than 10. (4) a HomopolymerRun (Hrun: Largest contiguous homopolymer run of the variant allele in either direction on the reference) less than 5. (5) a MappingQualityZero (MQ0: total count across all samples that had reads with a mapping quality of zero) less than 4.

De novo indel detection

A similar filtering protocol was applied to detect de novo indels: the proband should be heterozygous or homozygous at the indel site, whereas the both parents should be homozygous without indels at the same position [15]. Furthermore, other filters: DP > 10 at the variant position of each member of the trio families; PSQ more than 30 [15]; QD > 10; Hrun less than 5 and MQ0 less than 4.

Variant experimental validation

To validate genomic loci genotypes in offspring and parents was using the Sanger sequencing. A total of eight of nine in exon region putative DNMs sites were validated. In subsequent exonic analyses, we combine sites with complete validation data of de novo Mutations (Table S7).

Variant characterization

When assessing inherited variants, we removed those that are common in the population (i.e., found in dbSNP135). Deleterious variants obtained from nonsense, splice site and frameshift. We defined damaging missense SNVs as those that both Sorting Intolerant From Tolerant (SIFT) [16] and PolyPhen-2 [17] predicted to be damaging via the variant effect predictor(Table S8).

De novo CNA detection

The SegSeq [18] was used to detect copy number alterations (CNAs). Proband-specific CNAs, that did satisfy the following criteria would be discarded: (1) q0 (the fraction of reads with zero mapping quality in the called CNA regions) < 50%. (2) Corrected P-value < 0.05. (3) Overlap rate of repeat database (STR: simple tandem repeat, DGV, genomic SuperDups) < 50%. Putative de novo CNAs for each trio family were detected by the following filters: (1) the normalized average RD signal in the child was < 1.4 (for deletions) or >2.6 (for duplications); and (2) the normalized average RD signal in each parent was between1.6 and 2.4.

To infer recurrently amplified or deleted genomic regions, we re-implemented the GISTIC 2.0 algorithm [19], using copy numbers in 100-kb windows as markers. G-scores were calculated for genomic regions on the basis of the frequency and amplitude of amplifications or deletions affecting each gene. A significant CNA region was defined as having amplification/deletion with |G-score| > 0.1, corresponding to a P-value threshold of 0.01 from the permutation-derived null distribution.

De novo structural variations (SVs) detection

Based on the reliability of the algorithm, the Meerkat 0.185 [20] was used with suggested parameters for de novo SVs and breakpoints detection of 9 trio families. In summary, we mapped whole reads against the human reference genome (hg19) to find soft-clipped and unmapped reads. Then re-mapped these reads to identify discordant read pairs and searched for reads that were covered the candidate breakpoints, and refined precise breakpoints by local alignments. De novo SVs were generated by filtering out SVs events in parents, and appear only in the offspring. Only high confidence events were used for downstream analysis.

Network and pathway analysis

Projecting data onto network contexts is a powerful way to unravel patterns embedded in seemingly, scattered large data sets to amplify knowledge discovery related to complex diseases [21]. We conducted a comprehensive literature review of all variant genes with de novo SNVs variants and searched for the connections between the genes.

The Functional Interaction (FI) network of these genes was generated with the use of the application ReactomeFIViz [21] on the Cystoscope platform [22]. The genes focus of the FI Network by an spectral partition based network clustering algorithm [23]. Nodes in different network modules show in different colors. All white node genes of this Network were supported by at least one reference from the literature or in our study, other gray genes were through calculation of the algorithm, all of these genes were associated with a relevant pathway in the Pathway Database (Fig. S3).

Result

Genetic variants in WGSV

We performed whole genome next generation sequencing (WGS) on 9 trio families. The mean sequencing depth was around 50× (46.48× - 71.22×; Fig. S2) and the mean coverage was more than 99% (Table S2). Totally, we detected an average of 3,909,478 SNPs and 797,951 indels (Table S3) per sample. The mutation rate was approximately 1% [24], which is similar with the results of 1000 Genome Project [25], UK10 K project [26] and ASD project [7].

De novo mutations

We detected 632 de novo SNVs (70.2 per individual), 267 de novo indels (29.7 per individual), 28 de novo CNAs (3.1 per individual; Table S9) and 39 de novo SVs (4.4 per individual; Table S12). 9 de novo SNVs were located in exonic regions (1 mutation per individual; Table S13), and 7 of these were an amino acid change mutations (one was in Dbsnp135) and the other 2 were synonymous. Functional analysis with both PolyPhen-2 and SIFT predicted that 4 of the 6 missense mutations would be deleterious (Table S8).

267 highly confident de novo INDELs were retained for further analysis (Tables S4 and S5).

MRKH 6 had a missense mutation (G1535A; R512Q) in PIK3CD (OMIM: 602839; RefSeq mRNAs: NM_005026). MRKH 4 had a The c.1706C > G mutation in TNK2. MRKH 9 had a missense de novo mutation (C709 T; R237C) in SLC4A10 (OMIM: 605556; RefSeq mRNAs: NM_001178015).

De novo SNVs in non-coding regions

In non-coding regions, we identified 623 de novo SNVs, most of which were singletons and newly discovered based on the ENCODE database (Tables S4 and S5).We found 168 variants that mapped to DNase1 hypersensitive sites (DHS: 99), transcription factor binding sites (TFP: 52), enhancer elements [14], and Pseudogene [3]. 29 of 168 mutations were detected as recurrent mutations (Table S6). One mutation (chr1: 145291364A > T) presented in MRKH 1 was located in 3'UTR of NOTCH2NL. The region is also the promoter of NBPF10 (Table 2). A de novo mutations (chr1: 145291369 G > A) was identified in the 3'UTR of NOTCH2NL in MRKH 8.

Table 2

The noncoding. recurrent NBPF10 by whole-genome sequencing (WGS) analysis.

No.	Target gene	Noncoding. recurrent	Family ID
1	NBPF10 (Promoter)	TFP(BATF|chr1:145291202-145292202)	Trio01(chr1:145291364)	Trio08(chr1:145291369)	Trio09(chr1:145292171)
		TFP(NFKB1|chr1:145290201-145291965)	Trio01(chr1:145291364)	Trio08(chr1:145291369)	Trio10(chr1:145290279)

*TFP -transcription factor binding peak.

De novo SVs

A total of 39 de novo SVs were identified in all 9 trio families. Statistic data and a full list of de novo SVs are presented in Tables S10, S11 and S12. Six types of SVs were observed, including deletions insertions(del_ins), tandem duplications (tandem_dups), inversions_reverse(invers_r), inter-chromosomal translocations(transl_inter), deletions, insertions, and inversions_follow (invers_f). Of these SVs, 31% were deletions, about 15% were inter-chromosomal translocations(transl_inter). For these 39 SVs, fork stalling and template switching (FoSTeS), non-homologous end joining (NHEJ) and alternative end joining(alt-EJ) were the dominant mechanisms.

In 9 trio families’ genomes, we found that 15 SVs occurred in type of gene-gene and were predicted as a direct link to a gene such as DUSP16 (dual specificity phosphatase 16), KIF6(kinesin family member 6), PRKCZ(protein kinase C zeta), and ZNF665(zinc finger protein 665), and 19 genes contained a breakpoint in one of these trio families. Specifically, one breakpoint in chr17 is effect to the target gene: TBC1D16, (Intron &Promoter) (Table S10) (Fig. 1, inter-chromosomal translocations (chr17-chrX)).

Fig. 1

De novo structural variations.

Copy-number alterations (CNAs)

We profiled the 9 trio families for CNAs using WGS and found abnormalities of chromosomal regions that identified a total of 28 de novo CNAs. A full list of CNAs is presented in Table S9. These CNAs occurs on different chromosomes, including chr1, chr2, chr3, chr7, chr8, chr10, chr11, chr13, chr14, chr15, chr16, chr17, chr20, chr22 and the region of CNA in exonic is ranging from 4k to 200k.

Profiling of CNAs identified many putative disease related genes, that may be implicated in MRKH. We applied GSITIC analysis to WGS data to identify recurrent focal CNAs, and this analysis yielded 23 regions of focal amplification in proband-Mother and 18 regions of focal amplification in proband-Father (Table S14), several of which included CNAs previously detected in MRKH type 1, which encompassed genes such as OR2T2, LOC727924, OR4N4, OR4M2, PSG9.

A significant discovery, was the frequent amplification of PSG9, pregnancy specific beta-1 glycoprotein 9, which is gained in proband-mother (4/9) and lost in proband-father (5/9). (Figs. S4 and S5).

We identified 38 focal deletion regions in proband-Mother and 48 regions in proband-Father, which were also frequently truncated in MRKH. One of the most common focal deletions, detected in 9 trio familiess (90%), was a deleted region at 15q11.2 containing LOC727924, OR4N4 and OR4M2 (Fig. 2) and is a previously reported genomic alteration in MRKH.

Fig. 2

Landscape of CNV.

Inherited dominant mutations in autosomal genes

With analysis of Autosomal Genes with inherited dominant mutations, one gene in exon was yielded, and the gene named KRTAP5-10, Keratin Associated Protein 5-10. In non-coding regions, We found 3 recurrent variants that mapped to enhancer elements, with target gene that DUOX1(Dual Oxidase 1), MLLT4(Myeloid/Lymphoid Or Mixed-Lineage Leukemia; Translocated To, 4) and RPS6KA2(Ribosomal Protein S6 Kinase A2).These functional data indicated that these 3 genes might plays an important role in MRKH syndrome (Tables S15 and S16).

In this study, only genes interacting with one pathway gene show to minimize a diagram of clutter (Fig. S2). Newly added genes from pathway database color in grey, while original genes color in white. FIs extracted from pathways show as solid lines (for example, EGFR-PIK3CD, while those predicted based on NBC show as dashed lines (for example, TNK2-EGFR). These Extracted FIs involved in activation might play an important role possibly in MRKH syndrome.

Discussion

We detected eight de novo SNVs in nine families, in which three deleterious variants were validated. Our study also identified a number of previously unreported functional genomic regions, some of which include SV and CNV regions. Previous studies on MRKH patients using low resolution microarrays or gene panels alone identified very few potentially causative gene variants [27]. Therefore, WGS may complement conventional array and WES to capture the complete landscape of the genome in MRKH.

Multiple lines of evidences have reported that complex genetic factors are implicated in MRKH (CNAs, indels and SNVs) [2,10,[28], [29], [30], [31], [32], [33], [34]]. The clinical characterizations of MRKH syndrome are varied and complicated. As previously reported, the genes mutations linked with MRKH syndrome are believed to be associated with early embryonic developmental and might be located in several chromosomal loci [35]. Many individuals with MRKH have potential defects on several chromosomes, including 1, 4, 7, 8, 10, 11, 16, 17, 22 and X [33,[36], [37], [38], [39], [40]].

In these 9 trio families’ genomes, 15 SVs occurred in type of gene-gene and were predicted as a direct link to a gene. These results suggested that, different mutational mechanisms could act on the patients with MRKH syndrome, and the results indicated that MRKH syndrome have alterations affecting the short length of chromosome, such as chr1, chr13 and chr19. (Fig. 1 landscape of SV).

Among the significant MRKHS-associated genes identified in our study, BAZ2B, KLHL18, PIK3CD, SLC4A10 and TNK2 showed extensive genetic pleiotropy:

PIK3CD (OMIM: 602839; RefSeq mRNAs: NM_005026) generates phosphatidylinositol 3,4, 5-trisphosphate (PIP3). PIP3 plays a key role by recruiting PH domain-containing proteins to the membrane, including AKT1 and PDPK1, and activating signaling cascades involved in cell growth, survival, proliferation, motility and morphology.

The c.1706C > G mutation in TNK2 indicated that it should be correlated with Ack1 tyrosine kinase activation occurring on chromosome 3. Specifically, TNK2 was reported to be related to tumors, such as non-small cell lung cancer (PMID: 24461128) and bronchopulmonary dysplasia (PMID: 23897914). Diseases associated with this gene include infantile-onset mesial temporal lobe epilepsy with severe cognitive regression.

SLC4A10 is a membrane-bound transporter that facilitates is transfer of glucose and other sugars, bile salts, organic acids, metal ions and amine compounds into the cell. The role of these genes in the development of disease needs to be studied in model organism (e.g. mouse) (Table 1).

Table 1

Summary of the three damaged de novo SNVs by whole-genome sequencing (WGS) analysis.

No.	Genes	Cytoband	Genome_Change	cDNA_Change	Codon_Change	Protein_Change	Family ID
1	TNK2	3q29	g.chr3:195597011G > C	c.1517C > G	c.(1516-1518)cCc > cGc	p.P506R	trio04
2	PIK3CD	1p36.2	g.chr1:9780813G > A	c.1535 G > A	c.(1534-1536)cGg > cAg	p.R512Q	trio06
3	SLC4A10	2q24.2	g.chr2:162719515C > T	c.709C > T	c.(709-711)Cgt > Tgt	p.R237C	trio09

The de novo mutations (chr1: 145291369 G > A) was identified in the 3'UTR of NOTCH2NL, indicating the important influence of NOTCH2NL on MRKH processing. This gene is typically associated with Hajdu-Cheney syndrome and the region containing these genes involves the copy number alterations and chromosomal structural variation based on the literature.

However, this compound inheritance model gene dosage model is less likely to underlie the etiology of MRKH, Therefore, trait-specific genetic modifiers of the identified MRKHS-associated genes, second locus genomic modifiers in enhancers/promoters, and potential gene-environment interactions still demand further investigation.

In summary, the genetic diversity identified within our population further highlights the clinical heterogeneity of MRKH. BAZ2B, KLHL18, PIK3CD, SLC4A10 and TNK2 were identified as novel disease-associated genes. The ethnicity upon our rare-inheritance model of disease is unknown. The functions and mechanisms of all deleterious mutations in MRKH candidate genes of the trio families warrant further study.

Conclusion

We identified five de novo mutations in BAZ2B, KLHL18, PIK3CD, SLC4A10 and TNK2 by performing WGS, the functional involvement of all deleterious mutations in MRKH candidate genes of the trios warrant further study. WGS may complement conventional array to capture the complete landscape of the genome in MRKH.

Declaration of Competing Interest

A signed declaration of no conflict of interest has been provided by all authors involved in the study.