Hybridization alters maternal and paternal genome contributions to early plant embryogenesis

In animals, early embryonic development is supported by maternally deposited products and does not depend on zygotic transcription (Vastenhouw, Cao, and Lipshitz, 2019; O'Farrell, 2015). In plants, initial experiments to study the onset of transcription after fertilization in wheat, maize and tobacco found novel transcripts in cDNA libraries of late zygotes compared with the egg cell (Sauter et al., 1998; Okamoto et al., 2005; Sprunck et al., 2005; Ning et al., 2006; Zhao et al., 2011). RNA sequencing experiments comparing gametes with isogenic embryos in maize, rice and Arabidopsis have provided more evidence for the transcription of a portion of the zygotic genome within hours after fertilization (haf). In maize zygotes at 12 h after pollination (hap), 3605 genes were found to be upregulated compared with their expression in the egg and sperm (Chen et al., 2017). In rice, comparison of zygote and egg cell transcriptomes found 499 genes upregulated at least twofold in zygotes at 2.5 hap, 1981 genes at 5 hap, and 2485 genes at 9 hap (Anderson et al., 2017). Arabidopsis zygotes at 14 hap (approximately 6 haf) showed 2625 genes upregulated at least twofold compared with egg cells, whereas at 24 hap, 2951 genes were upregulated compared with egg cells (Zhao et al., 2019). Thus, transcriptomic analyses of isogenic zygotes of several species show transcriptional activity for about 10% of the genome. Immunostaining experiments against the active form of Pol II have detected signals for active transcription in some experiments on plant zygotes (Niedojadło et al., 2012; Kao and Nodine, 2019) but not in others (Pillot et al., 2010), perhaps because of differences in the species studied, the developmental stage of the zygotes analyzed, or immunological techniques. In summary, compared with animals, most of the evidence in plants indicates that transcriptional activation of at least a small portion of the genome occurs in the zygote.

Initial experiments on maternal and paternal contributions to early isogenic embryos relied on reporter genes and functional analysis of early-acting embryo defective (emb) mutants (herein, emb refers collectively to these mutants). Assays of maternally and paternally inherited reporter lines in Arabidopsis found almost 30 genes expressed primarily from the maternal allele (Vielle-Calzada et al., 2000; Autran et al., 2011; Del Toro-De León et al., 2014), as well as ten genes with more equal maternal and paternal expression or function (Weijers et al., 2001; Lukowitz et al., 2004; Xu et al., 2005; Andreuzza et al., 2009; Aw et al., 2010; Ueda et al., 2011; Guo et al., 2016; Yang et al., 2017). In Arabidopsis, the functions of hundreds of EMBRYO DEFECTIVE (EMB) genes are required in early embryogenesis (Meinke, 2020). Crosses of 49 emb/+ mothers to isogenic wild-type (WT) plants showed that paternal alleles for nine EMB genes immediately complemented the lack of maternal function, whereas paternal alleles for 40 EMB genes did not fully complement until 3-5 days after pollination (dap) (Del Toro-De León et al., 2014). Consistent with the RNA sequencing experiments mentioned above, these marker gene and functional genetic experiments in isogenic embryos show that zygotic activation varies between genes, with many paternal alleles displaying less activity than maternal alleles in the first 2-3 dap.

Maternal and paternal transcripts in embryos can be quantified using single-nucleotide polymorphisms generated by crossing polymorphic ecotypes. Initial allele-specific expression (ASE) RNA-sequencing (RNAseq) experiments in Arabidopsis found evidence for a transient maternal bias in early Landsberg erecta (Ler) × Columbia-0 (Col) embryos (the ‘×’ denotes that a cross in only one direction was used) (Autran et al., 2011), and for early and equal maternal and paternal contributions in Col/Cape Verde Islands-0 (Cvi) embryos (the ‘/’ denotes an experiment with crosses in both directions) (Nodine and Bartel, 2012). Recent ASE RNAseq experiments of hybrid zygotes have also included data for gametes, allowing for comparison of the egg, sperm and zygote transcriptomes. In rice, more than 97% of 14,049 genes had a strong maternal bias at both 2.5 and 9 hap, with maternal transcripts for most genes found in both the egg and zygote (Anderson et al., 2017). A second experiment with Arabidopsis Col/Ler hybrid zygotes had a moderate maternal bias at 14 hap (9.9% of 12,746 genes showed maternally biased transcripts), but no maternal bias by 24 hap (Zhao et al., 2019). Approximately 60% of genes upregulated at 24 hap were reported as biallelic, suggesting that maternal and paternal alleles of many genes are equally transcribed in late Col/Ler zygotes (Zhao et al., 2019). The different conclusions of these parent-of-origin studies have been ascribed to profiling of mRNA versus total RNA (Baroux et al., 2013), different hybrid combinations used for ASE experiments (Baroux et al., 2013; Del Toro-De León et al., 2014, 2016) and sporophytic maternal contamination (Schon and Nodine, 2017). Differences between Arabidopsis and rice may also reflect variability in the timing of parental contributions between different species. In any case, the extent of maternal and zygotic regulation of early embryogenesis requires further clarification (Baroux and Grossniklaus, 2015; Armenta-Medina and Gillmor, 2019; Dresselhaus and Jürgens, 2021).

To determine whether varying effects of hybridization between different combinations of ecotypes might explain the divergent conclusions on zygotic genome activation previously reached in Arabidopsis, we explored whether hybridization might alter the maternal effects of emb mutants. Twenty emb/+ mutants (listed in Table S1) in Col were crossed with fathers of the Cvi, Ler, Tsushima-1 (Tsu) and Burren-0 (Bur) ecotypes. emb/+×Cvi and emb/+×Ler hybrids showed decreased maternal effects compared with isogenic emb/+×Col crosses, whereas maternal effects in emb/+×Tsu and emb/+×Bur hybrids were indistinguishable from isogenic emb/+×Col. We then used consistent criteria to reanalyze published transcriptome data for Col/Cvi (Nodine and Bartel, 2012) and Col/Ler (Zhao et al., 2019, 2020) zygotes and embryos. Although our conclusions were consistent with the original analysis of the Col/Cvi hybrid (Nodine and Bartel, 2012), for the Col/Ler hybrid, we found that plotting of gene reads as the log₂(fold change) (log₂FC) greatly diminished the representation of monoallelic expression in the zygote and basal cell lineage samples, resulting in an incomplete picture of parent-of-origin expression. Genes in the zygote and basal lineage previously proposed to exhibit allele dominance (Zhao et al., 2019, 2020) were more accurately interpreted as showing monoallelic bias in one direction of the cross and biallelic expression in the other direction. A similar uniparental distribution was observed among reads mapping to intronic regions, strongly suggesting that transcription of the maternal and paternal alleles in the zygote and basal lineages of Col/Ler zygotes is unequal. Taken together, our findings show that the relative parental contributions to zygotic genome activation can be altered by hybridization, and that zygotic genome activation in the Col/Ler hybrid occurs in an asymmetric manner.

Maternal and paternal gametophytic effects occur when the phenotype of the embryo depends on the genotype of the egg or sperm, respectively (reviewed by Armenta-Medina and Gillmor, 2019). We previously showed that 11 emb mutants exhibited transient maternal effects during early embryogenesis, but no paternal effects, and that crossing these emb mutants in the Col background with the Cvi and Ler ecotypes decreased their maternal effects (Del Toro-De León et al., 2014). To extend this analysis, we tested 20 additional emb mutants for maternal and paternal effects in early embryogenesis (Fig. 1). Phenotypes of embryos derived from hand-selfed Col emb/+ plants, from emb/+ mothers crossed with WT Col fathers, from WT Col mothers crossed with emb/+ fathers, and control crosses of Col mothers with Col fathers, were scored at 2, 3, 5 and 14 dap.

For all 20 emb/+ × Col crosses, maternal effects were observed beginning at 2 or 3 dap, which typically decreased by 5 dap, and were absent or almost absent by 14 dap (Fig. 1B; Table S1; supplementary Dataset 1). To determine whether any of the 20 emb mutants also had a paternal effect, we used WT Col plants as mothers in crosses with emb/+ fathers. Paternal effects were less common than maternal effects (Fig. 1B; Table S1). At 2 and/or 3 dap, nse3, pect1, iyo and miro1 showed transient paternal effects that were statistically indistinguishable from their transient maternal effect, whereas nse1, zyg3, qqt2, gex1 and fac1 showed transient paternal effects weaker than their transient maternal effects. The remaining 11 emb mutants showed no paternal effect at all (Fig. 1B; Table S1). In summary, maternal effects were seen for all 20 emb mutants; maternal and paternal effects of a similar degree were seen for four out of 20 emb mutants; and a paternal effect weaker than the maternal effect was seen for five out of 20 emb mutants. These results show that for 16 out of 20 EMB genes, the maternal allele plays a more important functional role in early embryogenesis than the paternal allele.

To explore the effect of hybridization on the maternal effects of these 20 emb mutants, we crossed emb/+ mothers in Col with the Cvi, Ler, Tsu and Bur ecotypes as fathers, and scored mutant phenotypes at 2, 3 and 5 dap (Fig. 1C; Table S1). As controls, Col mothers were crossed with Cvi, Ler, Tsu and Bur and scored for morphologically abnormal phenotypes (Table S1). The percentage of phenotypically abnormal embryos scored for each emb/+ hybrid cross at each timepoint, normalized to the percentage of abnormal embryos for the corresponding emb/+×Col cross, is shown in Fig. 1C. Compared with isogenic (emb/+×Col) crosses, emb/+×Cvi and emb/+×Ler hybrids showed significantly smaller maternal effects at 2, 3 and 5 dap, whereas maternal effects in emb/+×Tsu and emb/+×Bur hybrids were statistically indistinguishable from isogenic crosses. The decreased number of maternal effects observed in the progeny of emb/+×Cvi and emb/+×Ler demonstrates that paternal alleles for many of these 20 EMB genes are active earlier in these two hybrids compared with isogenic Col.

The results above show that maternal effects of 20 emb mutants are lower in Col×Cvi and Col×Ler hybrids than in isogenic Col. Nevertheless, the hypothesis that genomic parent-of-origin contributions to zygotes and early embryos can differ between isogenic Col, Col/Cvi and Col/Ler hybrids needs to be explored further (Del Toro De León et al., 2016; Armenta-Medina and Gillmor, 2019; Dresselhaus and Jürgens, 2021). To compare Col/Cvi and Col/Ler ASE RNAseq datasets as directly as possible, we comprehensively remapped and reanalyzed data from Nodine and Bartel (2012) and Zhao et al. (2019, 2020).

ASE RNAseq data from Nodine and Bartel (2012) and Zhao et al. (2019, 2020) were mapped to their parent of origin (Table S2) and are plotted in Fig. 2 as line graphs in two different ways: as probability density versus log₂FC ratio of the parental ecotypes (as originally plotted in these three studies) (Fig. 2A-D); and as probability density versus maternal fraction of reads per gene (see Materials and Methods) (Fig. 2E-H). As previously shown, log₂FC representation of Col/Cvi datasets produced peaks centered around 0 (biallelic) for all stages (Fig. 2A) (Nodine and Bartel, 2012), whereas log₂FC representation of Col/Ler datasets showed peaks with slight maternal biases in 14 hap zygotes, which shifted to 0 (biallelic) in 24 hap zygotes (Fig. 2B) and stayed centered around 0 in subsequent cell types and developmental stages (Fig. 2C,D) (Zhao et al., 2019, 2020). Based on log₂FC ratios, Nodine and Bartel (2012) concluded that one- to two-, eight- and 32-cell Col/Cvi embryos had equal parental contributions to the transcriptome, and Zhao et al. (2019, 2020) concluded that Col/Ler zygotes had a maternal transcript bias at 14 hap, shifting to equal parental contributions in 24 hap zygotes and at later stages.

Fig. 2E-H shows the same datasets, represented as probability density versus maternal fraction of reads per gene. To avoid confounding effects by low-count genes that could easily accumulate in extreme peaks, strong filters were employed and the estimated maternal fraction for each gene was calculated across biological replicates using edgeR (see Materials and Methods). For the Col/Cvi dataset, plotting as maternal fraction of reads led to a similar interpretation as for log₂FC plots: small peaks near 0 and 1 emerged, representing genes that are exclusively paternal or exclusively maternal, but the majority of reads continued to be clustered around 0.5, i.e. genes with equal maternal and paternal contributions (Fig. 2E). By contrast, representation of Col/Ler datasets as maternal fraction of reads had a drastic effect on their interpretation, as a trimodal distribution emerged. At 14 hap, zygote samples had central peaks located at about 0.65 (Fig. 2F), whereas at 24 hap, this central peak was centered around 0.5 (Fig. 2F). At 24 hap, the Col×Ler central peak was considerably smaller than the corresponding Ler×Col peak, demonstrating the asymmetric effect of hybridization of these ecotypes. In addition to the central peaks, both Col×Ler and Ler×Col zygotes showed large peaks near 0 (completely paternal) and 1 (completely maternal) (Fig. 2F). In the Col×Ler cross, 13% (n=12,342) of genes at 14 hap and 24% (n=10,994) of genes at 24 hap showed more than 95% bias towards either parent. In the Ler×Col cross, 16% (n=12,122) of genes at 14 hap and 7% (n=11,000) of genes at 24 hap showed more than 95% bias towards either parent (values derived from Table S3). The strong filters used in our analysis indicated that the maternal and paternal peaks were not caused by sampling biases of genes with low counts.

The trimodal distribution observed in Col/Ler zygotes and embryos continued in the basal cell produced by asymmetric division of the zygote (Fig. 2G), and to a lesser extent in the suspensor of the globular embryo, which is produced by the basal cell (Fig. 2H). Interestingly, the apical cell resulting from the asymmetric division of the zygote had an almost monomodal distribution, with most reads clustered around 0.5, and only small peaks at 0 and 1 (Fig. 2G). This trend continued in the globular embryo, for which the peaks at 0 and 1 had essentially disappeared and almost all reads were clustered around 0.5 (equal parental contributions) (Fig. 2H). Thus, after the first division of the zygote, the asymmetric transcriptional state of the Col/Ler zygote is reset in the apical cell, yet continues in the basal cell lineage, decreasing during subsequent rounds of division of the suspensor lineage.

The apparent scarcity of genes with strong maternal and paternal transcript bias in the analyses by Zhao et al. (2019, 2020) can be explained by handling of biological replicates and log₂FC representation of data. Previous analyses of this data (Zhao et al., 2019, 2020) appeared to rely on adding the counts across biological replicates and then calculating log₂FC without properly handling zeroes (table S5 in Zhao et al., 2019; table S1 in Zhao et al., 2020). Zeros frequently arise due to lack of sampling or sequencing depth and are commonly dealt with by adding a small prior count to the numerator and denominator before calculating the fold change (Chen et al., 2016), although there are more sophisticated statistical frameworks (Love et al., 2014; Erhard, 2018). The methods of Zhao et al. resulted in a lack of consideration of biological variability across replicates and the generation of infinite log₂FC values when the counts came only from one parent, artificially excluding completely monoallelic genes in their log₂FC plots. Log₂FC representation also causes genes that tend towards unbiased behavior to cluster around 0, increasing their apparent frequency; and genes with biased parent-of-origin behaviors to spread out along the x-axis, minimizing their apparent frequency (Fig. 3). For example, a gene with one paternal and 100 maternal reads, and a gene with one paternal and 1000 maternal reads would both show almost the same maternal fraction (0.99 versus 0.999), but highly divergent log₂FC values (−6.6 versus −10). Thus, for the purpose of analyzing and interpreting parent-of-origin contributions on a gene-by-gene basis, we believe that representation of data as maternal fraction of reads is more useful than log₂FC representation.

To determine whether the allelic imbalance shown by gene transcripts in Fig. 2E-H was due to ecotype-dominant effects, parent-of-origin effects or a mixture of both, we used the exact test for negative binomial distributions in edgeR to calculate statistically significant deviations from the expected 1:1 ratio for each direction of the cross and for all samples (see Materials and Methods; Table S3). In agreement with Zhao et al., (2019), 16.3% of 13,540 genes in 14 hap zygotes showed maternal bias in both directions of the cross, and by 24 hap, this category had almost completely disappeared (Fig. 2M). In contrast to the analysis of Zhao et al. (2019, 2020) that concluded that 24 hap Col/Ler zygotes showed equal parental transcriptome contributions, in our analysis, thousands of genes showed maternal or paternal bias in one direction of the cross, and biallelic expression in the other direction, a category that we refer to as unidirectional bias. In 14 hap zygotes, 42.5% of 13,540 of genes showed unidirectional bias, whereas at 24 hap, 38.6% of 12,389 genes showed unidirectional bias, indicating that thousands of genes in Col/Ler zygotes contribute transcripts primarily from the maternal allele or paternal allele (Fig. 2M; Table S4). By contrast, less than 3.1% of genes showed Col or Ler dominant behavior at any stage (Fig. 2M; Table S4). It had previously been noted that a significant portion of genes showed monoallelic expression in the basal cell lineage of the Col/Ler hybrid (Zhao et al., 2020). These genes were previously described as ecotype dominant because their maternal-to-paternal ratio was found to be negatively correlated between reciprocal crosses. In contrast to the conclusion of Zhao et al. (2019) that 24 hap Col/Ler zygotes show equal parental contributions, our analysis shows that monoallelic expression in zygotes and the basal lineage is common and is primarily caused by genes that are not reciprocally biased towards any ecotype or parent but are rather unidirectionally biased, i.e. that parental contributions to the Col/Ler zygotic transcriptome vary greatly between genes and are thus unequal.

As shown in Fig. 2F-H, for a given cell type, a similar proportion of maternally or paternally biased genes exists in both directions of the cross in the zygote and basal cell lineage. Most of these genes are unidirectionally biased (Fig. 2M), i.e. they are only significantly biased in one direction. This suggests that the regulation of allelic expression in the zygote and embryo is strongly affected by an interaction between the ecotype and parent-of-origin effects. Considering that a comparable number of unidirectionally biased genes were discovered across zygote, basal cell and globular-suspensor samples (Table S4), we wondered whether cell identity also plays a role in regulating allele-specific expression. To test this, we examined the overlap of unidirectional maternally and paternally biased genes between 14 hap and 24 hap zygotes, the basal cell, and the suspensor, for each direction of the cross (Fig. 4A-D; Table S5). This analysis showed that allele bias for most genes is not preserved through developmentally adjacent samples. For example, of 4325 genes with maternal bias in the Col×Ler 14 hap zygote, only 692 also showed maternal bias in the 24 hap zygote (Fig. 4A). Of 4989 genes with maternal bias in the Ler×Col 14 hap zygote, only 506 also showed maternal bias in the 24 hap zygote (Fig. 4B). A similar lack of overlap of paternally biased genes was seen between 14 and 24 hap zygotes (Fig. 4C,D). Interestingly, although the total number of genes with maternal bias decreased from 14 to 24 hap, the number of paternally biased genes increased from 14 to 24 hap, perhaps reflecting a general increase in transcription of the paternal genome between 14 and 24 hap. Thus, during early development of zygotes and the basal cell lineage of Col/Ler embryos, allele-specific contributions are continuously restructured in a gene-specific manner. Taken together, these analyses contradict the conclusion that 24 hap zygotes and one-cell embryos show equal parental transcript contributions (Zhao et al., 2019, 2020), and show that the biological interactions between ecotypes that regulate parental contributions to Col×Ler and Ler×Col transcriptomes have previously been oversimplified.

After the first zygotic division, a significant portion of the transcriptome of the basal cell lineage remained monoallelic (Fig. 2G,H), even though parent-of-origin transcript bias was not maintained through cell divisions of the zygote and basal lineage (Fig. 4A,B). To test whether zygotic transcription of maternal or paternal alleles could explain the biases seen for thousands of gene transcripts in the basal cell lineage (Fig. 2F-H,M), we first quantified parent-of-origin transcripts for genes present in the zygote but absent in the egg cell. Because the majority of transcripts expressed in the egg cell are also found in the zygote (Zhao et al., 2019), excluding genes with egg-cell transcripts removed most zygotic genes, limiting the set of de novo parent-of-origin transcripts identified with this approach to a few hundred genes (Fig. S1; Table S6). In 14 hap Col×Ler zygotes, 179 of 251 genes showed paternal bias, with 80 of these having a maternal fraction of less than 0.1. In 24 hap Col×Ler zygotes, 199 of 326 genes showed paternal bias, with 66 having a maternal fraction less than 0.1. Although the paternal bias in de novo transcripts might indicate a paternal bias in transcription in the zygote, the lack of maternally biased genes could also be due to overlap in transcription of maternal alleles between the egg cell and zygote, limiting the conclusions that can be made from identification of de novo expressed genes.

As a more comprehensive method for characterizing de novo transcription, we then quantified reads mapping to variants in introns. Transcripts carried over from gametes should be spliced, whereas active de novo transcription leads to a small proportion of reads mapping to introns of immature intron-containing transcripts, which serve as a proxy for active transcription and can be detected in regular RNAseq experiments (Gaidatzis et al., 2015; Lee et al., 2013). Although reads derived from introns were less abundant than exon reads (as expected), our methods allowed us to identify statistically significant parent-of-origin behavior for intron reads for thousands of genes (Fig. 2I-L,N; Tables S7-S9).

In the 14 hap zygote samples, intron reads were heavily biased towards maternal alleles (Fig. 2J), suggesting that maternal alleles of most genes gain transcriptional competence before paternal alleles, and that active transcription is partially responsible for the maternal bias seen in whole-gene reads at this stage (Fig. 2F). In the 24 hap zygote, both parental genomes generated intron reads, but primarily at different loci, resulting in strong peaks of maternal and paternal gene transcripts when Col was used as a mother, and a more uniform distribution (with smaller peaks) when Ler was the mother (Fig. 2J). After the first asymmetric division, transcription in the basal lineage followed a pattern similar to the zygote, whereas transcription in the apical cell showed a strong central peak and was essentially biallelic by the globular embryo stage (Fig. 2K,L).

Compared with the Zhao et al. (2019, 2020) datasets, the Nodine and Bartel (2012) dataset had relatively few reads mapping to introns, but some trends were nonetheless evident. Except for the octant-stage Cvi×Col transcriptome, all other Col/Cvi datasets showed large peaks for the maternal and paternal fractions, with about 30% of genes showing equal maternal and paternal reads (Fig. 2I). Intriguingly, this was in contrast with the distribution seen using reads mapping to variants in whole genes, in which these monoallelic tendencies were much less frequent (Fig. 2E). Plots for active transcription in Col/Ler zygotes also showed increased bias towards one parent or the other (Fig. 2K,L) compared with total reads for the same stages (Fig. 2G,H). The difference between active transcription and steady-state transcripts suggests an important role for post-transcriptional regulation in parental mRNA contributions during early embryogenesis. It is also possible that gene expression occurs in waves of active transcription from each genome with allelic ratios partially balancing in the cytoplasm. Disparity between active transcription and cytoplasmic mRNA levels could also be affected by carryover of spliced mRNA from progenitor cells.

We performed a statistical analysis to examine whether the trends in intron reads were a mixture of parent-of-origin and ecotype effects (Table S8), as seen with whole-gene reads. In 14 hap zygotes, 60.9% of genes were biallelic, 34.6% showed unidirectional bias and 3% were maternal (n=4804), whereas in 24 hap zygotes, 47.3% were biallelic, 45% showed unidirectional bias and 1% were maternal (n=3113) (Fig. 2N; Table S9). After the first division of the zygote, 97.5% of genes in the apical cell were biallelic (n=1756) and only 1.8% of genes showed unidirectional bias, whereas in the basal cell, 85.1% of genes were biallelic and 14.3% of genes (n=2037) showed unidirectional bias. Except for maternally biased genes in 14 hap zygotes, which constituted 16.3% of genes using whole-gene reads and 3% using intron reads, most other classes of genes were very similar between whole-gene and intron reads at all stages and in all cell types, suggesting that the maternal reads in the 14 hap zygote are primarily inherited from the egg cell, with a smaller contribution from zygotic maternal transcription, and that active transcription is primarily responsible for determining transcript populations in the 24 hap zygote and at later stages. In the Col/Cvi hybrid, despite the marked bias towards uniparental expression (Fig. 2I), no statistically significant results were found (Tables S4 and S9). However, the lack of biological replicates for the Col/Cvi dataset and the general scarcity of the reads mapping to polymorphisms in introns makes it likely that this view is not complete. Taken together, this evidence suggests that active transcription is not equivalent for both parental genomes at early stages of embryogenesis, with the ecotypes used and the direction of the cross as key factors that determine allelic transcription.

Maternal effects occur when the phenotype of the zygote depends on the genotype of the mother (Armenta-Medina and Gillmor, 2019). Forty out of 49 emb mutants in the Col ecotype were previously shown to display transient maternal effects during embryogenesis, presumably due to early inactivity of the corresponding WT paternal allele (Del Toro De León et al., 2014). In this work, we confirmed transient maternal effects for 20 additional emb mutants in Col and in hybrids of Col/Cvi, Col/Ler, Col/Tsu and Col/Bur.

Comparison of maternal effects and transcriptome data for EMB genes in isogenic Col (Table S10) shed light on the molecular nature of the maternal effect for the MONOPTEROS (MP) and GAMETE EXPRESSED PROTEIN 1 (GEX1) genes. MP transcripts were not detected in the egg cell but were present at low levels in isogenic embryos beginning at the one-cell stage. Thus, the maternal effects seen for the mp mutant were likely to be zygotic maternal effects. By contrast, GEX1 transcripts were present at very high levels in the egg cell, decreased in the 14 hap zygote, and were almost absent later, suggesting that the maternal effects observed for gex1 were likely to be maternal gametophytic effects. Transcripts for the other EMB genes we tested were present in both egg cells and zygotes, so the maternal effects observed in these isogenic emb×Col crosses could be either gametophytic or zygotic.

For emb/+×Cvi and emb/+×Ler crosses, maternal effects were decreased compared with those of emb/+×Col, suggesting that activation of paternal alleles for these EMB genes occurs earlier in Col×Cvi and Col×Ler hybrids than in isogenic Col. Looking at parent-of-origin transcripts for individual EMB genes in Col×Ler zygotes and one-cell embryos, and in Col×Cvi one- to two-, eight- and 32-cell embryos (Table S10), all EMB genes for which there are data are biallelic at all stages, with the exception of ZYG3, which is paternal at 24 hap, and ATLA1 and IYO, which are maternal at 24 hap (Table S10). Biallelic zygotic and embryo transcriptomes in Col×Ler and Col×Cvi make sense with the reduced maternal effects for emb mutants in these hybrids compared with those in isogenic Col, and are consistent with the conclusions of Nodine and Bartel (2012) and Zhao et al. (2019) that transcriptional activation of paternal alleles in Col/Cvi occurs by the one- to two-cell embryo stage and in Col/Ler hybrids by the 24 hap zygote stage.

It should be noted that emb/+×Ler and emb/+×Cvi crosses do indeed show transient maternal effects, but of less magnitude than isogenic emb/+×Col and hybrid emb/+×Tsu and emb/+×Bur crosses. These maternal effects could originate from maternal transcripts carried over from the egg cell to the zygote and/or zygotic maternal transcripts. At the genome level, Col/Ler zygotes provide an example of how transient zygotic maternal effects could occur. As seen in Fig. 2J, 14 hap Col/Ler zygotes have a strong maternal bias in intron transcripts. This maternal zygotic transcriptional bias is short lived: after the first division of the zygote, the apical cell shows more equal maternal and paternal intron transcripts, whereas in the basal cell, most genes show either a strong maternal or paternal bias (Fig. 2K). The increased maternal effects of emb mutants in emb/+×Col compared with emb/+×Ler suggests that a maternal bias in zygotic transcription may last longer in isogenic Col than in Col/Ler zygotes. Another possibility is that gametophytic maternal effects are stronger in isogenic Col than in Col/Ler hybrids. In either case, variation in gene expression and transcript stability between isogenic and hybrid embryos may be due to differences in the timing of similar developmental mechanisms.

Previous studies with the Col/Ler hybrid found that 66% of informative reads in 14 hap zygotes mapped to the maternal genome, whereas in 24 hap zygotes and in both apical and basal cells, after the first zygotic division, maternal and paternal genome reads were balanced at about 50% each (Zhao et al., 2019, 2020). Genes upregulated at 14 hap in the zygote compared with those in the egg cell were primarily biallelic, whereas downregulated genes were most often maternally biased (Zhao et al., 2019). Zhao et al. (2019) concluded that a maternal effect exerted through mRNA carryover from the egg cell was present at 14 hap but was quickly abolished via mRNA degradation. Despite similar amounts of total reads mapping to the maternal and paternal genomes, 33% of genes in the apical cell and 20% of genes in the basal cell were interpreted as monoallelic (Zhao et al., 2020). The allelic bias for these genes was interpreted as negatively correlated with their bias in the reciprocal cross and, thus, these imbalances were attributed to allele-dominant effects, i.e. genes which are primarily expressed from the Col or Ler allele, regardless of the parent of origin.

In our analysis of parent-of-origin contributions in Col/Ler hybrids, it became clear that although the contributions of the maternal and paternal genomes considered as a whole are indeed similar, the parent-of-origin contributions of thousands of genes vary widely. For example, in 14 hap zygotes, 42.5% of genes showed parental bias in only one direction of the cross (Fig. 2M), whereas at subsequent developmental stages, most of these genes showed varying parent-of-origin contributions (Fig. 4A-D). A similar trend was observed in 24 hap zygotes (Fig. 2M, Fig. 4A-D). Thus, the assertion that the maternal and paternal genomes make equal contributions to 24 hap zygotes and the subsequent basal cell lineage (Zhao et al., 2019, 2020) applies only when reads for the entire genome are considered. These observations strongly suggest that, at least in the Col/Ler hybrid, parental contributions to the early Arabidopsis embryo are dynamic and shaped by several factors, including the parent of origin, the combination of parental ecotypes, and cell identity.

Given the unique epigenetic state of gametes and their dynamic remodeling after fertilization (Okada et al., 2005; Borg et al., 2020; Ingouff et al., 2017, 2010; Niedojadło et al., 2012), it would not be unexpected for the maternal and paternal genomes to undergo different or asynchronous reprogramming in the zygote. In isogenic Arabidopsis, maternally encoded histone 3 lysine 9 dimethylation (H3K9me2) as well as various components of the RNA-dependent DNA methylation pathway were shown to repress several paternal marker genes, whereas maternally encoded components of the chromatin assembly factor complex (CAF1) acted to promote the expression of paternal marker genes (Autran et al., 2011). In addition, POLYCOMB REPRESSIVE COMPLEX 2 (PRC2) has been shown to regulate the imprinted expression of several genes in Arabidopsis embryos (Raissig et al., 2013). In the case of hybrids, some Arabidopsis ecotype combinations may have divergent epigenetic configurations in their gametes, which could lead to changes in DNA methylation and gene expression due to siRNA-mediated interaction between the two genomes, as has been found in vegetative tissues of Arabidopsis (Groszmann et al., 2011; Greaves et al., 2012). Rice zygotes undergo significant changes in siRNA populations compared with gametes (Li et al., 2020, 2022) and, thus, varying epigenetic configurations between gametes could potentially alter the epigenetic pathways that may regulate parent-of-origin expression in the embryo, priming or silencing a different set of genes in each ecotype. Likewise, differences in parent-of-origin transcription between developmental stages (Fig. 4A-D) could be a consequence of an interaction between reconciliation of epigenetic marks in the gametes and zygote and developmental programs in the embryo.

After the first division of the zygote establishes the apical and basal lineages, parent-of-origin and ecotype effects almost disappear in the apical lineage (apical cell and globular embryo proper), which transitions to essentially equal maternal and paternal transcripts for most genes (Fig. 2G,H,M). However, strong maternal and paternal transcript peaks persist in the basal cell (Fig. 2G) and, to a lesser extent, in the suspensor of the globular embryo (Fig. 2H). In agreement with the analysis of gene reads, analysis of active transcription via intron reads suggests that many genes in the basal lineage are being transcribed mostly from either the maternal or paternal allele, whereas active transcription is mostly biallelic in the apical cell and globular embryo proper (Fig. 2K,L). Although thousands of genes in the suspensor showed an allelic ratio far from 1:1 or marked peaks of monoallelic expression, our analysis identified few genes in the suspensor with a statistically significant bias (Fig. 2M,N). This may be due to the increased variability in relative parental contributions for suspensor samples between biological replicates in comparison with other tissues (Fig. S2), a drawback of the conservative statistical approach we used to define parentally biased genes.

Why might asymmetric parent-of-origin expression persist longer in the basal lineage than in the apical lineage? One possibility is that inherent biological differences between the functions of the suspensor and the embryo cause differential parental genome activity in these two tissues. For example, the suspensor is a transient tissue that moves nutrients and hormones from the maternal tissues to the embryo proper (Peng and Sun, 2018) and, thus, the parental conflict theory would predict a high pressure to imprint genes in the suspensor (Pires and Grossniklaus, 2014). What is clear is that the basal cell and its lineage show much more widespread parental bias in transcript abundance than the apical cell lineage, and that the direction of the cross and ecotypes used strongly affect parental contributions to the embryo transcriptome.

Col/Cvi transcriptomes showed almost equal parental contributions starting at the one-to two-cell embryo stage, with only about 100 genes that show monoallelic behavior, suggesting quick degradation of maternally inherited products and synchronous activation of both parental genomes (this work; Nodine and Bartel, 2012). Thus, one- to two-cell and later embryos in the Col/Cvi hybrid almost completely lack the variable pattern of parental contributions seen in Col/Ler zygotes and one-cell embryos, such as unidirectionally biased transcripts in the zygote and monoallelic behaviors seen in the basal cell and globular suspensor. Although it is possible that sampling of whole one- to two-cell Col/Cvi embryos diluted any monoallelic patterns that might have existed in the basal lineage, the dramatic contrast between Col/Cvi and Col/Ler hybrids points to intrinsic differences between hybrid combinations in terms of the dynamics of zygotic genome activation (Baroux et al., 2013; Del Toro-De León et al., 2016).

One possibility to explain variation in parent-of-origin genome activation in Col/Ler and Col/Cvi hybrids might be the epigenetic differences between Cvi, Ler and Col gametes. Throughout male gametogenesis in Col, CpG methylation patterns remain stable (Calarco et al., 2012; Ingouff et al., 2017). However, the major CpG-DNA maintenance methyltransferase MET1 is barely detected in the mature egg cell of Col (Jullien et al., 2012), and the signal of a fluorescent reporter for CpG-methylation sharply decreased in the mature egg cell (Ingouff et al., 2017). Taking this into account, it is tempting to speculate that differences in CpG methylation in parental genomes might be partly responsible for the differential remodeling and transcriptional activation after fertilization, which could also explain how an accession like Cvi behaves so differently in these terms, considering that the Cvi ecotype has significantly lower CpG methylation in genes compared with Col and Ler (Schmitz et al., 2013; Pignatta et al., 2014; Kawakatsu et al., 2016). Compared to Col and Ler, Cvi also lacks HISTONE DEACETYLASE 6 (HDA6) function, shows decondensed chromatin and reduced DNA and H3K9 methylation at chromocenters (Tessadori et al., 2009), all of which are characteristic of active transcription and could impact genome activation in the zygote (Baroux et al., 2013).

In summary, our analysis shows that Col/Ler and Col/Cvi hybrids decrease the maternal effects of emb mutants compared with those of isogenic Col, pointing to Col/Tsu and Col/Bur hybrids as better proxies for understanding parent-of-origin behavior in Arabidopsis. Contrary to previous reports, the Col/Ler hybrid exhibits parent-of-origin transcriptome contributions to the zygote and basal lineage, which vary widely between genes and between crosses, demonstrating that on a gene-by-gene basis, the maternal and paternal genomes in Col/Ler hybrids make unequal contributions to the zygote at both 14 and 24 hap. By contrast, one- to two-, eight- and 32-cell Col/Cvi hybrid transcriptomes show almost no parent-of-origin behavior. The differing parent-of-origin transcriptome patterns and altered maternal effects of emb mutants in Col/Cvi and Col/Ler hybrids suggest that neither of these hybrid transcriptomes paints an accurate picture of parent-of-origin gene regulation in isogenic Columbia embryos. In the future, differences in zygotic genome activation between hybrids might be leveraged to better understand transcriptional regulation in plant zygotes and embryos. Although mechanisms that regulate parent-of-origin expression in the Arabidopsis embryo at the genome level remain to be described, studies of gene expression in the Arabidopsis endosperm have uncovered multiple epigenetic pathways regulating imprinting (Batista and Köhler, 2020). In embryos of the liverwort Marchantia polymorpha, polycomb-mediated histone 3 lysine 27 trimethylation is widely deposited on paternal chromosomes and is responsible for transcriptional silencing of hundreds of paternal alleles during the brief sporophytic phase of this lower plant (Montgomery et al., 2022). It will be interesting to determine whether epigenetic mechanisms for imprinting in the Arabidopsis endosperm or Marchantia embryo also regulate parent-of-origin gene expression in the Arabidopsis embryo.

In notation of crosses, the mother is always listed first. In referring to hybrid embryos, ‘Col×Ler’ refers to a hybrid cross in one direction, whereas ‘Col/Ler’ refers to hybrid crosses in both directions, i.e. Col×Ler and Ler×Col. emb/+ crosses were conducted as in Del Toro-De León et al. (2014). Seed stocks used for all crosses are listed in Table S1. For each cross at each timepoint for each emb mutant, seeds from at least three different siliques from three different plants were examined. Statistical significance was calculated using Fisher's two-tailed exact test for P<0.05. An emb mutant was considered to have a maternal effect when the proportion of mutants in the emb/+×Col cross was significantly different from the Col×Col control, and a paternal effect was called when the portion of mutants in the Col×emb cross was significantly different from the Col×Col control. Maternal and paternal effects for a particular emb mutant were labeled as indistinguishable when there was no statistically significant difference between the proportion of mutants observed between the emb/+×Col and Col×emb/+crosses (Table S1). For comparison of emb maternal effects between isogenic Col and hybrid crosses, mutant frequency for each gene was normalized by adding 0.02 and dividing by the frequency in isogenic Col for that gene plus 0.02. This addition was to stabilize mutants with low penetrance and was based on the mutant frequency observed in WT Col. A Wilcoxon signed rank test was then used with the null hypothesis that these values were centered at 1 (equivalent to isogenic Col) for each ecotype used for each developmental stage (2, 3, and 5 dap).

Custom genomes for Ler and Cvi were generated by replacing the variants from each ecotype (1001 Genomes Consortium, 2016) in the reference sequence of Col. The FASTA files were then each joined with the FASTA file from Col. These were used to map the reads. Reads were mapped using STAR v2.7.8.a (Dobin et al., 2013). The following flags were used: --alignIntronMax 900 --alignMatesGapMax 900 --outFilterMultimapNmax 1 --outFilterMultimapScoreRange 0. Two-pass mode was used with all libraries from each study used together for the first pass to better define the relevant splice junctions.

For the parent-of-origin-specific counting, for each hybrid, an annotation file that included only the variants in genes was generated by adjusting the Araport 11 annotation for indels, then each gene span was reduced to only the variant regions generating an entry per variant that overlapped a gene. These entries were then split into two categories: one annotated as located in the reference (Col) chromosome and one in the alternative (Cvi or Ler) chromosome. This annotation was put in the SAF format for use with featureCounts (see documentation for featureCounts; Liao et al., 2014).

Reads overlapping variants were counted using featureCounts v2.0.1 with the following flags: -O -p --fraction -F SAF. The fractionated counts option was added to avoid counting multiple times reads that overlapped more than one variant.

To analyze parental bias in intron reads, an annotation file was generated where only the segments of the genome that contained sequence variants and that are spliced out in all annotated isoforms in Araport 11 were considered. Reads were counted using this file with the same parameters as reads mapping to variants in whole genes.

To be included in differential expression analysis, genes were required to have at least 20 counts and at least 3 counts per million in all but one of the replicates available or in the only library for each Col/Cvi sample. For each library, the sum of parental counts was used to calculate library depths and normalization factors for edgeR. Then, for each library, maternal versus paternal counts were used to test for significant deviation from the expected 1:1 ratio using the exact test implemented in edgeR. Dispersion was calculated only among biological replicates with the tagwise option. As the Col/Cvi samples had no biological replicates, dispersion was set to 0.58 (the median value among samples). Fold changes were calculated using edgeR and were converted to maternal read fractions per gene. Estimated maternal fractions, false discovery rates (FDRs) and total variant counts are available in Table S3. The same data for intronic reads are available in Table S8.

To generate the logarithmic fold-change plots from Fig. 2, raw allelic counts for genes were filtered in the same way as for differential expression analysis. To replicate results from Zhao et al. (2019, 2020), logarithmic fold change was then calculated for the sum of counts across replicates. The maternal fraction of reads estimated using edgeR was used for the rest of the plots in Figs 2–4. To allow for comparison across samples in Fig. 2, the range of the yaxis was set to the maximum in each row. Probability densities for Fig. 2 were estimated with the geom_density() function in R, which produces a smoothed version of a frequency histogram, so that areas under all curves are equal to 1. Nodine and Bartel (2012) and Zhao et al. (2019, 2020) likely used a similar calculation for their density curves, although the method used to produce their curves, and the exact definition, were not specified.

All RNAseq data analyzed has been previously published. Col/Cvi Illumina data were from Nodine and Bartel (2012) [Gene Expression Omnibus (GEO) accession number GSE33713]; Col/Ler Illumina data for zygote samples were from Zhao et al., (2019), (GSE121003); and Col/Ler Illumina data for one-cell and globular stage embryos were from Zhao et al., (2020) (GSE107700).

Seeds of emb lines were obtained from the Arabidopsis Biological Resource Center (https://abrc.osu.edu Ohio State University, Columbus, OH, USA). We are grateful to David Meinke and colleagues for their years of effort assembling and maintaining the emb mutant collection, which served as the foundation of this work. We thank Daphné Autran, Sean Cutler and Raju Datla for critical reading and comments.

The peer review history is available online at https://journals.biologists.com/dev/lookup/doi/10.1242/dev.201025.reviewer-comments.pdf.