TRIM4 is associated with neural tube defects based on genome-wide DNA methylation analysis

Sample collection

DNA extraction and bisulfite conversion of genomic DNA

Genomic DNA was extracted from spinal cord tissues with 0.5% sodium dodecyl sulfate lysis buffer and protease K (1.5 mg/mL) for nuclear protein digestion at 56 °C for 1 h. Total genomic DNA was harvested using a QIAamp DNeasy Blood and Tissue Kit (QIAGEN, Hilden, Germany), followed by precipitation with 70% alcohol. For each sample, genomic DNA was quantified by NanoDrop ND-1000, and 500 ng genomic DNA was used for bisulfite conversion using an EpiTect Whole Bisulfitome Kit (Qiagen). Bisulfite-converted genomic DNA was eluted and stored at − 20 °C until use.

DNA methylation analysis

DNA from spinal cord tissues from the NTD group (three cases) and control group (three cases, matched by sex and gestational week) were used for the methylation array assay. Genomic DNA was extracted, purified, and bisulfite converted. Samples were randomized across 3MSA-4 plates for processing based on instructions for the Illumina Infinium HumanMethylation450 BeadChip [30]. To minimize batch effects, cases and controls were randomly distributed into different arrays. Raw intensity was read using Illumina GenomeStudio Software 2011.1, and background normalization was performed.

The raw data was first pretreated with R software minfi package. R software Illumina Methylation Analyzer (IMA) was then used to screen methylation levels and regions in the samples [31, 32]. Beta value (β value) was used as an indicator of methylation for each locus in each sample (β values range from 0 to 1, corresponding to completely unmethylated and fully methylated sites, respectively). Delta beta (Δβ) is defined as the difference in β values between the two groups in which the greater absolute value represents a greater degree of difference. DiffScore is the parameter and model measuring differences as provided by the Illumina Company in which a greater absolute value indicates a more significant difference. Delta beta and DiffScores were both used to differentially screen for methylated genes. Detection p values represent the confidence levels of chip signal values. Detection probes with a p value > 0.05 were not reliable and were therefore excluded from further analysis.

Differential methylation gene analysis

To identify differentially methylated genes between NTD cases and controls, the following five filters were applied: (i) absolute β value difference > 0.10, (ii) DiffScore > 13, (iii) detection p value < 0.05, (iv) participate and play an important role in multiple pathways, and (v) differentially methylated CpGs located in the transcriptional regulatory region of the gene, such as 5′ untranslated region (UTR), transcriptional start sites (TSS) 200, first exon, 3′UTR, and more.

To identify enriched gene functions associated with differentially methylated regions (DMRs), Gene Ontology (GO) analysis was performed using online gene enrichment tools from Shanghai Biotechnology Corporation (http://enrich.shbio.com/index/ga.asp). The Fisher accuracy test was used for enrichment analysis with the clusterProfiler package [33], which was released within the Bioconductor project. Genes associated with 450K CpG sites were ranked by unadjusted p values (from smallest to largest) and computed hypergeometric p values for overrepresentation of each biological process Gene Ontology (GO) category. A term containing 10 or more genes as well as a p value < 0.05 was considered a selection criterion for interested terms.

Real-time quantitative PCR

Total mRNA from spinal cord tissues from 14 cases and 28 controls (including 3 samples in microarray analysis) was extracted with TRIZOL (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instructions. Complementary DNA (cDNA) was synthesized using a PrimeScript™ RT reagent kit (TAKARA, Japan) with mRNA as the template. Transcript expression of candidate genes was evaluated using specific primers in a 20 μL reaction (2 μL of template cDNA, 1 μL of 10 μmol/L each primer, 10 μL of 2× SYBR Green Master Mix, 0.4 μL of ROXII, and 5.6 μL of ddH2O) on an ABI Prism 7500HT sequence detection system (Applied Biosystems) with the following cycling parameters: predenaturation at 95 °C for 30 s, 45 cycles of denaturation at 95 °C for 5 s, and annealing at 60 °C for 20 s. Relative levels of target gene expression were calculated according to the ΔΔCt method [34] and normalized to glyceraldehyde 3-phosphate dehydrogenase (GAPDH) gene expression. The primers for candidate genes and GAPDH are shown in Additional file 1: Table S1.

Pyrosequencing and sequence analysis

DNA from spinal cord tissues obtained from 14 cases and 14 controls (including 3 samples in microarray analysis) was used for pyrosequencing and sequence analysis. Pyrosequencing was performed using the PSQ Gold SQA reagent kit (Qiagen) on a Pyromark Q96 ID platform (Qiagen, Germany) according to the manufacturer’s instructions. Briefly, 50 μL biotinylated PCR products were mixed with 3 μL streptavidin sepharose beads (Amersham Biosciences AB, Sweden) and 47 μL binding buffer, followed by shaking at 2000 rpm for 10 min. The immobilized complex was captured by using a vacuum prep tool. To dissociate and discard the unbiotinylated strand, single strand purification was performed. The beads that bound to single biotinylated strands were released to a 96-well microtiter plate to which 49 μL annealing buffer, and 1 μL complementary sequencing primer (pBR-V1.AS) had been added. The processed mixture was loaded onto the PyroMark ID system for pyrosequencing set with 10 cycles of ATCG dispenses. The obtained sequences were directly used to search the database constructed. Results were identified as organisms with a 100% DNA sequence matching. Sequences for pyrosequencing primers are shown in Additional file 2: Table S2.

Next-generation sequencing and Sanger sequencing

Genomic structures of human TRIM4 genes were confirmed by NCBI GenBank (NM_033091.2 and NM_033091.3). Gene variants in TRIM4 were detected by next-generation sequencing. The Agilent SureSelect XT Custom enrichment system was used to construct genomic DNA-fragment libraries and target enrichment of TRIM4. Each fragmented genomic DNA library was connected to an index adapter, and the connected libraries were gel purified and PCR amplified (Phusion, Thermo Scientific). The Agilent 2100 biologic analyzer was used to determine the quantity and quality of the library. Forty-eight libraries were pooled in total and then hybridized to RNA library baits, after which the targeted sequences were purified and amplified (Herculase II fusion, Stratagene). Sequencing was performed on an Illumina HiSeq2000 DNA sequencer (version 3). Sequence alignment was performed with BWA software (Li and Durbin, 2009) based on the h19 database.

To confirm genotyping results for TRIM4 from next-generation sequencing, Sanger sequencing was designed to detect variants in all exons, 1 kb upstream of the transcription start site, and 3′UTRs in TRIM4. The primer design tool (https://www.ncbi.nlm.nih.gov/tools/primer-blast/) was used to design specific primers. Extracted DNA was amplified by polymerase chain reaction (PCR) using Pfu DNA Polymerase, MgCl₂ (25 mM), 10× PCR buffer (without Mg²⁺:100 mM Tris-HCl (pH 8.8, 25 °C), 500 mM KCl, 0.8% (v/v) Nonidet P40), and dNTP (10 mM) (Sangon Biotech, Shanghai, China). The optimal annealing temperature for PCR was 58 °C. Amplified PCR products were purified and recollected using a SanPrep Column DNA Gel Extraction Kit (Sangon Biotech, Shanghai, China). DNA sequencing was performed on an ABI 3730xl DNA sequencer (Applied Biosystem, Foster City, USA) following the manufacturer’s standard sequencing protocols. For data analysis, the dbSNP in NCBI, Genome 1000, and ExAC database was used as a reference dataset for rare variant allele frequency in a Chinese Han population.

Western blotting

Spinal cord tissues from 14 cases and 14 controls, including 3 samples in microarray analysis, were lysed in ice-cold RIPA buffer (Solarbio, R0010, China) with 1 mM phenylmethanesulfonyl fluoride (PMSF). After quantification using the BCA method, proteins were electrophoresed on 10% SDS-polyacrylamide gels and transferred to polyvinylidene difluoride (PVDF) membrane (Millipore, USA). The membranes were blocked with 5% non-fat milk in PBST containing 0.1% Tween-20 for 60 min and incubated with primary rabbit anti-RNF87 (TRIM4) (ab26300, 1:800, Abcam, Cambridge, UK) and a mouse anti-GAPDH antibody (sc-365062, 1:10000, Santa Cruz Biotechnology) overnight at 4 °C, followed by incubation with secondary polyclonal goat anti-rabbit HRP-conjugated (ZDR-5306, 1:2000, ZSGB-bio, Beijing, China) or polyclonal goat anti-mouse HRP-conjugated antibody (ZDR-5307, 1:2000, ZSGB-bio, Beijing, China) for 2 h at room temperature. Signals were visualized using enhanced chemiluminescence (ECL, Millipore, USA) reagents. Optical density values of each protein divided by the loading control (GAPDH) were regarded as the relative density of each protein.

MTHFR and MTRR polymorphism analysis

Polymorphism analysis was assessed at three loci in each case: at nucleotide 677C>T (rs1801133) and 1298A>C (rs1801131) in MTHFR and 66A>G (rs1801394) in methionine synthase reductase (MTRR). Human MTHFR and MTRR gene polymorphism detection kit (YZYMT-014, Wuhan, China) was used to determine the genotype of the sample by comparing changes in Ct values of wild-type and mutant primers by a real-time fluorescent TaqMan probe PCR method. PCR was performed on an ABI Prism 7500HT sequence detection system (Applied Biosystems) with the following cycling parameters: 42 °C for 5 min, 94 °C for 3 min, 45 cycles of denaturation at 94 °C for 15 s, and annealing at 60 °C for 60 s.

Statistical analysis

The genotype and alleles of MTHFR and MTRR were computed by direct counting. Alleles and genotypes distribution of MTHFR and MTRR were evaluated by Fisher’s exact test using the SPSS software package version 20.0. A non-parametric test was used to compare the differences in methylation and expression levels between NTDs and controls using GraphPad Prism 6 (GraphPad Software, California, USA). All data are represented as the mean ± standard error of the mean (SEM). p value less than 0.05 was considered statistically significant.

Screening for differentially methylated genes

To investigate whether DNA methylation was different between controls and NTD cases, we measured global DNA methylation in spinal cord tissues using the Illumina Infinium Human Methylation450 BeadChip, which examines over 485,000 CpG sites spanning the whole human genome.

Filtered by absolute Δβ difference > 0.10 and DiffScore > 13, a total of 4648 differentially methylated sites (DMSs) were screened out, of which 1987 sites were hypermethylated and 2661 were hypomethylated. After filtering out unreliable chip signal values (detection p value > 0.05), 461 DMSs were identified, including 330 hypomethylated sites and 131 hypermethylated sites. These DMSs were distributed in 243 genes, including 156 hypomethylated genes, 83 hypermethylated genes, and 4 genes with both hypermethylated and hypomethylated sites (ACTR3C, GALNT9, HLA-DQB1, HLA-DQB2). There were more hypomethylated sites in HLA-DQB1 and HLA-DQB2 genes than hypermethylated sites. There were more hypermethylated sites in the GALNT9 gene than hypomethylated sites. The number of hypomethylated and hypermethylated sites in the ACTR3C gene was the same. CpG distribution of DMSs was analyzed and compared to CpG distribution in the Illumina 450K array, which covers 99% of annotated RefSeq genes and exhibits a wide distribution of probes among CpG islands, shores (2 kb flanking the islands), shelves (2 kb flanking the shores), and sea (regions outside the previous three categories). The results showed that both hypomethylated and hypermethylated DMSs in NTD spinal cord tissues were preferentially situated in CpG islands rather than shores and shelves (Fig. 1). The distribution of DMSs in the gene region was also analyzed and compared to RefSeq genes. Interestingly, the global distribution of DMSs relative to RefSeq genes showed a significant enrichment in hypomethylated CpGs in the intergenic regions. In mainly translating regions, including TSS200-TSS1500, 5′UTR, exon 1, and 3′UTR, hypomethylated and hypermethylated DMSs were preferentially distributed in the TSS200-TSS1500 region (Fig. 2).

Bioinformatics analysis of candidate genes of interest

We performed Gene Ontology analysis of 243 methylated differential genes, of which only 131 hypomethylated and 78 hypermethylated genes were annotated in GO databases. After sorting the 131 hypomethylated genes according to the p value, we found 7 terms with more than 10 gene assemblies in biological process (BP) GO terms. They were GO:0001816, cytokine production (7.63%); GO:0045087, innate immune response (11.45%); GO:0002520, immune system development (8.40%); GO:0050776, regulation of immune response (9.92%); GO:0006955, immune response (14.5%); GO:0006952, defense response (15.27%); and GO:0010033, response to organic substance (21.37%) (Additional file 3: Table S3). Hypermethylated genes were detected in categories associated with lipid, cellular lipid, and organic acid metabolic processes; establishment of localization in the cell; cellular localization; response to organic substances; protein localization; transport; and transmembrane transport (Additional file 4: Table S4). A further selection of target genes with a significant number of DMSs in a transcriptional regulation domain identified 5 immune specific hypomethylated genes (TLR1, TRIM4, MAP2K2, CALCOCO2, GNAS) and 5 hypermethylated genes (SMPD3, EGFR, HSPB7, MAGT1, SCT) involved in metabolic processes, localization, and transportation pathways.

Candidate gene expression level analysis

To reveal the relationship between methylation status and gene expression levels of these target genes, real-time quantitative PCR was performed to detect mRNA expression levels. Of the hypomethylated differentially methylated genes (DMGs), TRIM4 and TLR1 expression was higher in NTDs, while MAP2K2, CALCOCO2, and GNAS expression showed no significant difference (Fig. 3a). No significant differences in expression levels of hypermethylated DMGs were found in NTDs (Fig. 3b). When evaluating NTD subgroups, we found that TRIM4 and TLR1 expression was increased, but MAP2K2 expression was decreased in cases with SB compared to controls. There were no significant differences in TRIM4, TLR1, and MAP2K2 expression in cases with CHC (Fig. 3a). No significant differences in expression levels of hypermethylated DMGs were found in both cases with SBs and CHCs (Fig. 3b). To mitigate any possible bias, we added a separated control group (control group 2) matched to NTD cases by fetal sex and gestational age for real-time PCR analysis of TRIM4. In line with the results of control group1, the expression of TRIM4 was significantly higher in the NTD group compared to controls group 2. Since the results of the two control groups were consistent, we combined the two control groups into one group (Fig. 3a).

Western blot analysis was performed to further verify changes in TRIM4 protein expression levels in NTDs. In accordance with the results obtained from real-time PCR, we found that TRIM4 protein levels were significantly increased in NTDs compared with controls (Fig. 4).

Validation of methylation levels in differentially expressed genes by pyrosequencing

Pyrosequencing was used to identify methylation levels in different transcription regulatory regions in TRIM4, TLR1, and MAP2K2 from 14 NTD samples and matched controls. However, pyrosequencing was only conducted for TRIM4 and TLR1 due to the low primer score for MAP2K2. CpG site information for pyrosequencing is shown in Additional file 5: S1. The regions targeted for pyrosequencing contain differential methylation sites obtained from microarray analysis as well as surrounding CG sites. We analyzed seven CpG sites in exon 1 of TRIM4 and found that the degree of methylation in four sites (POS1, POS3, POS4, POS6) was significantly decreased in NTDs compared with controls (Fig. 5a). Total methylation of the seven CpG sites in exon 1 between NTDs and controls was also significantly different (NTDs: mean ± SEM = 4.057 ± 0.4712, N = 14; controls: mean ± SEM = 5.977 ± 0.7482, N = 13. p = 0.0419) (Fig. 5b). None of the four CpG sites in the promoter region showed differences in methylation levels between NTDs and controls (Fig. 5c). Total methylation levels in these two regions showed no significant difference (NTDs: mean ± SEM = 4.656 ± 0.3877, N = 14; controls: mean ± SEM = 5.971 ± 0.6334, N = 13, p = 0.0539) (Fig. 5d). These results indicate that the exon 1 region in TRIM4 shows significant hypomethylation. The sequence information of the regions targeted for pyrosequencing in TRIM4 was shown in Additional file 6: S2. Pyrosequencing results of the TLR1 gene showed no difference in methylation levels between NTDs and controls (Fig. 6).

TRIM4 mutation analysis

To investigate whether SNPs contribute to increased TRIM4 gene expression in NTD cases, we reanalyzed the previous genomic DNA sequencing results obtained from a different population of 100 NTD cases that was independent of the 14 NTD fetus samples. Forty rare variants (MAF < 0.01) were identified, including 4 upstream variants, 32 intron variants, 1 downstream variant, 1 3′UTR variant, 1 missense variant, and 1 synonymous variant. We also detected an exon region around the DMS in the 14 spinal cord tissues and found 3 3′UTR variants, 3 missense variants, 2 synonymous variants, 1 stop-gain variant, and 6 novel variants. Of these, 3 3′UTR variants (rs2572009, rs1048705, rs2572010), 1 missense variant (rs76665876) and 2 synonymous variants (rs2247761, rs2247762) were identified in both sequencing methods. Only the missense variant (rs76665876) is rare and was predicted to have a potential pathogenic possibility by using Polyphen and AvSIFT programs (Polyphen score, 0.980; SIFT score, 0.14). Therefore, we also analyzed the association between mRNA expression and the rare missense variant (rs76665876) in the 14 spinal cord samples with NTDs. The rs76665876 variant was found in 3 of the 14 NTD cases, and TRIM4 gene expression was not significantly different. There were 8 cases with hypomethylation without the rs76665876 variant in NTDs, and TRIM4 mRNA expression was significantly increased in these cases. Collectively, the mRNA expression, DNA methylation, and sequencing data suggest that genetic variants in TRIM4 genes only slightly contribute to the pathogenesis of human NTD. The main factor affecting TRIM4 mRNA levels in NTDs is the changes in DNA methylation.

MTHFR and MTRR polymorphism analysis

To investigate whether the variants in MTHFR (C677T and A1298C) and MTRR (A66G) were associated with methylation levels of TRIM4, we compared TRIM4 methylation levels between NTDs with and without homozygote mutants in MTHFR and/or MTRR. No difference in methylation levels was found between the two groups (NTDs with homozygote mutants in MTHFR and/or MTRR: mean ± SEM = 4.617 ± 0.5224, N = 6; NTDs without homozygote mutants in MTHFR and/or MTRR: mean ± SEM = 3.638 ± 0.7189, N = 8, p = 0.3230).