The MADS-box transcription factor PHERES1 controls imprinting in the endosperm by binding to domesticated transposons

Abstract

MADS-box transcription factors (TFs) are ubiquitous in eukaryotic organisms and play major roles during plant development. Nevertheless, their function in seed development remains largely unknown. Here, we show that the imprinted Arabidopsis thaliana MADS-box TF PHERES1 (PHE1) is a master regulator of paternally expressed imprinted genes, as well as of non-imprinted key regulators of endosperm development. PHE1 binding sites show distinct epigenetic modifications on maternal and paternal alleles, correlating with parental-specific transcriptional activity. Importantly, we show that the CArG-box-like DNA-binding motifs that are bound by PHE1 have been distributed by RC/Helitron transposable elements. Our data provide an example of the molecular domestication of these elements which, by distributing PHE1 binding sites throughout the genome, have facilitated the recruitment of crucial endosperm regulators into a single transcriptional network.

Keywords: A. thaliana; MADS-box transcription factors; PHERES1; endosperm; genetics; genomics; imprinting; plant biology; seed; transposable elements.

Figure 1.. Identification and expression profile of PHE1 target genes.

(a) Spatial distribution of PHE1 binding sites around transcription start sites (TSS). The dotted pink lines indicate the spatial interval used to define PHE1 target genes. (b–c) Expression of PHE1 (b), and its target genes (c), across different stages of seed development. Gene expression is represented as a Log2-fold change between expression in the endosperm at the stages indicated on the x axis vs. expression in the pre-globular stage. A k-means clustering analysis was performed to group PHE1 targets that show similar expression trends across seed development. Gene expression data were retrieved from Belmonte et al. (2013). (d) Overlap between PHE1 targets, genes that were found to be significantly upregulated in osd1 when compared to wild-type (wt) seeds, and genes that were found to be significantly downregulated in phe1 phe2 osd1 when compared to osd1 seeds. P-values were determined using hypergeometric tests. (e) Expression of PHE1 targets that are significantly upregulated in osd1 seeds when compared to wt seeds. Genes marked in green are also significantly downregulated in phe1 phe2 osd1 seeds when compared to osd1 seeds.

Figure 1—figure supplement 1.. Examples of ChIP-seq read, peak, and motif distributions in PHE1 targets.

For each target gene, the read coverage for Input (as control) and PHE1–GFP ChIP samples of each replicate is shown. Called peaks in each replicate are indicated by the grey boxes. PHE1 binding sites, identified from the overlapped peak regions for replicates 1 and 2, are represented by the green boxes. PHE1 DNA-binding motifs are those defined in Figure 3a, and are represented here by the purple (Motif A) and pink boxes (Motif B). Transposable elements in the vicinity of the target gene are indicated. Red boxes correspond to TEs of the RC/Helitron superfamily, whereas white boxes represent TEs of other superfamilies. For this analysis, three different PHE1 targets with distinct imprinting statuses were selected: (a) AGL62 – non-imprinted; (b) MEDEA – a maternally expressed gene; (c) YUCCA10 – a paternally expressed gene.

Figure 1—figure supplement 2.. Schematic representation of the phe1 phe2 mutant.

CRISPR/Cas9 was used to generate a premature stop codon in PHE1, represented by the red asterisk. PHE1 mutagenesis was performed in the phe2 background because of the proximity and probable redundancy of the two PHE genes (see Materials and methods).

Figure 1—figure supplement 3.. Expression of PHE1 paralogs in osd1 vs. wt and phe1 phe2 osd1 vs. wt.

Paralogous genes were selected according to the phylogenetic analysis performed by Parenicová et al. (2003). All genes are significantly upregulated in osd1 vs. wt and in phe1 phe2 osd1 vs. wt.

Figure 2.. Transcription factor genes are enriched among PHE1 targets.

(a) Enriched biological processes associated with PHE1 target genes. Numbers on bars indicate number of PHE1 target genes within each GO term. (b) Enrichment of transcription factor (TF) families among PHE1 target genes (see Materials and methods). Numbers indicate the total number of Arabidopsis genes belonging to a certain TF family, and the total number of genes in that family targeted by PHE1. *, p-values <0.05. P-values were determined using the hypergeometric test.

Figure 3.. RC/Helitrons carry PHE1 DNA-binding motifs.

(a) CArG-box like DNA-binding motifs identified from PHE1 ChIP-seq data. (b) Fraction of PHE1 binding sites (green) that overlap transposable elements (TEs). Overlap is expressed as the percentage of total binding sites for which spatial intersection with features on the y-axis is observed. A set of random binding sites is used as control (blue). This control set was obtained by randomly shuffling the identified PHE1 binding sites within random A. thaliana gene promoters (see Materials and methods). P-values were determined using Monte Carlo permutation tests (see Materials and methods). Bars represent ± s.d. (n = 2494, for PHE1 binding sites and random binding sites). (c) Density of PHE1 DNA-binding motifs in different genomic regions of interest. P-values were determined using χ² tests. (d) Fold-change in the expression of genes flanked by RC/Helitron TEs. Fold-change was determined by comparing endosperm and embryo, or by comparing endosperm and seed coat. Genes were divided into four categories depending on their PHE1 target status, and the presence of RC/Helitrons with and without PHE1 DNA-binding motifs. Gene expression data were retrieved from Belmonte et al. (2013). Pre-globular seed stage was used in this analysis. P-values were determined using two-tailed Mann-Whitney tests (n = number represented below boxplots).

Figure 3—figure supplement 1.. PHE1 DNA-binding motifs show similarity to type II MADS-box CArG-boxes.

(a) Occurrence of PHE1-binding motifs in PHE1 binding sites. (b) Density of PHE1-binding motifs across PHE1 binding sites. (c) Alignment between PHE1-binding motifs and known motif matches. Both motifs match with previously characterized type II MADS-box CArG-boxes.

Figure 3—figure supplement 2.. RC/Helitron family overlap with PHE1 binding sites.

Fraction of PHE1 binding sites (green) overlapping with different RC/Helitron TE families. Overlap is expressed as the percentage of binding sites where spatial intersection with the families specified on the y-axis is observed. A set of random binding sites is used as control (blue). This control set was obtained by randomly shuffling the identified PHE1 binding sites within A. thaliana gene promoters (see Materials and methods). P-values were determined using Monte Carlo permutation tests (see Materials and methods). Bars represent ± s.d., (n = 2494, for PHE1 binding sites and random binding sites).

Figure 3—figure supplement 3.. Clustering of bound PHE1 DNA-binding motifs.

PHE1-DNA binding motifs and their flanking sequences were paired on the basis of sequence homology. The pairwise homologous sequences were then merged into higher-order clusters, which were based on shared elements in the homologous pairs (see Materials and methods). Motifs are those identified in Figure 3a.

Figure 3—figure supplement 4.. PHE1 DNA-binding motif densities in RC/Helitron consensus sequences.

(a) Density of perfect PHE1 DNA-binding motifs in the consensus sequences of different RC/Helitron families. Perfect motifs are those described in Figure 3a. As a reference, the density of perfect PHE1 DNA-binding motifs in the consensus sequences of all TE families is shown. (b) Density of nearly perfect PHE1 DNA-binding motifs in the consensus sequences of different RC/Helitron families. Nearly perfect motifs are those sequences where only one nucleotide substitution is required to generate a perfect PHE1 DNA-binding motif. As a reference, the density of nearly perfect PHE1 DNA-binding motifs in the consensus sequences of all TE families is shown. P-values were determined using χ² tests.

Figure 4.. Parental asymmetry of epigenetic marks in imprinted gene promoters conditions PHE1 binding.

(a) Fraction of non-imprinted and published imprinted genes targeted by PHE1. P-values were determined using the hypergeometric test. The list of published imprinted genes used for this analysis is detailed in Figure 4—source data 1. (b) Heatmap of endosperm H3K27me3 distribution along PHE1-binding sites. Each horizontal line represents one binding site. Clusters were defined on the basis of the pattern of H3K27me3 distribution (see Materials and methods) (c) Metagene plot of average maternal (♀, pink), paternal (♂, blue) and total (grey) endosperm H3K27me3 marks along PHE1 binding sites. (d) CG methylation levels in maternal (♀, upper panel) and paternal (♂, lower panel) alleles of PHE1 binding sites associated with MEGs (yellow), PEGs (green) and non-imprinted (grey) PHE1 targets. P-values were determined using two-tailed Mann-Whitney tests. (e) Sanger sequencing of imprinted and non-imprinted gene promoters bound by PHE1. SNPs for maternal (Ler) and paternal (Col) alleles are shown (n = 1 biological replicate). Maternal:total read ratios for imprinted genes are as follows: AT1G55650 – 1.0; AT2G28890 – 0.97; AT3G18550 – 0.44; and AT2G20160 – 0.32.

Figure 4—figure supplement 1.. Directionality of H3K27me3 distribution in Cluster 1 binding sites.

(a) Binding sites and their respective target genes can be present in the Watson (+) or Crick (–) strand of DNA. (b) Metagene plots of H3K27me3 distribution in PHE1binding sites associated with target genes located in the Watson (orange line) or Crick strands (blue lines). For this analysis, only H3K27me3 distributions in PHE1-binding sites associated with Cluster 1 were considered (Figure 4b). Total, maternal, and paternal H3K27me3 is represented in the top, middle, and bottom panels, respectively.

Figure 4—figure supplement 2.. Characterization of H3K27me3 clusters.

(a) Parental gene expression ratio of PHE1 targets associated with Cluster 1 and Cluster 2 binding sites, and all endosperm expressed genes. P-values were determined using two-tailed Mann-Whitney tests. (b) Fraction of Cluster 1 and Cluster 2 binding sites that target PEGs (Figure 4—source data 1) and putative PEGs (Moreno-Romero et al., 2019). (a, b) H3K27me3 clusters are those defined in Figure 4b.

Figure 4—figure supplement 3.. Parental-specific PHE1 ChIP.

qPCR of purified ChIP-DNA. Enrichment is shown as % of input DNA, in regions associated with MEGs (pink), PEGs (blue), and non-imprinted (grey) PHE1 target genes. Bars represent ± s.d. Data from one representative biological replicate are shown (n = 2 biological replicates).

Figure 5.. Ancestral RC/Helitron insertions are associated with gain of imprinting in the Brassicaceae.

(a–e) Phylogenetic analyses of PHE1-targeted PEGs and their homologs. Each panel represents a distinct target gene and its corresponding homologs in different species. The genes shown on a grey background have homologous RC/Helitron sequences in their promoter region. The arrow indicates the putative insertion of an ancestral RC/Helitron. The identity of the RC/Helitron identified in A. thaliana is indicated. These A. thaliana RC/Helitrons contain a PHE1 DNA-binding motif and are associated with a PHE1 binding site. The inset boxes represent the alignment between the A. thaliana PHE1 DNA-binding motif and similar DNA motifs contained in RC/Helitrons that are present in the promoter regions of orthologous genes. When available, the imprinting status of a given gene is indicated by the presence of ♂ (PEG) or ⚥ (non-imprinted), and reflects the original imprinting analyses done in the source publications (see Materials and methods). The maternal:total read ratio (M/T) for each gene is also indicated. §: potential contamination from maternal tissue. *: accession-biased expression. The scale bars represent the frequency of substitutions per site for the ML tree. The tree is unrooted. Gene identifier nomenclatures: AT, Arabidopsis thaliana; AL, Arabidopsis lyrata; Araha, Arabidopsis halleri; Bostr, Boechera stricta; Carubv, Capsella rubella; Cagra, Capsella grandiflora; Tp, Schrenkiella parvula; SI, Sisymbrium irio; Bol, Brassica oleracea; Brapa, Brassica rapa; Thhalv, Eutrema salsugineum; AA, Aethionema arabicum; THA, Tarenaya hassleriana; Cpa, Carica papaya; TCA, Theobroma cacao; Gorai, Gossypium raimondii; RCO, Ricinus communis; FVE, Fragaria vesca; and GSVIVG, Vitis vinifera.

Figure 5—figure supplement 1.. RC/Helitrons carrying PHE1 DNA-binding motifs are more prevalent in PEGs.

Fraction of PHE1 binding sites targeting MEGs, PEGs, or non-imprinted genes where a spatial overlap between a RC/Helitron and a PHE1 DNA-binding motif is observed. P-values were determined using the hypergeometric test.

Figure 5—figure supplement 2.. Analysis of PHE1-targeted PEG orthologs and upstream RC/Helitron sequences within Brassicaceae.

Alignments between A. thaliana RC/Helitrons associated with PEGs (Figure 5a–e) and homologous sequences in other Brassicaceae, where the putative homologous PHE1 DNA-binding motifs are present across all six species. For each aligned fragment, the percent sequence identity is indicated. The distances between the binding motifs and the TSS of PEGs are labelled. *, the distance to the TSS, in cases where the annotation for the 5′-UTR of the gene is not available.

Figure 6.. PHE1 establishes 3x seed inviability of paternal excess crosses.

(a) Target status of upregulated genes in paternal excess crosses. Highly upregulated genes in 3x seeds are more often targeted by PHE1 (p<2.2e–16, χ² test). n = numbers on the left. (b, c) Seed inviability phenotype (b) of paternal excess crosses in wild-type (wt), phe1 phe2, and phe1 complementation lines, with their respective seed germination rates (c). The maternal parent is always indicated first. Remaining control crosses are shown in Figure 6—figure supplement 1a–b. n = numbers on top of bars (seeds). (d) Parental expression ratio of imprinted genes in the endosperm of 2x (white) and 3x seeds (grey). Solid lines indicate the ratio thresholds for the definition of MEGs and PEGs in 2x and 3x seeds. (e) Accumulation of H3K27me3 across maternal (♀) and paternal (♂) gene bodies of PEGs in the endosperm of 2x and 3x seeds (white and grey, respectively). H3K27me3 accumulation in MEGs and non-imprinted genes is shown in Figure 6—figure supplement 2. P-value was determined using a two-tailed Mann-Whitney test.

Figure 6—figure supplement 1.. Rescue of 3x seed inviability in phe1 phe2.

(a, b) Seed inviability phenotype of wild-type (wt) and paternal excess crosses in wt, phe1 phe2, and phe1 complementation lines (a), with respective seed germination rates (b). (c, d) Status of endosperm cellularization in wt and paternal excess seeds. (a, b, d) n = numbers on top of bars (seeds).

Figure 6—figure supplement 2.. Distribution of H3K27me3 in 2x and 3x seeds.

Accumulation of H3K27me3 in maternal (♀) and paternal (♂) gene bodies of MEGs and non-imprinted genes in the endosperm of 2x and 3x seeds (white and grey, respectively).

Figure 7.. Control of imprinted gene expression by PHE1.

Schematic model of imprinted gene control by PHE1. Maternally expressed genes (left panel) show DNA hypermethylation of paternal (♂) PHE1 binding sites. This precludes PHE1 accessibility to paternal alleles, leading to predominant binding and transcription from maternal (♀) alleles. In paternally expressed genes (right panel), RC/Helitrons found in flanking regions carry PHE1 DNA-binding motifs, allowing PHE1 binding. The paternal PHE1 binding site is devoid of repressive H3K27me3, facilitating the binding of PHE1 and transcription of this allele. H3K27me3 accumulates at the flanks of maternal PHE1 binding sites, whereas the binding sites remain devoid of this repressive mark. PHE1 is able to bind maternal alleles, but fails to induce transcription. We hypothesize that the accessibility of maternal PHE1binding sites might be important for deposition of H3K27me3 during central cell development (possibly by another type I MADS-box transcription factor). It may also be required for the maintenance of H3K27me3 during endosperm proliferation.

Author response image 1.