Cryptorchidism, defined as the failure of one or both testes to descend into the scrotum, represents one of the most common human congenital anomalies. While it can be part of hundreds of congenital malformation syndromes, it mainly appears as an isolated condition (Barthold et al., Reference Barthold, Reinhardt and Thorup2016). The major clinical consequences of this abnormality include reduced fertility in adulthood and an increased risk for testicular cancer (Hutson et al., Reference Hutson, Balic, Nation and Southwell2010).
The etiology of non-syndromic cryptorchidism remains unknown in most cases and seems to be multifactorial. Familial clustering and increased prevalence in first-degree relatives suggests that genetics may play a role in the development of this disease (Elert et al., Reference Elert, Jahn, Heidenreich and Hofmann2003; Schnack et al., Reference Schnack, Zdravkovic, Myrup, Westergaard, Wohlfahrt and Melbye2008). Mutations in the gene encoding insulin-like factor 3, which regulates the transabdominal phase of testicular descent, and in its receptor, RXFP2, have been linked to approximately 5% of cryptorchidism cases (Foresta et al., Reference Foresta, Zuccarello, Garolla and Ferlin2008). However, analysis of the HOXA10 gene failed to show consistent results (Cheng et al., Reference Cheng, Wang, Xu, Pei, Liu, Huang and Lu2017; Kolon et al., Reference Kolon, Wiener, Lewitton, Roth, Gonzales and Lamb1999). Similarly, candidate gene studies showed conflicting associations of the androgen receptor (AR) and estrogen receptor 1 (ESR1) genes with cryptorchidism (Barthold et al., Reference Barthold, Reinhardt and Thorup2016). Recent genome-wide association studies (GWASs) have identified several loci in suggestive association with cryptorchidism, although none of these have shown genome-wide significance (Barthold et al., Reference Barthold, Wang, Kolon, Kollin, Nordenskjöld, Olivant Fisher and Devoto2015). The results of the previous studies are thus suggestive of the complex genetic architecture that is associated with cryptorchidism.
Maternal and environmental factors are also believed to contribute to the etiology of cryptorchidism. Maternal smoking and gestational diabetes increase the risk for this disease (Zhang et al., Reference Zhang, Wang, Zheng, Liu, Zhang, Zheng and Chen2015). Furthermore, pieces of evidence from animal models and human epidemiological studies suggest that exposure to endocrine-disrupting chemicals (EDCs), such as isphenol A, phthalates, and pesticides, may be causally associated with cryptorchidism (Virtanen & Adamsson, Reference Virtanen and Adamsson2012). EDCs act as antiandrogens or estrogens that disrupt normal hormonal stimulation during testicular descent. Because these exogenous factors have the potential to elicit additional epigenetic effects, it is highly likely that epigenetic alterations are also involved in the etiology of cryptorchidism (Guerrero-Bosagna & Skinner, Reference Guerrero-Bosagna and Skinner2014).
Monozygotic (MZ) twins offer the ideal model to dissect the relative contributions of genetic and environmental factors in the pathogenesis of multifactorial diseases (Baranzini et al., Reference Baranzini, Mudge, van Velkinburgh, Khankhanian, Khrebtukova, Miller and Kingsmore2010; Jin et al., Reference Jin, Zhu, Hu, Liu, Li, Li and Chen2014). It is generally considered that MZ twins are identical in genotype, and phenotypic discordance is ascribed to environmental factors. Previous studies have detected genetic and epigenetic differences between MZ twins (Bruder et al., Reference Bruder, Piotrowski, Gijsbers, Andersson, Erickson, Diaz de Ståhl and Dumanski2008; Yuan et al., Reference Yuan, Xia, Bell, Yet, Ferreira, Ward and Spector2014). Genetic differences may result from very early post-twinning mutational events, and epigenetic changes are driven by environmental factors (Heijmans et al., Reference Heijmans, Tobi, Stein, Putter, Blauw, Susser and Lumey2008; Kondo et al., Reference Kondo, Schutte, Richardson, Bjork, Knight, Watanabe and Murray2002). In the present study, genomic and epigenomic differences were analyzed in a pair of MZ twins discordant for cryptorchidism for the first time, aiming to identify the specific alterations responsible for cryptorchidism.
Materials and Methods
Case Report
The patient was an 11-month-old boy who had undergone orchidopexy at the Children's Hospital of Chongqing Medical University because of bilateral undescended testes, which were localized within the inguinal canal. He had a healthy MZ twin brother who showed normal testicular descent. The twins were delivered by the Cesarean section after 39 weeks of gestation to healthy, non-consanguineous parents. The patient had no family history and no other clinical symptoms except for cryptorchidism. A written informed consent was obtained from the parents of the twins. The study protocol was in accordance with the tenets of the Declaration of Helsinki and approved by the institutional ethnics committee of the Children's Hospital of Chongqing Medical University.
Genomic DNA Extraction
Approximately, 5 mL of peripheral blood was obtained from both twins using a heparin anticoagulation tube. Genomic DNA was extracted from leukocytes using a DNA isolation kit (Tiangen, Beijing, China) according to the protocol provided by the manufacturer. DNA samples were quantified using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific Inc., Waltham, MA). Purified genomic DNA was stored at −80 °C until use.
Zygosity Analysis
To confirm that the twins were MZ, amelogenin and 21 autosomal short tandem repeat (STR) profiles were obtained using AGCU EX22 STR Kit (AGCU ScienTech Inc., Wuxi, China) according to the manufacturer's protocols. All of the amplification reactions were performed in a GeneAmp PCR system 9700 (Perkin Elmer, Norwalk, CT). The amplified products were run on an ABI 3130 Genetic Analyzer (Applied Biosystems, Foster City, CA), and the data were analyzed using GeneMapper ID v3.2 software (ABI).
Whole Exome Sequencing (WES) and Data Analysis
Genomic DNA was randomly broken into manageable fragments and subjected to library preparation. Paired-end libraries with a span size of ~300 bp were adopted for exome sequencing. Exome enrichment was performed using the SureSelect Human All Exon 50 Mb kit (Agilent Technologies, Santa Clara, CA). The enriched libraries were sequenced using the Illumina HiSeq Xten sequencing system (Illumina, San Diego, CA), and 150-bp paired-end reads were generated. The sequencing data were mapped to the hg19 human reference sequence using the Burrows–Wheeler Aligner (BWA) program (Li & Durbin, Reference Li and Durbin2010). PCR duplicates were marked using the Picard program. Genome Analysis Toolkit (GATK) was then used for indel realignment and base quality recalibration (McKenna et al., Reference McKenna, Hanna, Banks, Sivachenko, Cibulskis, Kernytsky and DePristo2010). Single nucleotide variants (SNVs) and small insertions and deletions (InDels) were called and recalibrated using GATK. Variants were annotated using the ANNOVAR (Wang et al., Reference Wang, Li and Hakonarson2010).
To help identify candidate mutations, the following variants were filtered out: (1) synonymous substitutions, intergenic and intronic variants, and variants located within untranslated regions; (2) variants present in dbSNP142 and 1,000 Genomes databases; (3) variants with a minor allele frequency of >1%; and (4) variants with the phred-scaled quality score of <30 or a read depth of <4. In addition, only variants specific to the affected twin were further analyzed.
Validation of Discordant Variation
A total of 37 discordant variations identified by Whole Exome Sequencing (WES) were selected for amplification from genomic DNA by polymerase chain reaction (PCR) using specific primer pairs (Table S1). The resulting PCR products were purified and sequenced on an ABI 3730XL sequencer (ABI) using a BigDye Terminator Cycle Sequencing Kit (ABI).
Methylated DNA Immunoprecipitation Sequencing (MeDIP-Seq) and Bioinformatic Analyses
Genomic DNA samples were fragmented at sizes ranging from 200 bp to 900 bp with a Bioruptor sonicator (Diagenode, Liège, Belgium). Approximately, 800 ng of sonicated DNA was end-repaired, attached with polyA tails, and ligated to single-end adapters following the standard Illumina genomic DNA protocol. After agarose size-selection to remove unligated adapters, the adaptor-ligated DNA was used for immunoprecipitation using a monoclonal antibody against 5-methylcytidine (Diagenode). MeDIP enrichment was performed according to the previously described protocol (Tao et al., Reference Tao, Fang, Hu, Lu, Cao, Zhao and Pan2015). About 150 ng of the affected twin's sonicated and adaptor-ligated DNA that did not undergo immunoprecipitation was used as input. The input and antibody-enriched DNA was amplified by PCR using the single-end Illumina PCR primers and then purified by agarose gel electrophoresis. The products were then quantified on an Agilent 2100 Analyzer (Agilent) and subjected to qPCR qualification. Sequencing was performed on Illumina HiSeq 2000 using the TruSeq Rapid SBS Kit (Cat. # FC-402-4001, Illumina).
High-quality clean reads were aligned to the human reference genome (hg19) using SOAP2 (Version 2.21) with no more than 2 bp mismatches (Li et al., Reference Li, Yu, Li, Lam, Yiu, Kristiansen and Wang2009). Only uniquely mapped reads were further analyzed. The distribution of reads in chromosomes and in the different genome components was analyzed. Genome-wide methylation peak scanning was performed using MACS V1.4.2 (Zhang et al., Reference Zhang, Liu, Meyer, Eeckhoute, Johnson, Bernstein and Liu2008). The region with a p value of <10−5 was defined as a peak. The distribution of peaks in different genome components was also analyzed. The genomic differentially methylated regions (DMRs) between the twins were detected using a previously described method (Kim et al., Reference Kim, Kim, Seo and Kim2016), except that the significance cutoff of p value for the chi-square statistics and false discovery rate (FDR) was 0.01. The genes with at least one region overlapping with DMRs were considered differentially methylated genes. All differentially methylated genes were then subjected to Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis using the DAVID bioinformatics resources (Huang et al., Reference Huang, Sherman and Lempicki2009).
Results
Zygosity Analysis
STR analysis revealed that the twins had the identical alleles for all 21 microsatellite markers (Table 1), thereby confirming that they were indeed MZ.
TABLE 1 Short Tandem Repeat Typing Results of the Patient and Twin Brother
WES Analysis and Variant Validation
WES was performed on genomic DNA of the twins to identify the underlying genetic mutations. An average of 12.44 Gb of effective sequence data was obtained per individual. The mean sequencing depth of the target region of each individual was 160× and 161×, respectively. Furthermore, >98% of the target region in both samples showed at least 20× coverage (Table S2).
Bioinformatics analysis identified a total of 97,692 SNVs and 15,678 InDels in the affected twin, whereas 99,249 SNVs and 15,976 InDels were detected in the healthy twin (Tables S3 and S4). Given the discordant phenotypes of the twins, 1,368 SNVs and 917 InDels specific to the affected twin were identified. Non-synonymous variations, splice site variations, and InDels were selected for further analysis. After applying the stringent filtering criteria, approximately 35 discordant SNVs and 11 InDels were obtained (Figure 1; Table S5).
FIGURE 1 Illustration of the variant filtering procedure. Note: SNV = single nucleotide variation, InDels = small insertions and deletions; MAF = minor allele frequency.
All 35 discordant SNVs were selected to validate by the Sanger sequencing. However, none of these variations could be confirmed. While 34 of these SNVs were found to be false-positive calls, an SNV was turned out to exist in both twins. Two InDels in genes that were expressed in fetal testis according to the GeneCards database (http://www.genecards.org/) were selected for validation, neither of which was validated. And these InDels were proved to be false-positive calls.
Genome-Wide DNA Methylation Analysis
MeDIP-Seq was performed to detect epigenomic differences between the twins. A total of 33,033,859 clean reads were generated per sample, 69.9% of which could be uniquely aligned to the human genome (Table S6). In addition, the reads were detected in most chromosomal regions (Figure S1). The analysis of the read distribution in different genome components revealed that these were mostly located in intronic regions, which contained 39% of the unique reads (Figure S2A). Within genic regions, the number of reads dramatically decreased within the transcription start sites of both samples, whereas a gradual increase was observed in intragenic regions (Figure S2B).
Whole genome methylated peak scanning revealed approximately 51,647 peaks in the affected twin and 11,154 peaks in the healthy twin, covering 2.14% and 0.36% of the genome, respectively (Table S7). The majority of the peaks were present in intronic regions, followed by coding sequences (CDSs; Figure S3). The differentially methylated genes between the twins were then determined. A total of 5,410 genes were differentially hypomethylated, and 2,383 genes were differentially hypermethylated in the affected twin compared to the healthy twin (Supplementary file 2). Most of these genes were differentially methylated in the intronic regions, followed by the CDSs (Figure S4).
GO analysis was performed to determine biologically functional relationships among differentially methylated genes. In addition, GO terms enriched in the biological process were identified. The results showed that eight hypermethylated genes were related to microtubule-based movement (FDR < 0.05; Supplementary file 3). The hypomethylated genes were mainly correlated with 36 categories and 1,056 genes were involved (FDR < 0.05; Supplementary file 4). The top 15 categories included protein phosphorylation, cell migration, positive regulation of GTPase activity, cell adhesion, and protein autophosphorylation (Figure 2A).
FIGURE 2 Gene ontology and pathway analysis of hypomethylated genes. A. Top 15 biological processes involving hypomethylated genes. B. Top 15 significant pathways involving hypomethylated genes.
KEGG pathway analysis revealed the significant pathways that were associated with differentially methylated genes. The hypermethylated genes were not observed to tightly associate with any pathways. The hypomethylated genes were significantly enriched in 25 pathways and 408 genes were involved (FDR < 0.05; Supplementary file 5). The top 15 significant pathways included calcium signaling, Rap1 signaling, axon guidance, cAMP signaling, and glutamatergic synapse (Figure 2B).
Discussion
In the present study, no genetic differences were detected in MZ twins discordant for cryptorchidism. The results coincided with the findings of earlier similar studies, in which no genetic discordances were found in MZ twins discordant for several disorders, including renal agenesis (Jin et al., Reference Jin, Zhu, Hu, Liu, Li, Li and Chen2014), cataract (Wei et al., Reference Wei, Sun, Hu, Yang, Qiao and Yan2015), VACTERL association (Solomon et al., Reference Solomon, Pineda-Alvarez, Hadley, Hansen, Kamat and Donovan2013), hypothyroidism (Magne et al., Reference Magne, Serpa, Van Vliet, Samuels and Deladoëy2015), heart defects, and diaphragmatic hernia (Zhang et al., Reference Zhang, Thiele, Bartmann, Hilger, Berg, Herberg and Reutter2016). Our findings provide additional evidence that genetic mutations are not a common cause for phenotypic discordance among MZ twins.
MeDIP-seq is a cost-effective method in establishing the global DNA methylation patterns of the human genome. Using this methylome technology, epigenetic changes were detected in the current twin pair, thereby suggesting that the discordant phenotype might have resulted from alterations in DNA methylation patterns. Although MZ twins share a single uterus, they might be exposed to different intrauterine environments such as uneven blood supply and different nutritional conditions (Czyz et al., Reference Czyz, Morahan, Ebers and Ramagopalan2012). The epigenetic changes triggered by non-shared early exposures might thus affect the expression of regulatory genes during embryogenesis, thereby leading to a particular disease phenotype.
The epigenetic alterations underlying the etiology of cryptorchidism are poorly understood to date. Our comparison of DNA methylation profiles between the twins revealed that differentially methylated genes were enriched with several significant biological processes and pathways. Among these, the actin cytoskeleton pathways, which were identified by both GO and KEGG pathway analysis (Figure 2), require further investigation. In previous GWAS on non-syndromic cryptorchidism, pathway analysis of suggestive signals also identified various genetic loci that were related to cytoskeleton-dependent functions (Barthold et al., Reference Barthold, Wang, Kolon, Kollin, Nordenskjöld, Olivant Fisher and Devoto2015). Moreover, alterations in the expression of cytoskeleton-related genes were observed in a rat strain with congenital cryptorchidism (Barthold et al., Reference Barthold, McCahan, Singh, Knudsen, Si, Campion and Akins2008). Therefore, despite the need to confirm our methylation results in additional samples, the present study offers the clues that aberrations in the methylation patterns of cytoskeleton-associated genes may contribute to the pathogenesis of cryptorchidism.
To the best of our knowledge, this is the first genome-wide study to investigate the genomic variations and epigenetic differences in discordant MZ twins with cryptorchidism. Previous studies have detected gene mutations in only a fraction of cryptorchidism cases. It has also been suggested that the pathogenesis of cryptorchidism is multifactorial, which include genetic and environmental factors. In the present study, no genomic mutations but epigenetic differences were detected in discordant MZ twins. Our findings provide evidence for the multifactorial nature of cryptorchidism, particularly non-genetic factors, thereby providing new insights into the pathogenesis of cryptorchidism.
The present study has two limitations. First, only one pair of MZ twins was investigated. More pairs of discordant MZ twins should be recruited in the future to confirm the findings of this study. Second, DNA methylation is tissue specific. The methylation patterns in leukocytes were examined in the present study, whereas cryptorchidism is a developmental defect that involves the gubernaculum and testis, which are generally difficult to obtain from human subjects. As a result, animal models should be used to verify the epigenetic alterations identified in this study.
In summary, we have performed WES and MeDIP-Seq analyses in a pair of MZ twins discordant for cryptorchidism. No genomic mutations were detected, whereas epigenetic alterations were observed, thereby suggesting that non-genetic factors may contribute to the pathogenesis of cryptorchidism.
Disclosure of Interest
None.
Supplementary material
To view supplementary material for this article, please visit https://doi.org/10.1017/thg.2017.33
Acknowledgments
This work was supported by the National Natural Science Foundation of China (grant number 81300463) and the Natural Science Foundation Project of Chongqing (grant number cstc2012jjjq10001).
