VIM proteins regulate transcription exclusively through the MET1 cytosine methylation pathway

Abstract

In Arabidopsis, variant in methylation (VIM) proteins are required for the maintenance of DNA methylation in the CpG dinucleotide context. VIM1 acts as a cofactor of DNA methyltransferase 1 (MET1), although the mechanism for this co-regulation remains unclear. In this study, we used RNA-seq analysis to profile the transcriptomes of vim1, vim1 vim2 vim3, and met1 null mutants. Consistent with previous studies indicating functional redundancy between these VIM proteins, we found no transcripts that were significantly misregulated in vim1 mutants. However, we identified a large set of VIM protein regulatory targets through analysis of vim1 vim2 vim3 mutants, and we observed that this set is essentially identical to that regulated by MET1. Log₂ fold changes in gene expression relative to wild type are strongly correlated between vim1 vim2 vim3 and met1 mutants. While the largest subset of these transcripts is upregulated and enriched with transposable elements, we also found small subsets of downregulated genes in each mutant, which are enriched with protein-coding genes. Together, these results expand on previous studies that profiled cytosine methylation in the vim1 vim2 vim3 mutant, and show that VIM proteins function in transcriptional regulation via their roles in the MET1 DNA methylation pathway.

Keywords: :

• DNA methylation
• RNA-seq
• epigenetics
• gene expression
• transposable element

Introduction

Cytosine methylation is a mechanism used by many eukaryotes to mediate transcriptional repression. In many organisms, such as mammals, cytosine methylation exists predominantly in the CG dinucleotide context.^1,2 In plants such as Arabidopsis, methylation can additionally occur in the CHG or CHH contexts (where H represents A, C, or T), but CG dinucleotides are still the most abundant methylation sequence context.^3,4 Cytosine methylation in Arabidopsis functions primarily to silence transposons and other repetitive elements,^5,6 although the bodies of expressed genes can also be methylation targets.^7,8 Methylation is established through the activity of the Domains Rearranged Methyltransferase 2 (DRM2), an ortholog of mammalian DNMT3 methyltransferases that can act on all three sequence contexts.⁹ Maintenance of CG and CHG methylation are performed by DNA methyltransferase 1 (MET1) and chromomethylase 3 (CMT3), respectively.^10,11 MET1 is an ortholog of mammalian DNMT1, which is present at DNA replication forks and actively methylates CG sites on the newly synthesized daughter strands of DNA.^12,13

The maintenance of CG methylation in plants also requires the activity of additional proteins, such as the Variant In Methylation (VIM/ORTH) family.^14-16 In Arabidopsis, this family is composed of five highly similar proteins, some of which (VIM1, 2 and 3) are expressed in vegetative tissues, while VIM5 is more highly expressed in endosperm.^17,18 Its first identified member, VIM1, was found through a screen for trans-acting mutations that affect DNA methylation in natural Arabidopsis strains.¹⁴ Later, it was found that VIM1 functions redundantly along with VIM2 and VIM3 to maintain CG methylation at centromeres and other repetitive loci.¹⁷ VIM1 has been shown to bind DNA containing methylcytosine in vitro, with a particular specificity for hemimethylated cytosines.¹⁹ Because of this activity, it was proposed that VIM proteins function to recruit MET1 to CG sites through the activity of its methylcytosine binding SET- and RING-associated (SRA) domain, as MET1 lacks a methylcytosine binding module.

Studies on the VIM animal homolog UHRF1 have shown that it fills a very similar biological role,²⁰ with at least two potential mechanisms having been proposed for its function. One group suggested that a direct interaction occurs between the UHRF1 SRA domain and a region of DNMT1, supporting their model with a yeast two-hybrid assay and in vitro approaches.²¹ More recently, a separate set of results was published showing that UHRF1 ubiquitylates histone H3 in frog oocyte extracts through the activity of its RING domain.²² This study found that DNMT1 was recruited to target loci by binding the ubiquitylated histone H3. In both of the proposed mechanisms, UHRF1 and DNMT1 coordinately regulate cytosine methylation at target loci.

The recent emergence of whole-genome bisulfite sequencing (BS-seq) has allowed DNA methylation profiling in Arabidopsis across a variety of conditions.^3,4 Notably, methylomes were recently mapped in a large number of gene silencing mutants.²³ That study concluded that VIM1, VIM2, and VIM3 function entirely redundantly, and noted that the methylation profile of a vim1 vim2 vim3 mutant strongly resembled that of a met1 null mutant. However, these profiling results did not address the question of whether VIM proteins serve additional roles in transcriptional regulation outside of their functions in maintaining DNA methylation. Uncertainty remains over whether VIM proteins are involved in other chromatin modification pathways, raising the possibility that VIM proteins transcriptionally regulate additional loci outside the set of MET1 targets. Alternatively, VIM proteins might function exclusively as MET1 cofactors, in which case there would be a complete overlap in the sets of genes regulated by VIM proteins and MET1.

In this study, we investigate the transcriptional regulatory function of VIM proteins using RNA-seq. We sequenced whole transcriptomes of met1, vim1, and vim1 vim2 vim3 mutants to compare the sets of regulatory targets. Strikingly, the met1 and vim1 vim2 vim3 target sets are almost identical, suggesting that VIM proteins do not have additional transcriptional regulatory functions outside their role in maintaining DNA methylation. We also found no significant differentially expressed targets in the vim1 mutant, confirming studies that have shown functional redundancy between VIM1 and other VIM proteins.¹⁷ We have generated an extensive list of regulated targets, including transposable elements, protein coding genes, and noncoding RNA transcripts. Comparable to previous studies in CG methylation mutants,^4,7 we found that most of the regulated targets are transposable elements. Together, these analyses provide a comprehensive picture of the targets regulated by CG methylation in Arabidopsis.

Results

VIM protein and MET1 target gene sets overlap extensively

We extracted RNA from vim1, vim1 vim2 vim3, met1, and wild type plants, then prepared strand-specific Illumina RNA-seq libraries in biological triplicates according to a protocol developed by Wang et al.²⁴ After final library cleanup and amplification, the average fragment size of each library ranged between 330 and 370 base pairs. Sequencing data was generated on the Illumina Hi-Seq 2000 platform, and the resulting reads were mapped to the Arabidopsis TAIR10 genome build using TopHat. Because transposable elements are the predominant targets of DNA methylation in Arabidopsis, we included the TAIR10 annotated transposable elements in our transcriptome reference file. To increase confidence in our results, we included only reads that mapped uniquely to a single genomic location in our differential expression analysis by filtering the TopHat output files for a MAPQ value of 10 or greater. Reads were aligned with no more than two mismatches. Summary statistics of the data are presented in Table 1.

Table 1. Summary statistics for the RNA-seq data

	Number of reads	Mapped reads		Uniquely mapped reads	% Uniquely
	(in millions)	(in millions)	% Mapped reads	(in millions)	mapped reads
vim1 #1	22.3	21.9	98.2	19.1	85.4
vim1 #2	33.7	33.0	97.9	27.9	82.9
vim1 #3	34.2	33.6	98.2	27.3	79.8
vim1 vim2 vim3 #1	28.5	27.8	97.6	23.1	81.3
vim1 vim2 vim3 #2	35.5	34.8	97.9	29.3	82.5
vim1 vim2 vim3 #3	33.2	32.4	97.4	26.7	80.3
met1 #1	70.6	69.0	97.8	57.7	83.6
met1 #2	34.8	33.9	97.5	28.4	81.7
met1 #3	40.2	39.5	98.2	31.7	78.8
wild type #1	30.7	30.0	97.8	23.6	76.8
wild type #2	37.3	36.4	97.8	28.0	75.1
wild type #3	37.0	36.4	98.4	32.1	86.7

We analyzed differential expression between mutant and wild type samples using the R package DESeq.²⁵ No transcripts were significantly differentially expressed in the vim1 mutant (Fig. 1A), reflecting the functional redundancy that our group and others have observed between VIM1, VIM2, and VIM3.^17,23 The analysis identified 2387 transcripts differentially expressed between vim1 vim2 vim3 mutant and wild type, and 2748 transcripts differentially expressed between met1 and wild type (Fig. 1B and C; Tables S1 and S2). Owing partially to advances in high-throughput sequencing technology and an additional focus on transposable elements in our analysis, the number of targets identified in the met1 mutant is substantially greater than that found by Lister et al. (2008).⁴ As expected, the distribution of p-values for each differential expression analysis is uniform, except for the low p-values associated with differentially expressed transcripts and the high p-values associated with transcripts having low read counts (Fig. S1).

Figure 1.Identification of differentially expressed transcripts in vim1 vim2 vim3 and met1 mutants. MA-plots of log₂ fold change transcript expression vs. mean normalized counts in vim1 (A), vim1 vim2 vim3 (B) and met1 (C) mutants relative to wild type. Red points represent differentially expressed transcripts meeting the adjusted P value cut-off of 0.05 or lower. Black points represent all other transcripts not called differentially expressed. Notably, no transcripts, other than VIM1, were identified as differentially expressed in the vim1 mutant, while over 2000 transcripts were differentially expressed in vim1 vim2 vim3 and met1 relative to wild type. The majority of the red points fall above the gray horizontal line at zero, indicating that most differentially expressed genes are upregulated in the vim1 vim2 vim3 and met1 mutants.

Remarkably, the sets of genes differentially expressed in vim1 vim2 vim3 and in met1 overlapped to a very large extent (Fig. 2A). To determine whether the magnitude of up- or downregulation was well correlated between the two data sets, we generated a scatter plot based on the expression profiles of the vim1 vim2 vim3 and met1 differentially expressed genes (Fig. 2B). Based on the paired log₂ fold change data for the transcripts shared between the two data sets, the Pearson’s correlation coefficient r = 0.984, indicating a very strong correlation exists between the expression profiles. When we examined the transcripts identified as misregulated in vim1 vim2 vim3 only or in met1 only, we also found a correlation between the two mutants in log₂ fold change relative to wild type (r = 0.938 and r = 0.920, respectively). Despite similar fold changes in expression for both mutants, many of these transcripts were called differentially expressed in only one mutant due to our adjusted p-value filtering.

Figure 2. VIM proteins and MET1 regulate extensively overlapping sets of transcripts. (A) Venn diagram illustrating the large overlap between the targets differentially expressed in met1 mutants (green circle) and the targets differentially expressed in vim1 vim2 vim3 mutants (orange circle). Numbers of targets specific to each mutant and shared between the two mutants are indicated. (B) Scatter plot showing high correlation between the log₂ fold change expression values in met1 relative to wild type (y-axis) and vim1 vim2 vim3 relative to wild type (x-axis). Green data points represent transcripts differentially expressed only in met1; orange data points correspond to transcripts differentially expressed only in vim1 vim2 vim3; and gray data points represent transcripts differentially expressed in both mutants relative to wild type. As a confirmation of the vim1 vim2 vim3 mutant genotype, the VIM1, VIM2, and VIM3 transcripts are represented with purple diamonds labeled 1, 2, and 3, respectively. Transcripts with an expression value of zero in the mutant or wild type sample were excluded from this analysis. A black line representing y = x is shown for reference.

Small subsets of differentially regulated transcripts were identified as specific to the vim1 vim2 vim3 or met1 mutant plants. We chose six targets each from the vim1 vim2 vim3-specific set and the met1-specific set for follow-up expression analysis. However, none of these targets could be confirmed independently through reverse transcription quantitative PCR (RT-qPCR) (Fig. S2). On the other hand, several of the transcripts found to be differentially expressed in both vim1 vim2 vim3 mutants and met1 mutants relative to wild type were validated through further expression analysis (Fig. 3). One possible explanation for this discrepancy is that the identified vim1 vim2 and vim3-specific and met1-specific transcripts are primarily false-positive results. In accordance with this interpretation, many of the false discovery rate (FDR) adjusted p-values of the transcripts in these subsets are very near the cutoff of 0.05 that we established for significance, making them marginal differential expression calls. We conclude that the VIM proteins and MET1 regulate an essentially identical set of targets.

Figure 3. Validation of selected transcripts regulated by VIM proteins and MET1. RT-qPCR relative expression analysis of four upregulated targets and two downregulated targets from the set of transcripts co-regulated by VIM proteins and MET1 is shown. Similar fold changes in expression are seen in the vim1 vim2 vim3 mutant (gray bars) and met1 mutant (white bars) relative to wild type (black bars). Expression values represent the averages calculated from three independent biological replicates. Error bars represent the standard error of the mean. For comparison of RT-qPCR with the RNA-seq results, fold changes relative to wild type from the RNA-seq data set are indicated above each bar.

Transcriptional regulation of transposable elements and protein-coding genes by VIM proteins

Consistent with the established role of CG methylation in transcriptional repression, most of the transcripts identified as differentially expressed relative to wild type in the vim1 vim2 vim3 and met1 mutants were upregulated. We observed that a great majority of the upregulated transcripts were transposable elements (Fig. 4), with a smaller number of protein-coding genes and other transcripts. This result matches previous studies that have established CG methylation as a mechanism for transposon silencing.⁵

Figure 4. Comparison of upregulated and downregulated targets differentially expressed in vim1 vim2 vim3 and met1 mutants. Pie charts show the proportions of transposable elements (black), protein-coding genes (white) and other transcripts such as microRNAs or unclassified RNAs (gray) among the genes upregulated (left) and downregulated (right) in both vim1 vim2 vim3 and met1 mutants. Transposable elements are highly represented in the upregulated set, while most of the downregulated genes are protein-coding. Numbers of targets in each category are indicated on the pie charts, with total numbers of targets in the upregulated and downregulated sets indicated in parentheses.

In each mutant, a smaller set of transcripts was downregulated relative to wild type, with protein-coding genes highly represented in this category (Fig. 4). Although protein-coding genes are enriched in the sets of downregulated genes, it is notable that the numbers of upregulated protein-coding genes is still at least double the number of downregulated protein-coding genes in each mutant. This bias may reflect that there are multiple possible modes of misregulation for protein-coding genes in the mutants. It is well established that protein-coding genes can be transcriptionally repressed when adjacent to a methylated transposable element.^26-28 One well-known target of this regulatory mechanism is the FWA locus, which has a promoter region comprised of repeats derived from a SINE retrotransposon.²⁶ Protein-coding genes might also be up- or downregulated in vim1 vim2 vim3 and met1 mutants due to secondary regulatory effects from de-repressed transcription factors.

Discussion

The functional interface between VIM proteins and MET1 is a subject of ongoing research. Until now, it has been unclear whether VIM proteins have transcriptional regulatory functions outside their roles in the maintenance of CG methylation. Here, we show that VIM proteins and MET1 regulate very similar sets of transcriptional targets, strongly indicating that the influence of VIM proteins on gene expression is exclusively via their function as MET1 cofactors. This conclusion expands upon previous studies that showed highly similar DNA methylation profiles between vim1 vim2 vim3 mutants and met1 mutants.^23,29

Beyond showing that the targets regulated by VIM proteins and MET1 are essentially the same, we demonstrated a high correlation between gene expression levels in vim1 vim2 vim3 and met1 mutants. This further reflects the mechanistic interplay between VIM proteins and MET1 in regulating CG methylation. Although the details of this process are unconfirmed in Arabidopsis, research on animal homologs suggests that VIM/UHRF proteins either recruit MET1 to target loci through direct interactions,²¹ or by depositing an ubiquitylated histone mark that serves as a binding substrate for MET1.²² Our results showing high correlation between the vim1 vim2 vim3 and met1 expression profiles are consistent with either model. Further research is necessary to elucidate how much similarity exists between the plant and animal systems for maintaining CG methylation.

Our previous research indicated that VIM proteins have a role in the regulation of chromatin structure aside from maintenance of DNA methylation. A comparison of centromeric structure between vim1 mutants and wild type plants showed that centromere repeats were relatively decondensed in the vim1 mutant. Additionally, the centromeric histone H3 variant HTR12 was mislocalized relative to the centromeric repeats in vim1.¹⁴ This subnuclear phenotype contrasted with that of a ddm1 mutant also deficient in centromeric CG methylation, which more closely resembled the wild type centromere organization. This observation indicated that VIM1, in particular, might have a role in the maintenance of chromatin structure unconnected from its function in DNA methylation maintenance. Although the results of our transcriptional profiling analysis do not exclude this possibility, it does suggest that any additional role VIM1 might have in chromatin modulation does not impact gene expression. We propose that VIM1 has a centromere-specific function in the maintenance of chromatin structure that does not affect gene or transposon mRNA expression. VIM1 may interact with centromere-specific epigenetic marks or proteins to produce cytosine methylation-independent chromatin modifications at centromeres.

Our comparison of gene expression between vim1 mutant and wild type plants showed that no transcripts were significantly differentially expressed in the mutant. This result is consistent with previous studies that show high functional redundancy between VIM1, VIM2, and VIM3 in vegetative tissues.^17,23 However, we demonstrated that the vim1 mutation was sufficient to produce a hypomethylation phenotype at centromeric repeats, while vim2 and vim3 mutants did not have a significant loss of methylation at centromeres.^14,17 It is surprising that vim1 centromeres are hypomethylated, yet no effect on mRNA transcription was detected in vim1 plants. It is possible that transcriptional changes are occurring within the vim1 centromere, as other epigenetic silencing mutants derepress transcription of 180-bp centromere repeats,³⁰ but our enrichment for poly-adenylated transcripts does not allow us to effectively assess the expression of these repeats. Further research designed to detect such transcripts could be undertaken in vim1 samples to determine whether VIM1 is required to silence centromere repeat transcription.

Materials and Methods

Plant materials

The met1-3 and vim1-2 vim2-2 vim3-1 mutants used in this study have been previously described.^17,31 Since met1-3 homozygous mutants are sterile, we obtained met1-3 homozygotes from the segregating first-generation progeny of a met1–3 heterozygote provided by Jerzy Paszkowski (University of Geneva). The homozygotes were confirmed through genotyping PCR. Plants were grown under a 16 h photoperiod at 22 °C in controlled growth chambers.

RNA-seq library preparation

Floral bud tissue was collected from wild type, vim1 vim2 vim3, and met1 plants in three independent biological replicates, and total RNA was isolated using a Qiagen RNEasy Plant Miniprep Kit. Illumina RNA-seq libraries were prepared according to a strand-specific protocol.²⁴ Briefly, poly-A mRNA was selected using oligo dT Dynabeads (Life Technologies 61002) and fragmented in SuperScript II (Life Technologies 18064-014) First Strand Synthesis Buffer at a final concentration of 5 mM MgCl₂ for 10 min at 95 °C. Double-stranded cDNA was then synthesized from the mRNA, with the second strand labeled with dUTP. After end repair and adenylation, adapters were ligated to the cDNA, and libraries were size selected using a double-sided SPRI method with Ampure XP beads (Beckman Coulter A63880). The libraries were PCR amplified in 15 cycles using barcoded primers. Samples were multiplexed and sequenced on an Illumina Hi-Seq 2000 machine by the Genomics Resources Core Facility at Weill Medical College. Raw read data has been deposited in the NCBI Sequence Read Archive (SRX471361).

Sequence read alignment and differential expression analysis

Reads were aligned to the TAIR10 assembly of the Arabidopsis thaliana genome using TopHat version 2.0.4.³² Since transposons are known targets of regulation by CG methylation, we included the TAIR10 transposable element annotation based on Quesneville et al.³³ in our GTF reference file. Since there is some redundancy between this annotation and the AGI-annotated transcriptome, we filtered out all AGI-annotated transcripts that overlapped with the transposable element annotation. The output bam files were sorted with SAMtools,³⁴ and differential expression analysis was performed using the R package DESeq.²⁵ To define differentially expressed genes, an adjusted P value cutoff of 0.05 and FDR threshold of 0.05 were used. No fold-change cutoff for expression was applied.

Reverse transcription quantitative PCR

Total RNA was extracted from inflorescence buds and used as a template for cDNA synthesis with SuperScript III reverse transcriptase (Life Technologies 18080-051). Quantitative PCR analysis was performed on a Bio-Rad iCycler machine with a home-made SYBR Green master mix. Actin2 (At3g18780) was used as a constitutively expressed transcript for normalization. The reported results are averages of at least three independent biological replicates.

Disclosure of Potential Conflicts of Interest

No potential conflicts of interest were disclosed.

Acknowledgments

We acknowledge Natalie Henkhaus, Erika Hughes, and Haiyi Wang for contributing helpful insights on the results. We thank Jerzy Paszkowski (University of Geneva) for providing the met1-3 line. This research was funded by the National Institutes of Health (GM078256 to E.J.R.), the Department of Defense (BC093067 to M.S.S.), and institutional funds from the Boyce Thompson Institute.

10.4161/epi.28906