Boechera seeds (sexual and asexual) were generous donations of Bitty Roy (University of Oregon, USA) and Eric Schranz (Wageningen University, USA); seeds were originally isolated from the native habitat of Boechera in the USA (Supplementary Figure S1), and they represent green-house or growth-chamber grown bulked seeds when donated. The progeny were further screened by bulked and single seed flow cytometry to remove odd genotypes (4n, 6n, aneuploids), and 3n and 2n plants of Apo-1 and Apo-2 were identified and characterized by cytology (unpublished data). Subsequent generation of seeds were surface-sterilized, stratified at 4°C and grown on half-strength MS agar plates for 17–18 days. For reproductive development, they were transferred subsequently to soil (Profi-substrate, Einheitserde, Germany) and quartz sand in a 4:1 combination. Seedlings and plants were grown in dedicated plant growth chambers with long day photoperiod (16-h light/8-h darkness cycle), under cool white light with ca. 120 μmol m^-2s^-1 intensity, 60–65% air humidity and 22°C. Plants were ploidy genotyped either in a Partec flow cytometer (Sysmex, Germany) or in a Guava^® easyCyte flow cytometer (Millipore, USA).

Seventeen to eighteen days old seedlings from Sex-1 and Apo-1 accessions were subjected to mRNA sequencing. Seedlings were snap-frozen and kept at -80°C until RNA extraction; tissue was ground using a tissue lyser (Retsch, Germany). RNA extraction was performed using QIAGEN RNeasy Mini Kit and DNase–treated according to manufacturer’s instructions (Qiagen, USA). The RNA quality was further validated in a Bioanalyzer (Life Technologies, USA). mRNA transcripts purification, cDNA synthesis, library preparation and NGS sequencing were performed as per the routine pipeline established at the facility of Fasteris, Geneva, Switzerland. Illumina HiSeq2000 was used to carry out paired-end sequencing generating a minimum of 100 bp per read.

All reads from Illumina sequencing (ArrayExpress accession: E-MTAB-4972) were first quality-assessed using FastQC (Andrews, 2010) and mapped to Arabidopsis (TAIR-10) genome using Tophat-2.0.9 (Kim et al., 2013). The percentage of reads mapped from Sex-1 library were 33.0% and from Apo-1 library 34.9%, respectively. To quantify transcript abundances and differential expression analysis, Cufflinks-2.1.1 and Cuffdiff-2.1.1 tools were used with upper quartile normalization and multi-read correction was applied to procure fragments per kilobase of transcript per million mapped reads (FPKM) values (Trapnell et al., 2010, 2012). In parallel to obtaining a list of differentially expressed genes based on corrected p-values from cuffdiff (N = 87), FPKM values from cufflinks file were used to generate the second list of differentially expressed candidate genes based on twofold changes between expression values across both genotypes. It must be noted that mapping Boechera reads to Arabidopsis genome (80% similarity) is of heterologous nature, therefore, we relaxed the filtering criteria and categorized genes showing mis-regulation by twofold change. Based on this list, we found significant number of genes (∼4100) to be differentially expressed across both genotypes. The list of differentially expressed genes was then subjected to GO enrichment analysis using BiNGO tool from Cytoscape after Benjamini and Hochberg false discovery rate (FDR) correction with p-value cut-off of ≤0.05 (Maere et al., 2005). Gene ontology and annotations were improved with data acquired from the following websites: www.arabidopsis.org; http://planttfdb.cbi.pku.edu.cn; http://ahd.cbi.pku.edu.cn, where appropriate.

In vitro osmotic stress treatment was carried out using half-strength Murashige and Skoog (½MS) agar plates infused with PEG-8000 (polyethylene glycol). Once the MS media was solidified, PEG was overlaid and allowed to equilibrate overnight. Two concentrations of PEG (400 and 550 g/L) were tested on Boechera seedlings. Prior to transferring seedlings, PEG solution was poured according to (Verslues et al., 2006). Since 550 g/L PEG almost completely inhibited the root growth, 400 g/L PEG was used for final experiments. For salt (osmotic + ionic) stress, ½MS agar plates supplemented with 50 mM NaCl were used as this concentration causes a similar effect on primary root elongation in Arabidopsis (Verslues et al., 2006). Individual seedlings (n = 20 seedlings each for mock and PEG stress treatment in four replicates, and n = 8 for salt stress treatment) grown for 11–12 days on mock ½MS agar plates in vertical position were then transferred to both stress treatment and mock plates, arranged in a manner to have the same start point for further growth of roots and grown for next 7 days. Each plate contained seedlings from all three genotypes, and the experiments were replicated at least four times to avoid plate effect. Following the seventh day of treatment, images were captured and primary root length measurements were carried out using the software ImageJ (Schneider et al., 2012). Drought treatment of flowering plants pre-grown at control conditions was performed by stopping watering for 9 days until soil in pots was almost dry and some plants showed first signs of wilting, and then resuming watering to maintain the acquired weight of the pots and prevent plants drying out.

Zebularine, a chemical analog of the cytosine nucleoside that efficiently inhibit DNA methyltransferase activity (Baubec et al., 2009) was used as a DNA methylation inhibitor. Different concentrations of zebularine (50–175 μM) were first tested for their growth effect on Boechera seedlings in vitro. For final experiments, we used 175 μM zebularine treatment, as the seedlings still retained chlorophyll but their growth was inhibited. Seedlings (n = 10 each for control and for treated) grown on ½MS agar plates for 11–12 days were transferred to ½MS mock agar plates or supplemented with 175 μM zebularine. Arrangement of seedlings in the plates was similar to that of PEG treatment. Images were recorded and samples were snap-frozen for RNA allowing 7 days of growth upon exposure to the nucleoside.

Seedlings grown on vertical ½MS agar plates, and premeiotic buds from plants grown on soil in control conditions were collected and immediately snap-frozen for RNA extraction. Developmental staging of pre-meiotic flower buds used for gene expression validation was determined based on clearing of ovule samples by DIC microscopy. The bud size and images of corresponding ovule development for each genotype are furnished in Supplementary Figure S3. RNA isolation and cDNA synthesis were performed as described previously (Johnston et al., 2008, 2010). In brief, frozen tissues were ground similarly to the RNA-seq experiment; RNA was extracted using Trizol method (Chomczynski, 1993) and treated with DNase I; reverse transcription was done with SuperScript III First-Strand Synthesis System, according to the manufacturer instructions (Thermofisher, USA). To validate RNA-seq results, qRT-PCR of selected genes with strongly increased or decreased expression in the apomict was carried out. Three biological replicates of three seedlings each in three technical replicates were used for all control seedlings; similar experimental setup was also used for qRT experiments on premeiotic buds and zebularine-treated seedling experiments. The list of candidates included genes from meiotic and stress machinery. Primer sequences are listed in Supplementary Table S5. SYBR Green assays were performed with StepOnePlus Real-Time-PCR System (Applied Biosystems, USA) and normalized as described in (Pellino et al., 2011). Calculations of relative expression levels and generation of graphs was carried out using Microsoft Excel 2010.

Pre- and post-meiotic ovules were fixed in ice-cold 9:1 ethanol:acetic acid fixative and incubated overnight at 4°C, then re-hydrated in 90, 80, and 70% ethanol series. Ethanol was replaced with clearing solution (8:2:1 chloralhydrate:glycerol:water) for 12–16 h at 4°C, and ovules were dissected out on microscopy slides with help of insulin needles in a drop of clearing solution, and mounted under cover slips. Differential interference contrast (DIC) imaging was performed on a Leica DMI6000 microscope.

Venn diagrams were made with online Venn diagram generator available at http://bioinformatics.psb.ugent.be/webtools/Venn. Images were treated in Adobe suite of software such as Illustrator or Photoshop CS for quality of representation only.

Quantitative data were analyzed using one-tailed unpaired Student’s t-test for comparisons between groups. In addition, significance of differential gene expression response to treatments was evaluated using two-way ANOVA. Significance levels used were α ≤ 0.01 (^∗∗) and α ≤ 0.05 (^∗). Calculations were carried out using Microsoft Excel 2010 and http://vassarstats.net/anova2u.html.

We conducted a transcriptome-wide differential gene expression analyses between a sexual diploid, and an apomictic triploid Boechera that originated presumably through hybridization and polyploidization in the same stress-prone habitat (Figure 1, Supplementary Figure S1). We performed comparative RNA-seq analyses between the apomictic seedlings versus the sexual seedlings, in order to uncover the molecular signatures associated with the evolutionary history that might have formed the basis of evolution and maintenance of apomixis. A little over 3000 Boechera genes were found differentially expressed between Apo-1 and Sex-1 seedlings, accounting for about 15% of the total number of corresponding Arabidopsis loci qualified for differential analyses upon heterologous mapping. From the analysis of gene ontology (GO) enrichment of these genes we noted that response to stimulus and stress categories predominated, ranging up to 23% of the total (Figure 1, Supplementary Table S1). Both biotic and abiotic GO stress categories were present; and abiotic stress categories were frequented by response to osmotic, drought and salt stress categories. Other primary categories included response to wounding, hormones, ethylene, oxidative stress, metabolic processes etc. Regardless of lists of genes upregulated or down-regulated in the apomictic seedlings, the overall annotation enrichment helped us to pinpoint the essence of stress being an important category to distinguish the vegetative transcriptome of the apomict over the corresponding sexual (see the bar chart in Figure 1; Supplementary Figures S1 and S2). Intriguingly, a total of 12 genes associated with meiosis were found differentially expressed in apomictic seedlings (Table 1). We validated our transcriptome data by a two-fold approach: (a) comparison of real-time qPCR-determined expression levels of representative candidate genes not only between Apo-1 and Sex-1, but also in another apomict, Apo-2, which is distinct from Apo-1 in terms of ploidy and evolutionary history, and (b) testing deregulation of meiotic and stress response genes also in reproductive organs at onset of (apo)meiosis.

Consistent with our RNA-seq data (Table 1, Supplementary Table S3, Figures 2A,D), we could confirm that homologs of ASYNAPTIC 1 (ASY1) and MULTI-POLAR SPINDLE 1 (MPS1) were significantly down-regulated in seedlings of Apo-1 compared to Sex-1 (Figures 2B,E) as well as at pre-meiotic stage in flower buds (Figures 2C,F, Supplementary Figure S3). Between Apo-1 and Apo-2 seedlings, ASY1 level was lower in Apo-1 than Apo-2 while MPS1 showed insignificant differences between the two apomicts. Similar trend was also observed in premeiotic flower buds. In order to subtract ploidy from reproduction mode-related effects, we compared the gene expression in groups of the diploids with triploid (2n vs. 3n) and the sexual with apomicts (Sex vs. Apo). When mean seedling expression value of both apomicts was compared to that of the sexual, MPS1 expression was found significantly decreased, and ASY1 displayed a similar downregulation trend. However, ploidy groups also exhibited significant difference in gene expression, and MPS1 expression in seedlings was tested as ploidy-independent (Figure 2E). Meiotic gene deregulation was further accompanied by upregulation of gene homolog of an abiotic (light) stress response GATA TRANSCRIPTION FACTOR 24 (GATA24) (Supplementary Table S2), otherwise annotated as ZIM-LIKE 1 (ZML1) in Arabidopsis, both in apomictic seedling tissues and flower buds of Apo-1 and Apo-2 (Figures 2G–I). However, ZML1 levels were much higher in Apo-1 compared to Apo-2 in both the stages and probably reflected ploidy-specific gene expression pattern. Similarly, a biotic stress-specific gene JASMONATE-ZIM-DOMAIN PROTEIN 9 (JAZ9) was found activated in Apo-1 only (Figures 2J,K, Supplementary Table S2). ZML1 showed both significant effect between reproductive and ploidy groups, but JAZ9 expression seemed to be rather ploidy-dependent. Together, a subset of meiotic and stress-specific genes was deregulated in the apomictic seedlings and premeiotic buds, thus confirming our RNA-seq data and a possible association with the reproductive stage. This finding let us to propose that the identified transcriptional modifications in an early stage could likely follow throughout apomictic development. Accession-specific differences between the two apomicts could also be noticed pointing at the role of ploidy and/or evolutionary history in gene expression, yet suggesting their role preferentially in the apomicts. The above validations support that our dataset of about 700 stress- and stimulus-associated genes differentially expressed in the asexual over sexual seedlings could be resourceful for connecting stress, ploidy and reproductive nodes of gene regulatory networks in development.

A careful examination of our dataset indicated that transcription factor deregulation is indeed one of the major components of the apomict/stress-specific transcriptome accounting for about 15% (Supplementary Tables S1 and S2). Whereas members of the AP2/EREBP (APETALA2/ethylene-responsive element binding proteins) transcription factor family predominated equally in both upregulated and down-regulated genes in apomictic seedlings, transcription factors belonging to the NAC [NAC domain: NAM (for NO APICAL MERISTEM), ATAF 1 and 2, and CUC2 (for CUP-SHAPED COTYLEDON 2)], bZIP (BASIC LEUCINE ZIPPER) and MYB (R2R3-type) families were enriched in the upregulated gene lists (Table 2, Supplementary Tables S2 and S4). bHLH (BASIC HELIX-LOOP-HELIX), WRKY (WRKY DNA-BINDING PROTEIN) and C2H2-type Zn-Finger transcription factor families were exclusively abundant in the down-regulated category (Supplementary Table S2). Since the stimulus, hormonal control and transcription factor regulation categories are tightly inter-connected; we examined if stress-related hormonal regulatory features could be identified as specific for Apo-1 seedlings transcriptome. We noticed that nearly one-third of differentially expressed genes (26 and 36% of down-regulated and upregulated, respectively) respond to hormones such as abscisic acid (ABA), auxin, jasmonic acid (JA), ethylene, salicylic acid (SA), gibberellic acid (GA) and cytokinin (Supplementary Table S2). It is interesting to note that genes of the auxin biosynthesis and response [e.g., YUCCA8, CYTOCHROME P450 CYP79B2, INDOLE-3-ACETIC ACID INDUCIBLE 5 (IAA5), IAA19, IAA32], cytokinin response [e.g., CYTOKININ INDEPENDENT 2 (CKI2), TYPE-A RESPONSE REGULATOR (ARR4), ARR5, Arabidopsis HISTIDINE-CONTAINING PHOSPHOTRANSFER FACTOR 5 (AHP5)] regulatory members are generally down-regulated, in contrast to those that belong to biosynthesis and response of ABA [e.g., NINE-CIS-EPOXYCAROTENOID DIOXYGENASE 3 (NCED3); RESPONSIVE TO ABA 18 (RAB18; Figures 2C,D); ABSCISIC ACID RESPONSIVE ELEMENTS-BINDING FACTOR 2 (ABF2), DRE-BINDING PROTEIN 2A (DREB2A), RD26 NAC] and the JA pathway [LIPOXYGENASE 2 (LOX2), FATTY ACID DESATURASE 8 (FAD8), ALLENE OXIDE SYNTHASE (AOS), JAZ9] were upregulated in the apomictic Boechera seedlings (Table 2, Supplementary Table S2). The antagonistic situation between different hormone responses, enrichment of specific hormone regulation categories of ABA and JA, and upregulation of several genes coding for the LATE EMBRYOGENESIS ABUNDANT (LEA)-type storage proteins [e.g., LEA30, RAB18 and LEA28 (Figures 3J–N)] suggest that the apomictic seedlings likely have better stress acclimation than the sexual seedlings, perhaps in a ploidy-dependent manner.

Among the several stress-responsive transcription factors identified to be deregulated in the apomict (discussed above), we chose to examine the NAC transcription factors further. Table 2 enlists all the NAC factors uncovered by our RNA-seq approach. Particularly six Boechera protein homologs of the Arabidopsis NACs, NAC019, NAC032, NAC055, RD026, NAP and ATAF1 (Table 2, Supplementary Table S4) that presumably might function as well during apomeiotic transition, could be mapped within a probable regulatory network in Arabidopsis, as pinpointed by the String database (Figure 1) ¹. We examined if some of these NAC homologs were indeed differentially expressed. Boechera homologs of NAC019 and NAC055 showed significant upregulation both in the apomictic seedlings and in the apomeiotic flower buds of Apo-1 (Figures 3A–F). NAC032 confirmed the trend from transcriptome data and showed upregulation in Apo-1 seedlings as well. Although these genes did not show significant upregulation in the seedlings of Apo-2, the corresponding apomeiotic floral buds of Apo-2 showed stronger enrichment than Sex-1 and Apo-1 (Figures 3C,F) exhibiting a significant expression increase of the average between apomicts. The expression changes for NAC019, NAC055 and NAC032 between Apo-1 and Apo-2 were also found to be significantly different. Although there were certain differences in NAC expression profiles between the apomicts at early developmental stages, which could concern rather difference in ploidy level or genetic background, NAC055 and NAC019 expression in flower buds was not different between ploidy groups but only significantly reproductive-mode-dependent. Together, NAC gene family upregulation may represent a key unifying background for early stress response, in apomictic Boechera species.

To understand stress sensitivity of Boechera species, we carried out in vitro osmotic stress treatments in line with the output from transcriptome analysis of Boechera seedlings featuring stress signatures such as drought, salinity and osmotic stress tolerance. We treated Boechera seedlings with PEG-8000 (polyethylene glycol 8000), using plate assays, generating changes in osmotic balance and therefore decreasing water availability as a major component of drought stress (Verslues et al., 2006). We tested multiple concentrations of PEG for Boechera growth and found that PEG 400 g/L would impose moderate osmotic stress; PEG 550 g/L seemed to be too harsh for the seedlings as most of the exposed plants did not survive well, regardless of accessions (data not shown). Alteration in osmotic conditions had a noticeable root growth effect in Sex-1 roots but not in the apomicts. Primary root elongation of the Sex-1 seedlings treated with PEG 400g/L was found to be smaller by ∼25% in comparison to mock treated Sex-1 seedlings (Figures 4A–C). However, for both Apo-1 and Apo-2 seedlings, no significant difference was observed in primary root lengths in either of the treatments, suggesting the prevalence of a better osmotic stress tolerance mechanism in both apomicts than in the sexual. This finding was further corroborated by the differential response of the sexual and apomicts to salt concentrations that cause both osmotic and ionic stress using 50 mM NaCl which caused a very similar effect decreasing the primary root growth exclusively in the sexual (Figure 4C). Whereas primary root extension was similar in the control treatment across all Boechera accessions, only the sexual exhibited stress-induced effect, and both osmotic and salt stress caused significantly different effect between reproductive groups but that was not ploidy-dependent. In addition to the in vitro assays, we enforced moderate drought on soil grown plants in standard growth room conditions at 22°C to understand if osmotic stress induced by drought impact the overall growth and reproduction of the three genotypes. We analyzed the ratio of reproductive modes in ovules upon drought stress, asking if stress impacted meiosis versus apomeiosis. The proportion of rare sexual-like triads and the routine meiotic tetrads in relation to apomeiotic dyads did not detectably change upon drought exposure in both apomicts (Figure 4). In sum, apomictic seedlings seem to be better acclimated to osmotic and ionic stress in seedlings, but stress does not significantly shift obligate apomeiotic mode to sexual mode during the reproductive phase.

Our transcriptome approach uncovered the parallel alterations between stress and meiosis regulation in the apomictic Boechera. Unpublished data from our lab and meta-analysis of our transcriptome data to large scale databases indicated definite role of methylation machinery in apomixis specific gene expression changes. We indeed observed that the Arabidopsis homologs of several of our stress-specific candidate transcription factors and meiotic genes are likely epigenetically regulated in seedlings via DNA methylation and/or small RNAs, or repressive histone methylation marks (Supplementary Tables S2–S4). In the absence of full genome and methylome sequences of all three species studied here, we conducted a DNA methylation inhibitor assay using a demethylation agent zebularine (Baubec et al., 2009), to identify gene expression changes associated with global DNA methylation. Zebularine exposure resulted in significant effects on all genotypes causing strong decrease in root length (Figure 5A, compare to Figure 4A), due to the genome-wide demethylation effect of this nucleoside combined with its redox-active nature, similar to growth retardation of Arabidopsis when applied in vitro. It must be noted, however, that there was small but significant difference in response upon zebularine treatment; both apomicts retained a slightly better root elongation than the sexual (Figure 5B). We examined the effect of this DNA methylation inhibitor on selected candidate genes by real-time qRT-PCR comparing zebularine-treated versus mock-treated seedling gene expression of meiotic genes: ASY1 and MPS1; and a subset of stress-responsive transcription factors: ZML1, NAC019, NAC032 and NAC055. The first four candidates seem to be regulated by DNA methylation as their expression was affected by zebularine exposure (Figures 5C–E), and both NAC032 and NAC055 did not show similar regulation (data not shown). The qRT-PCR results showed some interesting changes in expression of the candidate genes. MPS1 seems to be repressed by DNA methylation in all species regardless of the reproductive mode, as there was no significant difference in response to zebularine application between the Boechera accessions used (Figure 5B). ASY1, ZML1 and NAC019 expression upon zebularine treatment was either unaltered or down-regulated in the sexual seedlings but they were strongly activated in Apo-1, and ASY1 upregulation seem to be common in both apomicts and different from the sexual (Figures 5C–F). ASY1 and NAC019 expression responses to zebularine were significantly different between species based on two-way ANOVA test, while ZML1 exhibited rather ploidy-dependent effect. Our data suggest locus-specific gene expression changes upon induced global DNA demethylation coupled with cellular stress distinguish between both ploidy groups and apomictic and sexual genotypes.

Boechera is one of the primary plant model systems to study evolution and regulation of apomixis associated traits such as meiotic circumvention (apomeiosis), apospory, parthenogenesis, polyploidy, hybridity and environment (Roy, 1995; Schranz et al., 2005; Rushworth et al., 2011). Boechera genus in particular comprises of numerous independently emerged asexual (sub)species indicating that this genus may harbor a specifically pronounced (epi)genetic background favorable for evolution of apomixis. Along with polyploid apomicts, Boechera genus also includes diploid apomictic populations that were presumably sired from apomictic polyploids or aneuploids, or occurred via hybridization (Kantama et al., 2007; Lovell et al., 2013). Apomixis in diploid condition is extremely rare, and polyploidy and hybridity seems to have been the springboard for evolution and maintenance of apomixis (Comai, 2005; Hand and Koltunow, 2014; Hojsgaard et al., 2014b). Unlike the polyploid apomictic parents, the rare apomict diploids might display segregation of apomictic traits (Noyes and Wagner, 2014) and fitness trade-offs, hence apomictic diploids do not seem to survive in most apomictic genera. In triploids, where meiotic synapsis of chromosomes fails causing so-called “triploid block,” apomixis remains the only possible mean of reproduction through seeds; to some extent this concerns other polyploids too. However, as most diploid apomicts in Boechera seem to represent products of hybridization (Beck et al., 2012), a similar “genomic block” for meiotic bivalent formation may take place. Yet another route for evolution of a triploid apomict might have been endo-duplication due to triggered stress response following exposure to harsh environmental conditions (Scholes and Paige, 2015). This scenario might have led to diploid (sexual) gamete in one genotype, which was then fertilized by a haploid sexual gamete, thus resulting in triploid genotype. We reported similar endo-duplication events in an Arabidopsis mutant (Johnston et al., 2010), and we believe it was due to genomic stress in the female gamete mediated by loss of a crucial cell cycle regulator. Especially triploids could have paved way for evolution of apomixis by fixing and maintaining heterozygosity. Spontaneous apospory reported in Ranunculus triploids (Hojsgaard et al., 2014a) might lend additional support along this line. In addition, asynchronous gene expression from different genomes is thought to facilitate asexual re-programming in hybrids (Carman, 1997; Hojsgaard et al., 2014a). As hybridity leads to “genomic shock” thus “genomic stress,” we may not be able to tease apart cause and consequence.

Understanding apomixis traits has been mostly restricted to transcriptome analyses of apomixis-specific organs or cell-types, and population genetic studies. For instance, Sharbel et al. (2010) reported heterochronic changes in gene expression during initiation of apomixis along with decrease of their transcript levels in the apomict by comparing expression tags across four stages of ovule development between diploid apomictic and diploid sexual Boechera species. Another study analyzed apomictic cell type-specific (e.g., MMC) gene expression profiles of a triploid apomictic Boechera similar to the one we used in this work, and compared to the corresponding profiles of a distantly related (10 million years of divergence) sexual species A. thaliana (Schmidt et al., 2014). Although number of interesting gene expression changes related to meiotic and mitotic cell cycle regulation, hormonal pathways, signal transduction, stress and epigenetic regulatory pathways were alluded to be specific for the apomeiotic MMC, comparative analyses seem to be constrained due to the lack of appropriate sexual control. It was also noted that all major core meiotic genes were expressed in the apomictic MMC in this triploid apomict, and the authors proposed that diplosporous apomeiosis likely result from modification of the core meiotic and cell division pathway in the MMC. Although some candidate loci underlying apomixis have been isolated or reported, functional analyses of apomixis in Boechera are limited primarily owing to the inherent difficulties in plant transformation of these species.

Environment has always remained a crucial factor with the highest impact on homeostasis and plant life cycle. Whereas most accessions of the model plant A. thaliana are drought-sensitive, the Boechera clade in general and a number of Boechera species in particular seem to have acquired an asset of biotic and abiotic stress-tolerance mechanisms, perhaps due to their predominant ability to frequently hybridize, forming novel species that occupy distinct ecological niches primarily in their native habitat ranges; and due to their repeated exposure to extreme and fluctuating environmental conditions involving major abiotic stress factors like drought, cold and several biotic factors as depicted in Figure 6 (Rushworth et al., 2011). Apomixis is absent in the Arabidopsis clade but prevalent in the Boechera; therefore, it is possible that general stress acclimation could also be an inherent characteristic background for evolution of apomicts, as it might likely provide an efficient (epi)genomic constitution to efficiently tolerate environmental hazards and to buffer mutation loads in the absence of recombination. In parallel, several stress-response genes have been deregulated in transcriptomes of apomictic organs and cell-types of Boechera (our unpublished data; Sharbel et al., 2010; Schmidt et al., 2014). In sum, Boechera represents an excellent model system for identifying molecular cues to gain in-depth understanding of apomixis regulation and its correlation with environmental change and stress response in a diverse genetic background linked to hybridity and polyploidy.

Apomictic Boechera in general seem to be more stress tolerant, at least in terms of osmotic stress acclimation, than the corresponding sexual. The apomicts grown in normal conditions exhibit elevated levels of a considerable number of stress-response genes when compared to the sexual (Supplementary Tables S1 and S4, Table 2, Figure 3). These genes include prominent drought and ABA response factors such as transcription factors NACs, bZIPs, DREB2A, MYBs, and LEA proteins (Nuruzzaman et al., 2013). For example, proper function of bZIP transcription factors AREB1, AREB2 and AREB3 has been shown to be important for drought tolerance (Amir Hossain et al., 2010; Golldack et al., 2011). When overexpressed, Arabidopsis NAC019, NAC055, and NAC072 elevated drought tolerance (Tran et al., 2004).

We noticed though that expression levels of some genes tested were significantly different between the two apomictic accessions used in this study. Both Apo-1 and Apo-2 are obligate diplosporous apomicts, and the main difference between them is their evolutionary history including ploidy level (Apo-1 being triploid and Apo-2 diploid), ecological habitats with different abiotic stress severity, distant geographical location and, therefore, almost certainly distinct parental genotypes. In order to disentangle the ploidy effect on the gene expression, in the absence of available diploid versus polyploid Boechera transcriptome data, we compared our results with published datasets of the Arabidopsis clade. About 4% of the Apo-1-specific upregulated genes showed the same behavior in Arabidopsis allo-tetraploids in an inter-specific hybridization study (Wang et al., 2004), and for stress-related genes it was only 1.8%. Here, three genes NAC019, NAC055 and ZML1 from our validated list were upregulated. Similar changes are not discernible in Arabidopsis auto-tetraploids compared to the corresponding diploid transcriptomes (Pignatta et al., 2010). Therefore, ploidy might not have been the sole cause of genotype-specific transcription differences and shaped up the entire transcriptome of the Apo-1. As Boechera apomicts were presumably sexual hybrids by origin that were eventually fixed by apomixis, hybridity between different genomes, thus genome collisions could have played a role for activation of stress-response genes and initiation of apomixis, as evident from upregulation of genes like NAC019 and NAC055 homologs in both apomictic floral buds. Upregulation of the homologs of stress-related genes in drought-tolerant apomictic Boechera correlates with their crucial function in stress adaptation shown in other plants and support the hypothesis that this constitutive overexpression confers better stress amelioration, in particular osmotic stress tolerance. Our transcriptome data show an intricate network between hormonal regulation and stress response; for instance, preferential activation of several ABA and JA-specific biosynthesis genes and responsive transcription factors in Apo-1 pinpoints that the apomict has an ability to withstand multiple stresses including osmotic stress and drought (Nguyen et al., 2016). Activation of mobile mRNAs of ascorbic peroxidases (APX1 and APX2) that are pivotal for maintaining cellular redox homeostasis by scavenging hydrogen peroxide and ROS, consequent down-regulation of genes of the auxin biosynthesis and response genes, giving thus a possibility to perturb auxin homeostasis, collectively suggest that the apomict Apo-1 has an efficient auxin-dependent growth control in order to cope-up with stress in comparison to Sex-1 (Tognetti et al., 2012; Naser and Shani, 2016), and in turn might activate ABA-specific responses required for stress acclimation. Intriguingly, whereas specific auxin-response factors (ARFs) seem to be down-regulated in seedlings (this study) and many ARFs down-regulated in ovules (Shah et al., unpublished data), they were reported to be enriched in the MMC cells of Apo-1 (Schmidt et al., 2014), perhaps due to the underlying hormonal homeostasis shifts between the seedling, nucellus stages, and the MMC prior apomeiosis. The latter data, however, is based on comparison between MMC cells of Arabidopsis versus Boechera, therefore, further validation would be required. The expression of most genes tested showed significant difference between the sexual and apomicts, with ploidy effect being significant at many instances; yet the effect of stress on in vitro growth was clearly dependent on the reproductive mode of the Boechera species but not the ploidy level. In sum, whether it is Apo-1, or the evolutionarily distinct Apo-2, there is a clear trend that the Boechera apomictic seedlings confer higher level of stress acclimation and underlying gene response in seedlings than the sexual counterpart, and likely play a pivotal role in sustaining apomixis through to apomeiotic transition.

Apomeiotic switch has been perceived as a developmental reprogramming event connected to meiotic shunt, stress, epigenetic and transcriptional changes, transposon control etc. (reviewed in Grimanelli, 2012; Horandl and Hadacek, 2013; Mercier et al., 2015). Previous work in Arabidopsis has indicated the importance of meiotic mutations in producing unreduced gametes by apomeiosis-like program. For instance, DYAD/SWITCH1 (SWI1) is an important candidate gene in Arabidopsis reported for switching the meiosis program toward apomeiosis (Ravi et al., 2008). In MiMe-1 and MiMe-2 triple mutants, a diplospory-like division also leads to the formation of unreduced gametes. By combining different meiotic mutant combinations [(MiMe-1: sporulation11-1 (spo11-1), omission of second division1 (osd1) and recombination8 (rec8); MiMe-2: spo11-1, osd1 and cyc1;2/tardy asynchronous meiosis (tam)], diplospory-like meiotic restitution have been successfully recreated in Arabidopsis (d’Erfurth et al., 2009, 2010). In maize, a mutation in an effector of small RNA pathway has been shown to be important for recreating diplospory-like event (Singh et al., 2011). Whereas role of these genes or mutations underlying apomixis in a natural apomict like Boechera has not yet been documented, Arabidopsis versus Boechera (Apo-1) MMC comparisons suggested possible down-regulation of certain meiotic genes: PARTING DANCERS, SPO11-2, TAM, SWI1, SOLO DANCERS and RAD50 (Schmidt et al., 2014). Whereas this data need real-time validation, none of these genes were deregulated in our transcriptome of Apo-1 versus Sex-1 seedlings, but deregulation of two other meiotic genes were noteworthy. Acute down-regulation of core meiotic homologs such as ASY1 and MPS1 both in seedlings and apomeiotic buds gives a hint that a global alteration in expression of these genes are likely necessary in order to make the apomeiotic switch. ASY1 is a Horma-domain containing DNA-binding protein component of the meiotic axis involved in synapsis and crossover formation, and it is required for DMC1-mediated interhomolog recombination during female and male meiosis (Armstrong and Jones, 2001; Sanchez-Moran et al., 2007). MPS1 is a coiled-coil protein involved in meiotic spindle formation and it is required for faithful chromosome segregation during female and male meiosis (Jiang et al., 2009). Interestingly, none of the meiotic genes we identified in our study (Table 1, Supplementary Table S3) showed transcriptional changes in transcriptomes of Arabidopsis polyploids or hybrids (Wang et al., 2004; Yu et al., 2010). Could evolutionary history traits in Boechera different from that of Arabidopsis influence expression of these genes? It is also likely that this transcriptional homeostasis reflects a genome-wide change in apomictic Boechera and might involve heritable epigenetic silencing.

Abiotic stress has been shown to affect meiotic progression, whereas different type, severity and duration of the stress factors may have different effects on reproduction (reviewed in De Storme and Geelen, 2014). In addition, how stress affects meiosis could be different between male and female tissues. As female tissues are more somatically dependent than the male organs, for instance, female meiosis seems to progress during water deficit but male meiosis abort in Arabidopsis (Saini, 1997). In barley, heat stress leads to precocious onset of meiotic chromosome condensation during male meiosis and a concomitant ASY1 upregulation (Oshino et al., 2007). In Arabidopsis seedlings, ASY1 seems to be marginally down-regulated during drought, and ASY1 expression levels seem to fluctuate inconsistently between time points during prolonged stress condition (Kilian et al., 2007). The same is true for MPS1 upon light stress, and for DYAD during cold, drought, and heat; therefore, a link between stress response and meiotic machinery may exist much earlier than onset of reproduction itself. However, direct connections between stress response and regulation of meiotic genes in apomicts have not been elaborated earlier. The stress-co-adapted and highly heterozygous apomictic genomes, similarly to polyploids (Chandra and Dubey, 2010; Coate et al., 2013; Hao et al., 2013; Syamaladevi et al., 2016), generally are considered to be a haven for buffering mutation load in the absence of meiotic recombination (Lynch, 1984). Therefore, not only those apomictic plants are well prepared for stress, they also feature “apomeiotic” mutations or down-regulation of genes such as MPS1 and ASY1 right from the very growth of the apomictic plant. The broadly adapted apomicts may more efficiently buffer mild abiotic stress than the sexuals, thus limiting cellular oxidative stress and formation of ROS and resulting DNA damage, hence, there is no need for meiosis and recombination in the apomict [sensu. FAR theory (Hadany and Beker, 2003a,b) and oxidative damage initiation hypothesis or meiosis (Bernstein and Bernstein, 2013; Horandl and Hadacek, 2013)]. In agreement with this, in our experimental setting, we observed neither reproducible transcriptional changes of Boechera homologs of MPS1 and ASY1 in seedlings upon osmotic stress (data not shown) nor reversion of apomeiosis to sexual-like meiosis. Recent preliminary reports indicate that a facultative apomictic Boechera species might incline toward sexual meiosis in a substantial fraction of ovules under drought and heat (Carman et al., 2015). These results are also reminiscent of similar observations in a facultative apomict Ranunculus under extended photoperiod (Klatt et al., 2016), and collectively hint that residual sexuality contained in the facultative apomicts could be triggered upon environmental cues. In light of our observation that obligate apomictic Boechera accessions exposed to mild drought did not exhibit any significant reversion to sexuality in our conditions, we propose that additional stress may not impinge on the reproductive mode of obligate apomicts that were already co-adapted for stress acclimation and meiotic circumvention.

Among several transcription factor groups that are predominantly transcriptionally activated in the stress-tolerant apomict, we have noticed upregulation of the family of NAC transcription factors, which have previously been reported to work in tandem with each other (Hickman et al., 2013). This reproduction mode-specific activation extends beyond seedling to reproductive phase, proposing the prime importance of this regulatory network during apomixis. Stress-response and developmental role of NACs have previously been discussed in various plant systems (Meng et al., 2009; Nuruzzaman et al., 2013; You et al., 2015). It is interesting to note that NAC genes chosen as candidates from our RNA-seq (NAC019, NAC032, NAC055) are also known to be expressed in specific cell-types and ovules in Arabidopsis (Yu et al., 2005; Johnston et al., 2007; Wuest et al., 2010; Schmidt et al., 2011; Schmid et al., 2012). Although the micro-dissected RNA-seq did not give consistent presence calls for several NACs in laser-captured single cell transcriptomes presumably due to experimental noise, NAC032, NAC002 and NAC103 homologs of Boechera were consistently present in Apo-1 apomeiotic cell types such as the MMC and surrounding ovular cells (Schmidt et al., 2014) (Supplementary Table S4). A recent work identified novel role of NAC family members in ovule separation (Goncalves et al., 2015), thus in addition to stress regulation, NACs might govern early reproductive development. A number of NACs are also known for their non-cell-autonomous movement via movement of their mobile mRNA (Thieme et al., 2015); therefore, one possibility may be that they maternally control decisions on sexual or apomictic reproductive development.

Epigenetic regulation involving DNA methylation and histone modifications has long been recognized as a potent regulator of stress-response in plants (reviewed in Probst and Mittelsten Scheid, 2015; Bilichak and Kovalchuk, 2016). Whereas the underlying regulation is extremely complex, a two-step model proposed for transcriptional activation of stress-response genes includes first DNA methylation in CG context, which then directs H3K9 methylation recruiting non-CG methylation (reviewed in Du et al., 2015; Probst and Mittelsten Scheid, 2015). Whereas DNA hypermethylation of stress-inducible genes might lead to enrichment of H3K9me2 and a concomitant depletion of histone acetylation, as shown for Arabidopsis seedlings under saline stress (Bilichak et al., 2012), drought activated genes are associated with elevation of H3K4 trimethylation and H3K9 acetylation (Kim et al., 2008). However, epigenetic specificity between DNA and histone methylation under stress might vary between genes, as shown for salt-responsive transcription factors in soybean (Song et al., 2012). A review of the previous data showing locus-specific DNA and histone methylation marks, and small RNAs in Arabidopsis (Supplementary Table S4) suggests that the corresponding Boechera loci might be under complex epigenetic control. Specifically, NACs, MYBs, and several other genes coding for stress-related transcription factors in different plant species have been shown to be targets of miRNAs (Rhoades et al., 2002; Guo et al., 2005), a class of small RNAs involved in post-transcriptional gene silencing and, in many of these cases, induced by stress. In light of these previous findings, our data of Boechera NAC019 being significantly deregulated in response to global demethylation, supports an underlying epigenetic regulation in the apomicts. The Arabidopsis NAC019 was previously shown to be activated upon allo-polyploidization in large-scale transcriptomic datasets (Wang et al., 2004), supporting the hybridization effect also in Boechera. As Boechera apomicts studied here likely represent allopolyploid and diploid hybrids, it is important to consider the importance of epigenetic landscape of the entire genome to be successful apomicts. Global changes in DNA methylation, RdDM-dependent transcriptional changes etc. would be important for the hybrid/polyploid genome to perpetuate apomixis in plant populations (reviewed in Verhoeven and Preite, 2014). Besides, the heterozygous genome and epigenetic landscape instigated by hybridity or polyploidy is compatible with the efficient stress ameliorative advantages in the apomicts, as discussed above.

In addition to the previous data that apomixis is largely epigenetically controlled via pathways of DNA methylation in several apomictic plant systems (Garcia-Aguilar et al., 2010; Zappacosta et al., 2014; reviewed in Hand and Koltunow, 2014), and also by observing the fact that DNA methylation pathway members are deregulated in Boechera apomictic cell types (Schmidt et al., 2014) and generally at different developmental stages of Boechera apomicts (unpublished data), our data indicating control of MPS1 and ASY1 by DNA methylation highlights the importance of an epigenetic angle underlying meiotic deregulation. In this work, we have provided indirect evidence that global DNA methylation might possibly control the apomict-specific transcriptional regulation of stress and meiotic genes in the seedlings in a ploidy-independent manner, though this effect could be coupled to a rather marginal DNA-damaging effect of zebularine (Liu et al., 2015). It is important to note that correlation of DNA methylation status with expression is complex and modifications specific to each cytosine will have to be mapped in the future by locus-specific bisulfite sequencing. Nonetheless, our work provides a base to further investigate the connection between genome-wide and locus-specific DNA methylation and gene expression alterations in Boechera species and its possible correlation with the apomixis program. It will also be interesting to examine how epigenetics, the redox status and DNA breaks play a role in regulation of meiotic genes like ASY1. Are there connection between apomict-dependent ASY1 regulation, and SPO11 (Horandl and Hadacek, 2013), which are critical for double strand break repair and chromosome pairing, and may their expression be induced by DSB-dependent signaling network at different developmental stages? The decreased levels of Boechera MPS1 and ASY1 expression starting from seedling stage along with constitutive stress response in the apomicts allow us to propose that apomictic initiation is pre-programmed from early plant development onward, perhaps due to global (epi)genome-wide deregulation and resulting changes in the transcriptome. Taken together, our work forms one of the first glimpses to view apomixis as a possible read-out of stress-related genome- and transcriptome-wide alterations. Whereas previous studies in Boechera have given importance to gene expression changes that occurred during the apomeiotic shift or parthenogenetic reprogramming during the reproductive phase, our work for the first time examines the earliest molecular events during the vegetative phase of the apomicts, and provides a transcriptomic resource for targeted work. The dynamic changes connecting these factors in evolutionary timeline have possibly resulted in genome-wide transcriptome adjustments uncovered in Boechera. The work has hinted the importance of the evolutionary history, constitutive stress response elevation, meiotic gene deregulation for priming apomixis in the apomictic species, and we anticipate that this phenomenon will likely be a common background for most, if not all apomicts, similar to the proposals in the apomeiotic dandelion (Verhoeven and van Gurp, 2012). Not only in plants, stress and/or epigenetics are important denominators that impinge on asexual reproductive mode in animals too, for instance in insects and lizards (Kearney et al., 2006; Jedlicka et al., 2015) and in mice (Kono et al., 2004). Nature must have dictated co-evolution of these mechanisms across kingdoms.