Pistil Transcriptome Analysis to Disclose Genes and Gene Products Related to Aposporous Apomixis in Hypericum perforatum L.

Unlike sexual reproduction, apomixis encompasses a number of reproductive strategies, which permit maternal genome inheritance without genetic recombination and syngamy. The key biological features of apomixis are the circumvention of meiosis (i.e., apomeiosis), the differentiation of unreduced embryo sacs and egg cells, and their autonomous development in functional embryos through parthenogenesis, and the formation of viable endosperm either via fertilization-independent means or following fertilization with a sperm cell. Despite the importance of apomixis for breeding of crop plants and although much research has been conducted to study this process, the genetic control of apomixis is still not well understood. Hypericum perforatum is becoming an attractive model system for the study of aposporous apomixis. Here we report results from a global gene expression analysis of H. perforatum pistils collected from sexual and aposporous plant accessions for the purpose of identifying genes, biological processes and molecular functions associated with the aposporous apomixis pathway. Across two developmental stages corresponding to the expression of aposporous apomeiosis and parthenogenesis in ovules, a total of 224 and 973 unigenes were found to be significantly up- and down-regulated with a fold change ≥ 2 in at least one comparison, respectively. Differentially expressed genes were enriched for multiple gene ontology (GO) terms, including cell cycle, DNA metabolic process, and single-organism cellular process. For molecular functions, the highest scores were recorded for GO terms associated with DNA binding, DNA (cytosine-5-)-methyltransferase activity and heterocyclic compound binding. As deregulation of single components of the sexual developmental pathway is believed to be a trigger of the apomictic reproductive program, all genes involved in sporogenesis, gametogenesis and response to hormonal stimuli were analyzed in great detail. Overall, our data suggest that phenotypic expression of apospory is concomitant with the modulation of key genes involved in the sexual reproductive pathway. Furthermore, based on gene annotation and co-expression, we underline a putative role of hormones and key actors playing in the RNA-directed DNA methylation pathway in regulating the developmental changes occurring during aposporous apomixis in H. perforatum.

Unlike sexual reproduction, apomixis encompasses a number of reproductive strategies, which permit maternal genome inheritance without genetic recombination and syngamy.
The key biological features of apomixis are the circumvention of meiosis (i.e., apomeiosis), the differentiation of unreduced embryo sacs and egg cells, and their autonomous development in functional embryos through parthenogenesis, and the formation of viable endosperm either via fertilization-independent means or following fertilization with a sperm cell. Despite the importance of apomixis for breeding of crop plants and although much research has been conducted to study this process, the genetic control of apomixis is still not well understood. Hypericum perforatum is becoming an attractive model system for the study of aposporous apomixis. Here we report results from a global gene expression analysis of H. perforatum pistils collected from sexual and aposporous plant accessions for the purpose of identifying genes, biological processes and molecular functions associated with the aposporous apomixis pathway. Across two developmental stages corresponding to the expression of aposporous apomeiosis and parthenogenesis in ovules, a total of 224 and 973 unigenes were found to be significantly up-and down-regulated with a fold change ≥ 2 in at least one comparison, respectively. Differentially expressed genes were enriched for multiple gene ontology (GO) terms, including cell cycle, DNA metabolic process, and single-organism cellular process. For molecular functions, the highest scores were recorded for GO terms associated with DNA binding, DNA (cytosine-5-)-methyltransferase activity and heterocyclic compound binding. As deregulation of single components of the sexual developmental pathway is believed to be a trigger of the apomictic reproductive program, all genes involved in sporogenesis, gametogenesis and response to hormonal stimuli were analyzed in great detail. Overall, our data suggest that phenotypic expression of apospory is concomitant with the modulation of key genes involved in the sexual reproductive pathway. Furthermore, based on gene annotation and co-expression, we underline a putative role of hormones and key actors playing in the RNA-directed DNA methylation pathway in regulating the developmental changes occurring during aposporous apomixis in H. perforatum.
Keywords: Hypericum perforatum, sexual reproduction, aposporous apomixis, microarray INTRODUCTION Apomixis defines a number of reproductive strategies, which, unlike sexual reproduction, permit the inheritance of the maternal genome without genetic recombination and syngamy. This strategy of asexual reproduction is documented in more than 120 angiosperm genera (Carman, 1997); however, no major seed crop species are apomictic, and attempts to introduce the apomixis trait into crop plants from apomictic relatives via conventional breeding schemes have been largely unsuccessful (Barcaccia and Albertini, 2013). This asexual mode of seed formation is a trait with enormous economic and social potential in agriculture, as genetically fixing highly complex genotypes, including hybrid cultivars, through apomixis would have tremendous advantages in plant breeding, biomass and seed production (Calzada et al., 1996).
The key biological features of apomixis are the circumvention of meiosis (i.e., apomeiosis), the differentiation of unreduced embryo sacs and egg cells, and their autonomous development in functional embryos without fertilization (i.e., parthenogenesis), and the formation of viable endosperm either via fertilizationindependent means or following fertilization with a sperm cell (Koltunow and Grossniklaus, 2003). In particular, gametophytic apomixis occurs with the parthenogenic development of unreduced egg cells from apomeiotic embryo sacs that can originate from a nucellar somatic cell (i.e., apospory) or a megaspore mother cell with no, or modified, meiosis (i.e., diplospory). These features deviate from sexuality as the capability to develop an embryo sac is strictly restricted to the reduced functional megaspore deriving from meiosis, and failure of the meiotic programme naturally results in the abortion of the ovule. While polyploidy does not appear to be absolutely required for the expression of apomixis, plant species reproducing through gametophytic apomixis are very often polyploids and/or hybrids in nature (Bicknell and Koltunow, 2004). Hence, it has been suggested that apomixis might be a consequence of epigenomic shock resulting from interspecific hybridization and polyploidization, which induce transcriptome changes and the deregulation of reproductive development (Garcia-Aguilar et al., 2010). Several lines of evidence suggest that such dramatic shifts in the reproductive process could rely on spatial and/or temporal changes in the expression of sexual pathway-related genes (Grimanelli et al., 2003;Sharbel et al., 2010).
Plants reproducing via aposporous apomixis (Nogler, 1984) possess the ability to differentiate viable unreduced embryo sacs from somatic cells of the ovule (i.e., aposporous initial cells). In principle, these cells are somatic cells belonging to the nucellus, which change their fate by being able to mitotically divide and develop functional embryo sacs by mimicking sexual gametogenesis development (Barcaccia et al., 2007;Galla et al., 2011). Whether the competence to develop functional embryo sacs is retained by all nucellar cells or is restricted to a few nucellar cells physically connected to the meiocytes or the meiospores in the ovule is not yet fully understood.
In recent years, expression studies and the characterization of mutants that are defective in their ability to set seeds has led to the identification of many genes whose expression patterns are necessary for proper ovule development and correct progression of sporogenesis and gametogenesis in model plant species (Christensen et al., 1997;Drews et al., 1998;Drews and Koltunow, 2011). The ability to undergo meiosis and gametogenesis is apparently restricted to particular cell types and is determined by the expression and interaction of mRNAs and small RNAs which are necessary for the proper specification of cell identity within the ovule as well, as for the onset and progression of sporogenesis and gametogenesis (Olmedo-Monfil et al., 2010;Armenta-Medina et al., 2011;Singh et al., 2011;Tucker et al., 2012). Interestingly, Garcia-Aguilar et al. (2010) reported that the DNA methylation pathway active during reproduction is essential for gametophyte development, and observed that lossof-function mutants for two genes involved in the maintenance of DNA-methylation status produced phenotypes that are strikingly reminiscent of apomictic development, suggesting that DNA methylation in the maize ovule might regulate the transcriptional expression of genes involved in the differentiation between apomixis and sexual reproduction.
Cellular levels of hormones including cytokinins, auxins and polyamines, among others, are also very important for fine-tuning ovule and embryo sac development (Bencivenga et al., 2012;Ceccato et al., 2013;Schmidt et al., 2014). Thus, while cytokinin negatively regulates cell proliferation in the sporophytic tissues surrounding the developing embryo sac (Cheng et al., 2013), auxin biosynthetic and influx genes are necessary for correct growth of the female gametophyte and proper specification of its cell types (Panoli et al., 2015).
Hypericum perforatum is an invasive perennial herb that is widely distributed in a variety of habitats and is regarded as a serious weed in many countries (Robson, 2002;Nürk et al., 2013). In recent years, H. perforatum has been studied to identify potential genes involved in the biosynthesis of active metabolites (He et al., 2012;Hofrichter et al., 2013). H. perforatum reproduces via aposporous apomixis (Noack, 1939;Davis, 1967;Mártonfi et al., 1996;Matzk et al., 2001;Barcaccia et al., 2007;Galla et al., 2011Galla et al., , 2015. As with other asexual plant complexes, apomixis and hybridization are closely linked in H. perforatum (Koch et al., 2013), and interestingly, the dosage of genetic factors has been proposed to influence the penetrance of apomixis, as tetraploid and hexaploid genotypes tend to be more apomictic and sexual, respectively, regardless of geographic origin (Molins et al., 2014) Recently, next generation sequencing (NGS) technologies have been used to sequence the flower transcriptomes of obligate sexual and unrelated apomictic H. perforatum genotypes (Galla et al., 2015). Here we expand on this work by performing a gene expression analysis of H. perforatum pistils collected from sexual and aposporous genotypes for the purpose of identifying genes and gene products, along with biological processes and biosynthetic pathways/networks potentially associated with aposporous apomixis.

Plant Materials
Hypericum perforatum L. plants from 2 naturally occurring tetraploids (2n = 4x = 32) and 3 induced tetraploids (2n = 4x = 32) were used for the microarray approaches (Table 1). For the production of induced tetraploids, plants of the diploid sexual line R1 (reselected from the tetraploid apomictic cultivar "Topaz") were converted to auto-polyploids by colchicine application, as reported by Schallau et al. (2010). Briefly, seeds were imbibed in water for 24 h, placed on filter paper soaked with an aqueous solution of 0.2% colchicine for 24 h and then planted in soil for germination. The C0 plants that survived this treatment were self-pollinated, and their progenies (C1) were screened for tetraploid plants. Sexual tetraploid C1 and C2 plants were used for further crosses with tetraploid apomictic pollinators and were screened for their reproductive phenotype (Schallau et al., 2010;Molins et al., 2014). The reproductive mode of all H. perforatum accessions was estimated by flow cytometric screening of 48 single seeds as described by Matzk et al. (2001).
Pistils were collected from three apomictic and three sexual plant accessions. In both cases, the developmental stages considered for sampling were bud lengths of approximately 4.0 mm, corresponding to Arabidopsis flower stage 11, and bud lengths of approximately 10.0 mm, corresponding to Arabidopsis flower developmental stage 14 (Galla et al., 2011).

Array Design, cDNA Synthesis and Hybridization
Sequences used for the design of 4x180 arrays (Agilent) were previously deposited at DDBJ/EMBL/GenBank under the accession GBXG00000000. The 60,594 tentative consensus sequences were used to design 60-mer oligonucleotide probes for a microarray analysis. Sequences corresponding to transcripts <100 nt in size were removed to reduce the set of target transcripts to 55,682 sequences. The custom microarray with the design ID 065680 was created using eArray, a free Agilent webbased application that enables the creation of custom microarray designs and oligo libraries (https://earray.chem.agilent.com/ earray/). For probe design, the Base Composition Methodology option was employed, and no linker sequence was used. Two different probes were designed for each sequence. A total of 111,200 probes were obtained from eArray and were used for creation of the 4x180K array format. Slides with printed arrays were ordered directly from Agilent Technologies (Santa Clara, CA, USA).

Statistical Analysis
Scanned images were transformed into quantified figures using the Agilent Features extraction software (v11.5), and expression data were normalized based on the 75th percentile. Principal component analyses were carried out using T-MeV v4.9.0 software (http://mev.tm4.org) and by using Manhattan distances. Statistical analyses (t-test, unpaired; p < 0.01) were carried out using T-MeV v4.9.0 software (http://mev.tm4.org) to define significantly modulated genes between groups of samples. Differentially Expressed Genes (DEGs) were clustered according to their expression patterns in eight groups as reported in Table 2. For each gene and comparison, the fold change (FC) was calculated as the ratio between the average of the normalized expression data measured in three replicates for test and reference samples. For each gene and comparison, the average of the normalized expression data measured for the sexual samples or the average of the normalized expression data measured for the samples corresponding to the earlier developmental stage (floral stage 11) was set as reference ( Table 2).

Bioinformatics-Array Annotation
Annotation of sequences used for the design of 4x180 arrays (Agilent) is described in Galla et al. (2015). Briefly, to annotate all assembled unigenes, a BLASTX-based approach was used to compare the Hypericum sequences to the nr database downloaded from the NCBI (http://www.ncbi.nlm.nih.gov/). The following approach was adopted to improve our efficacy in annotating 18,386 unigenes that did not reach significance for the BLAST results in the nr database. Unigenes were aligned to the Hypericum genome draft HEID 914233 v1 (Altschmied, personal communication) using BLAT (match ≥ 50 base pairs, fident ≥ 0.97). Additionally, high-quality genomic top hits were selected by preventing the existence of long non-aligned segments on the same side for each unigene and genomic contig. Briefly, a genomic top hits was retained only when the alignment fulfilled the two following criteria (hereafter described for the 5 ′ end of the alignment): (i) 5 ′ unaligned segment of the genomic contig equal or larger than the 5 ′ unaligned segment of the unigene only when the 5 ′ unaligned segment of the unigene was equal or smaller than 10 bp (allowing for an eventually untrimmed adapter sequence at the 5 ′ end of the unigene) (ii) 3 ′ unaligned segment of the genomic contig equal or larger than the 3 ′ unaligned segment of the unigene was allowed only when 3 ′ unaligned segment of the unigene was equal or smaller than 10 bp (allowing for an eventually untrimmed adapter sequence at the 3 ′ end of the unigene) or in cases where the unaligned segment of the genomic contig at the 3 ′ side was smaller than 3 bp. The same criteria were adopted for the 3 ′ ends of each alignment.
Next, high quality top hits were identified, and genomic sequences upstream and/or downstream of the aligned segments of the genomic top hits with min. 100 and max. 1001 bases of genomic sequence were extracted. All extracted sequences were Apospory expressed as a percentage was determined by flow cytometric screening of 48 single seeds. For details on the origin and composition of experimental populations, please refer to Schallau et al. (2010). Accessions marked with • were used for the microarray expression analysis, while those marked with * were used for Real-Time qPCR assays. + The apomictic accessions 222EU was selected for the flower developmental stage 11. ++ The apomictic accessions 3348EU was selected for the flower developmental stage 14. For each comparison the table reports on the main biological process occurring within ovules and the number of identified differentially expressed unigenes (DEGs). Three biological replicates were adopted for both test (t) and reference (r) samples. For each comparison the number of up-and down-regulated unigenes is also indicated.

Expression Correlation Network
An expression correlation network (Pearson correlation coefficient −0.98 <> +0.98) was created from all differentially expressed unigenes (fold change ≥ 2.0) using the plugin ExpressionCorrelation (http://apps.cytoscape.org/apps/ expressioncorrelation) implemented in Cytoscape 3.4.0 (http:// www.cytoscape.org/). Sub-networks were generated from the overall expression correlation network (Pearson correlation coefficient −0.98 <> +0.98; 714 nodes, 8367 edges) by selecting nodes (DEGs) according to their annotations. Next, all adjacent edges and nodes connected to these edges were selected and used to infer the subnetwork. For each sub-network, GO enrichment analyses were performed with the plugin BinGO with default parameters (p < 0.05, Benjamini & Hochberg's FDR correction with significance level 0.05).

Validation of Sequencing Data by Quantitative Real-Time PCRs
Plant materials were selected according to the genetic and cyto-histological bases of apospory recently described for H. perforatum (Schallau et al., 2010;Galla et al., 2011). Samples were collected separately from a minimum of five plant accessions ( Table 1). RNA extractions were performed using a Spectrum TM Plant Total RNA Kit (Sigma-Aldrich). cDNA synthesis was performed using the RevertAid First Strand cDNA Synthesis Kit (Thermo Scientific), following the instructions of the supplier. Primers used in the Real-Time RT-PCR experiments are reported in Table S6. Expression analyses were performed using StepOne thermal cyclers and the 7300 Real-Time PCR System (Applied Biosystems), equipped with 96-and 384-well plate systems, respectively, with SYBR green PCR Master Mix reagent (Applied Biosystems). The amplification efficiency was calculated from raw data using OneStep Analysis software (Life Technologies). Amplification performance expressed as fold change was calculated with the Ct method using HpTIP4 as a housekeeping gene (Pfaffl, 2001). Error bars indicate the standard error observed among the five biological replicates.

Data Availability
Raw sequences files were made available for download from SRA with the following accession numbers: SRR1646951, SRR1646953, SRR1646955, SRR1646956, SRR1647632, SRR1647633, SRR1647673, SRR1647674, SRR1647677, SRR1647678, SRR1647713, and SRR1647714. Sequences of unigenes investigated by Real-Time qPCR along with unigenes aligned to the BAC sequence HM061166.1 were deposited under the Transcriptome Shotgun Assembly project at DDBJ/EMBL/GenBank under the accession GBXG00000000. The version described in this paper is the second version, GBXG02000000.
The expression data discussed in this publication have been deposited in NCBI's Gene Expression Omnibus (Edgar et al., 2002) and are accessible through GEO Series accession number GSE84768 (https://www.ncbi.nlm.nih.gov/geo/query/ acc.cgi?acc=GSE84768).

RESULTS
A transcriptome analysis was performed on pistils collected from sexual and aposporous genotypes ( Table 1) at flower developmental stages spanning female sporogenesis and gametogenesis, i.e., floral stages 11-14 according to Galla et al. (2011). In sexual ovules, these stages span a number of critical developmental steps: the meiotic division of the megaspore mother cell, the degeneration of three reduced megaspores and the selection of a functional megaspore, and three successive mitotic divisions giving rise to a functional eight-nucleate embryo sac with an egg cell apparatus (Galla et al., 2011). Conversely, aposporous ovules are frequently characterized by the failure of the meiotic division followed by developmental program of one or multiple unreduced embryo sacs from aposporous initials differentiating from somatic cells of the nucellus (Galla et al., 2011).

Expression Analysis of Sexual and Apomictic H. perforatum Pistils Reveals Highly Differentiated Transcriptome Signatures at Early Stages of Ovule Development
A principal component analysis (PCA) computed by using DEGs transcriptomic data allowed us to separate the sexual samples from the apomictic samples, with the two first principal components explaining approximately 51.5% of the overall variance ( Figure S1). This analysis revealed that samples at flower developmental stages varying from meiosis (Sex1A, Sex1B, Sex1C) to gametogenesis (Sex2A, Sex2B, Sex2C) were closely related for the sexual genotypes whereas the apomictic genotypes revealed higher level of variation among samples of both early and late developmental stages ( Figure S1). When considering only the expression data related to DEGs annotated as reproductive-related, centroids of the sexual and apomictic samples were plotted and clustered according to a different pattern, with the two first principal components explaining approximately 95% of the overall variance (Figure 1). This analysis revealed very high uniformity of biological replicates for the apomictic samples at both early and late developmental stages, whereas the sexual samples were characterized by higher variability among replicates, especially at the late developmental stages (Figure 1). In particular, apomictic samples collected at the late developmental stages, corresponding to gametogenesis, could be clearly separated from the apomictic samples collected at the early developmental stages, spanning apomeiosis and the differentiation of the aposporous initials. Altogether, the separation among sexual and apomictic samples collected at developmental stages spanning meiosis was very well marked, whereas sexual and apomictic samples collected during gametogenesis were closely related, indicating a lower level of variation in this stage (Figure 1).
Overall, across the two developmental stages, 224 and 973 unigenes were found to be significantly up-and down-regulated with a fold change ≥ 2 in at least one comparison, respectively. For convenience, clusters of unigenes that were up-and downregulated in each comparison were designated C1-C8 ( Table 2). A differential expression analysis was performed to identify genes modulated throughout pistil development in both sexual and apomictic genotypes ( Table 2, clusters C1-C4), while at the same time providing the opportunity to compare sexual and apomictic transcriptomes at comparable pistil developmental stages spanning meiosis, apomeiosis and gametogenesis ( Table 2, clusters C5-C8).
In total, 33 and 137 unigenes were up-regulated between the two selected developmental stages of pistils, in sexual (  genotypes, respectively. Conversely, the numbers of downregulated unigenes in pistils collected at the same developmental stages were 53 and 182, respectively, for sexual (Table 2, cluster C2) and apomictic (Table 2, cluster C4) genotypes.
The comparison between sexual and apomictic pistils spanning early flower developmental stages, which are characterized by the circumvention of the meiotic program and the onset of cellular and biological modifications leading to the differentiation of aposporous initials, underlined as many as 760 unigenes that were differentially expressed ( Table 2, clusters C5 and C6). This comparison provided the highest number of DEGs. We also noticed a marked imbalance in the number of DEGs between sexual and apomictic pistils at early flower developmental stages, with 36 up-and 724 down-regulated, respectively ( Table 2, clusters C5 and C6). Interestingly, the comparison between sexual and apomictic pistils collected at late flower developmental stages (floral stage 14, clusters C7 and C8) was characterized by 32 DEGs, 18 and 14 of which were up-and down-regulated, respectively (Table 2, cluster C8).
Differences in the expression of transcripts encoded by genes located on the BAC clone encompassing the genomic locus associated with apospory (Schallau et al., 2010) were not significant (considering a fold change >2) in the pair-wise comparisons between reproductive strategies and developmental stages examined in this study ( Figure S2). However, the unigene HHKO36P02G7QAY, mapping on the BAC region encoding for the protein HnRNP, was significantly differentially expressed between sexual and apomictic pistils at floral stage 11 (fold change is equal to-1.89).
Venn diagrams were calculated to estimate the overlaps between unigenes differentially expressed in apomictic samples compared to the sexual samples and unigenes which, based on annotations of their homologs, are likely associated with sporogenesis and/or gametogenesis. In doing so, we discovered that 13 unigenes belonging to the clusters C1 and C2 were also included in clusters C3 and C4 (Figure 2A). Likewise, 10 unigenes belonging to clusters C1 or C2 were also included in the clusters C5 and C6. Notably, only the unigene isotig20049, sharing high similarity with the A. thaliana ribosome maturation protein SBDS encoded by the gene AT1G43860, appeared to be upregulated in apomictic samples across the early and late developmental stages (clusters C5/C7).
The ontological annotation of DEGs was performed according to two vocabularies: biological processes and molecular functions. Of the eight datasets, the numbers of differentially expressed genes that were annotated with at least one GO term were 63 (0.73% of clusters C1/C2), 265 (0.83% of clusters C3/C4), 474 (0.62% of clusters C5/C6) and 17 (0.53% of clusters C7/C8). For biological processes, the highest scores were recorded for terms associated with the sensing and response to stresses and stimuli as well as for terms associated with the metabolism of carbohydrates, lipids and secondary compounds ( Table 3). The annotations of anatomical structure morphogenesis (hits: 47, GO:0009653), cellular component organization (hits: 116, GO:0016043), cell cycle (hits: 71, GO:0007049), DNA metabolic process (hits: 80, GO:0006259) and signal transduction (hits: 45, GO:0007165) were also identified to a relatively high extent ( Table 3).
The annotation of DEGs according to the putative molecular functions of corresponding gene products revealed a number of terms assigned with higher frequency (Table 4). In this case, the highest scores were recorded for GO terms associated with nucleotide binding (GO:0000166), DNA binding (GO:0003677), transporter activity (GO:0005215) and protein binding (GO:0005515). Overall, 54 DEGs were annotated as having sequence-specific DNA binding transcription factor activity (GO:0003700). Among DEGs belonging to the clusters C5/C6, which are thus likely reflective of transcriptional changes associated with apomeiosis and cell fate changes needed for apospory initiation, several unigenes were annotated as chromatin binding (hits: 8, GO:0003682), RNA binding (hits: 27, GO:0003723) and nuclease activity (hits: 15, GO:0004518). A GO enrichment analysis allowed the identification of GO terms significantly over-represented in each cluster of DEGs (p <0.05, FDR < 0.05). While clusters 2, 7, and 8 were not enriched for any ontological terms (Table S1), DEGs included in clusters C5 and C6 were enriched for multiple terms, including   (Table S1). For molecular functions, the highest scores were recorded for terms associated with DNA binding (GO:0003677, FDR: 1.31E-10), DNA (cytosine-5-)-methyltransferase activity (GO:0003886, FDR: 3.52E-08) and heterocyclic compound binding (GO:1901363, FDR: 8.46E-07), among others (Table  S1). Notably, despite the abundance of unigenes included in clusters C5 and C6, very few ontological terms appeared to be significantly enriched for annotations associated with the included unigenes. Thus, transcription corepressor activity (GO:0003714, FDR: 2.84E-02) and regulation of jasmonic acid mediated signaling pathway (GO:2000022, FDR: 2. 84E-02) were the only enriched ontological terms in cluster 5, which included 36 unigenes. No ontological annotations related to molecular functions and biological processes were found to be significant in unigenes included in cluster C6. Summarizing, the ontological annotation of DEGs provided clear indications concerning the biological process and molecular functions that are overrepresented in concomitance to the early events of the aposporous developmental pathway. Real-Time RT-qPCR assays were used to corroborate the microarray expression data for a number of key genes, which were selected on the bases of their expression pattern or annotation. Nine genes were selected among those differentially expressed at early developmental stages in apomictic vs. sexual samples ( Figure S3). The piezo-type mechanosensitive ion channel encoded by the unigene isotig06122 exhibited a significant reduction in its expression in apomictic samples at floral stage 11, while no significant difference in its expression was recorded at floral stage 14. The same expression pattern was observed for the unigene isotig06679, orthologous to the Arabidopsis gene IDN2, which was found to be down regulated in apomictic pistils collected at developmental stages encompassing apomeiosis and induction of the aposporous developmental program (i.e., floral stage 11). The expression of unigene isotig06027 (similar to AT5G04390) decreased in apomictic samples in correspondence of both developmental stages. A decreased expression level in aposporous pistils was also observed for the three unigenes: isotig23509 (similar to the C2H2-type zinc finger protein encoded by AT5G04390), isotig16989 (similar to the mithocondrial translocase encoded by AT5G43970) and isotig07554 (similar to the the alpha-expansin 11 encoded by At1g20190) at floral stage 14 ( Figure S3). The two unigenes: isotig12158 and isotig17364, encoding for a MYB7like transcription factor and an ankyrin repeat family protein, showed the same expression pattern in sexual and apomictic samples and appeared to be modulated were to be modulated during pistil development. The unigene isotig17171, encoding for a protein of unknown function, was up-regulated in apomictic pistils in correspondence of both developmental stages.

Phenotypic Expression of Apospory Is Concomitant with the Modulation of Key Genes Involved in the Sexual Reproductive Pathway and the Response to Hormonal Stimuli
We investigated the hypothesis of deregulation of the sexual developmental pathway leading to apomixis (Grimanelli et al., 2001) by identifying DEGs associated with reproduction based on an annotated H. perforatum flower transcriptome (Galla et al., 2015). In doing so, we could annotate 2512 unigenes matching 1042 Arabidopsis gene models as related to plant reproduction, which include a number of DEGs (Table S2).
Six unigenes matching five Arabidopsis gene models were annotated as sporophytic mutants with gametophytic defects (Bencivenga et al., 2012). Despite the fact that gene products are mostly expressed in Arabidopsis sporophytic tissues, these genes are known to lead to gametophytic defects when mutated (Table S3). Additionally, 34 and 25 unigenes were annotated as sporogenesis-and gametogenesis-related, respectively (Table S3). Remarkably, the vast majority of these unigenes was included in clusters C4 and C6 and was down-regulated in aposporous pistils at floral stage 11.
Based on the expression patterns of these genes in sexual and apomictic pistils, 7 and 38 reproductive-related DEGs were included in clusters C1/C2 and C3/C4, respectively ( Figure 2B, Table S3). Furthermore, 43 unigenes that, based on the annotation of their Arabidopsis homologs were annotated with GO terms linked to cell fate commitment, sporogenesis, gametogenesis or embryogenesis, were differentially expressed (p < 0.01, FC ≥ 2.0) in the C5-C8 comparison (Table S3). Of these 38 were down-regulated in pistils of aposporous genotypes (C6/C8), whereas five were significantly up-regulated at both floral developmental stages 11 and 14 (C5/C7).
Next, in the attempt to identify additional levels of regulation acting during sporogenesis and gametogenesis, the H. perforatum flower transcriptome was annotated with terms for hormonal regulation, whereby annotations were provided by TAIR. Within the flower transcriptome, 2436 unigenes were annotated as responsive to hormonal stimuli (Table S2).
Further categorization of DEGs based on their commitment to hormonal biosynthesis, metabolism and signaling pathways is reported in Table 5. Among these genes, a restricted set of DEGs were involved in cytokinin biosynthesis, such as the tRNA dimethylallyltransferase 2 HPIPT2, and auxin For each unigene the table reports on its involvement in hormonal biosynthesis and/or perception. The Arabidopsis gene sharing the highest similarity, the gene name, the protein name, annotations related to reproduction and differential expression (gene expression cluster and fold change) are also indicated. homeostasis, such as HPPIN8, HPIAR3 and HPILL6. Similarly, we identified a number of genes that, based on their annotations, could be recognized as regulators of hormonal stimuli (Table 5). Remarkably, the two cytokinin receptorlike proteins HPCRE1 and HPACR4 were down-regulated in aposporous pistils at floral stage 11 (cluster C6, Table 5).
The overlaps existing between DEGs differentially expressed in aposporous vs. apomictic samples (clusters C5-C8), those annotated as reproduction-related (namely, "reproductiverelated") and unigenes thought to be responsive to hormonal stimuli (namely, "response to hormonal stimuli") are displayed in Figure 2C. While the majority of DEGs were not annotated according to the abovementioned categories, 7 DEGs matching 6 A. thaliana gene models were included in both datasets (Table 5, Figure 2C). These unigenes were annotated as sporophytic mutants with gametophytic defects (HPACR4) and regulators of megasporogenesis (HPTRN2) as well as regulators of megagametogenesis progression and gametophyte development (HPBRI1, HPJAZ1). The Hypericum homologs of HPFER and HPARF2, the former of which is involved in pollen tube perception and the latter of which is involved in gametophyte and seed development in A. thaliana, were also among these genes ( Table 5).
Real-Time qPCR assays showed that the unigene HPACR4 (F43FOUR02DAHH0), encoding for a protein involved in epidermal cell differentiation, exhibited a significant reduction in its expression in apomictic samples in correspondence of both floral stages ( Figure S4). The unigene isotig28540, orthologous to the Arabidopsis gene PIN8, exhibited an opposite expression pattern as it was significantly upregulated in apomictic pistils in correspondence of both floral stages ( Figure S4). Conversely, no significant differences were recorded by Real-Time RT-qPCR assays for HPIPT2 (isotig31470), HPCRE1 (F6Z56EY03HBNM2) and HPAMP1 (HHKO36P01AWLS2) (Table S7, Figure S4). Summarising, our expression analysis provides line of evidences suggesting that early events of the aposporous developmental pathway are concomitant with the modulation of multiple genes involved in the sexual reproductive pathway and hormonal homeostasis.

The Transcriptome of Sexual and Aposporous Pistils Is Characterized by Distinct Expression Correlation Networks
An expression correlation network (Pearson correlation coefficient −0.98 <> +0.98) was created starting from all differentially expressed unigenes (fold change ≥ 2.0). Subnetworks were generated from the original network by selecting nodes (i.e., unigenes) modulated in single comparisons, whereby highly correlated unigenes were connected to each other (Table S4).
An expression correlation network (Pearson correlation coefficient −0.98 <> +0.98) was created starting from unigenes (nodes) included in clusters C1, C2, C5, and C6, thus possibly reflecting a subset of expression changes concomitant to the failure of the meiotic program and the induction of aposporous apomixis. This approach was selected to scale down the complexity of data, thus allowing closer investigations of selected nodes. The resulting expression correlation network contains two sub-networks for a total of 26 nodes (unigenes) and 29 edges (Figure 3, Table 6). All edges included in the network represent positive correlations between unigenes that were found down-regulated in aposporous pistils. For these genes, fold changes varied between a minimum value of 2.01 for the unigene of unknown function HHKO36P01BR1EE and a maximum fold change value equal to 13.68 for the mitochondrial ribosome-associated GTPase 1 encoded by the unigene F6Z56EY02DHSUQ. In addition to being differentially expressed during early stages of pistil development between sexual and apomictic genotypes, five unigenes included in the network were significantly down-regulated during sexual development (cluster C2) ( Table 6).
Based on their annotations, unigenes included in the network could be further classified into the following ontological categories: (i) transcriptional regulators, including SPL (SQUAMOSA promoter-binding protein-like transcription factor), BSD domain (BTF2-like transcription factors), the nuclear transcription factor Y subunit B-4-like and a chromatin remodeling complex subunit (CHC1); (ii) translation regulators, including the PAB2 Poly(A)-binding protein and a mitochondrial ribosome-associated GTPase 1; (iii) cell cycle regulators, such as SKP-like 2 and mitotic-like cyclin 3B; (iv) proteins involved in signal perception and transduction (the Serine/threonineprotein kinase (BRI1) was included in this group ( Table 6); v) unigenes responsive to internal or external stresses and stimuli, such as glutathione S-transferase TAU 19 (GSTU19), cytochrome P450-like gene CYP71B and ATP synthase D chain (ATPQ); and (vi) genes encoding for proteins involved in cellular modification processes, such as the pectinesterase PME1 and glucan synthase-like 4 (GSL04). Based on annotations available at UniProt, this latter gene could also be included in the group of genes that are known to be involved in plant reproduction together with BRI1 (Table 6).
Remarkably, the nodes sharing the highest numbers of interactions within the sub-network were isotig18376 (HPBAT1) and HHKO36P02G29PD, this latter encoding for a putative retroelement pol polyprotein. We noticed that transcription factors included in this picture were localized in marginal positions. The only exception to this was represented by the SQUAMOSA promoter-binding protein-like transcription factor 6 (HPSPL6), which is connected to both central nodes (HPBAT1 and HHKO36P02G29PD). Interestingly, physical interaction between the two proteins ASK2 and CYC3B has been reported in A. thaliana (http://bar.utoronto.ca/interactions).
An enrichment analysis of GO terms associated with the nodes included in the network did not reveal over-or underrepresented ontological terms (p < 0.05, Benjamini & Hochberg's FDR correction with significance level: 0.05).
Sixteen genes were selected among those included in the expression correlation network (displayed in Figure 3 and  all nodes (DEGs) included in clusters (C1 OR C2) AND (C5 OR C6) were selected. Next, all adjacent edges and nodes connected to these edges were selected and used to infer the sub-network. For each unigene the table reports on the putative Arabidopsis orthologous, the gene name and the description of the gene. The relevant biological processes attributed to each unigene are also indicated. Table 6) and their expression tested by Real-Time RT-qPCR assays. From this group of DEGs, 12 unigenes displayed a similar expression pattern as they were significantly downregulated in aposporous pistils in correspondence of both developmental stages (Figure 4). In particular, this expression pattern was observed for the unigene HHKO36P01A0D2H, which shares some similarity with a MORC family CWtype zinc finger protein 4, the unigene isotig25064, which is predicted to encode a protein that belongs to the chromodomain remodeling complex (HPCHC1), and three unigenes encoding for transcription factors belonging to the families NF-YB4, BTF2-like, and HIT-type Zinc finger family protein ( Table 6, Figure 4). The same expression pattern was observed for the two unigenes: HPPTR2 and HHK036P02G29PD (Figure 4), the latter of which encodes for putative retroelement pol protein ( Table 6).
In contrast, the unigene HPBRI1 was down-regulated in aposporous pistils only in the terminal developmental stage ( Figure S5). Finally, Real-Time RT-qPCR assays could not confirm the differential expression detected with the microarray analysis for the unigenes: HPALADIN (isotig13015), HPBAT1 (isotig18376), and HPCYC3B (HHKO36P02GL54X) ( Figure S5 and Table S7).

Expression Analysis of Sexual and Apomictic H. perforatum Genotypes Reveals Highly Differentiated Transcriptome Signatures at Early Stages of Pistil Development
From the reproductive point of view, H. perforatum is an aposporous apomictic species whereby the alternative differentiation of a somatic cell gives rise to a functional embryo sac containing an unreduced egg cell. In principle, aposporous initial cells in apomictic plants are somatic cells belonging to the nucellus, which change their fate through the ability to mitotically divide and develop functional embryo sacs by mimicking sexual gametogenesis (Barcaccia et al., 2007;Galla et al., 2011). While the unreduced egg cell has the ability to undergo parthenogenesis, fertilization of the binucleated central cell of an aposporous embryo sac by a pollen nucleus is normally required for endosperm formation (Matzk et al., 2001). Expression analyses were performed on Hypericum pistils collected from genotypes exhibiting opposing reproductive forms at flower developmental stages spanning megasporogenesis and megagametogenesis. These developmental stages were selected to investigate expression changes associated with the transition from sporogenesis to gametogenesis in sexual pistils, concomitant to apomeiosis and development of aposporous initials in aposporous pistils.
Overall, considering the expression of all DEGs ( Figure  S1), apomictic samples were characterized by high variability in terms of the expression levels of multiple genes, leading to high differentiation of samples at both stages. Conversely, a principal component analysis revealed higher uniformity in sexual samples, whereby biological replicates clustered more closely together. Noticeably, restricting the analysis to unigenes which, based on their annotations, are potentially related to reproduction (Figure 1) allowed for a closer clustering of apomictic samples in correspondence of both developmental stages. High variability in overall gene expression observed for the aposporous pistils collected at stages corresponding to apomeiosis and development of the aposporous initials ( Figure  S1) might be reflective of the asynchronous development of these processes. The comparison between apomictic and sexual pistils at stage 14 underlined the relatively low number of genes that were significantly differentially expressed. Conversely, highly differentiated transcriptome signatures characterized sexual and apomictic pistils collected at flower stage 11. This flower developmental stage in apomictic genotypes encompasses the failure of the meiotic program and changes in cell fate required for the induction of aposporous embryo sac development. We are aware that the complexity of the analyzed tissue is such that multiple genes could be modulated in response to processes not directly connected to ovule and gamete development. Nevertheless, the switch from sexuality to aposporous apomixes, that is phenotypically characterized by dramatic changes in ovule and gamete development, was reflected in the modulation of the expression of 760 genes between sexual and aposporous samples. Remarkably, approximately 95% of genes modulated in this developmental stage were down-regulated in aposporous samples. This pervasive decrease in the expression of many genes could be reflective of a change in the rate or timing of development of some cell lines in the ovule relative to others (Galla et al., 2011). Temporal shifts in the expression of genes involved in the sexual developmental program have been already proposed in both diplosporous (Grimanelli et al., 2003) and aposporous (Galla et al., 2015) plant species. Interestingly, a global decrease in gene expression within the apomict compared to the sexual ovules during the early stages of development was also reported for the diplosporious Boechera divaricarpa (Sharbel et al., 2010). The extent to which de-regulation of genes required for sexual gametes formation is shared among aposporous and diplosporous species is currently not known. Nevertheless, this finding pinpoints an intriguing similarity in overall transcriptome changes occurring in concomitance with embryo sac development from somatic cells and unreduced functional megaspores, as observed in both aposporous and diplosporous species.
Ontological annotation of DEGs identified in each pairwise comparison led to the identification of several classes of annotation associated with flower development, the regulation of cell cycle progress and DNA metabolic processes. Along with these annotations, multiple terms were related to responses to biotic and abiotic stimuli and carbohydrate metabolic processes. Similarly, Schmidt et al. (2014) reported that microdissected aposporous initial cells of Boechera gunnisoniana were enriched for terms associated with carbohydrate transport, membrane lipid biosynthetic process and terms associated with the biosynthesis and transport of amino acids. Modulation of the genes involved in responses to biotic and abiotic stimuli as well as carbohydrate and lipid metabolisms in sexual vs. apomictic accessions was previously documented by Galla et al. (2015). Interestingly, although phenotypic changes associated with apomeiosis of the apospory-type are apparently restricted to a limited number of cells within the ovule, these highly localized developmental alterations are accompanied by the modulation of large transcriptome portions ( Table 2), involving the sensing (i.e., response) of external (exogenous biotic or abiotic) and internal (i.e., endogenous) stimulus (Table 3) and the generation or catabolism of precursor metabolites and energy ( Table 3). More recently, Klatt et al. (2016) observed that photoperiod extension, which could be regarded as a moderate abiotic stress, has a significant effect on the rate of megasporogenesis in the facultative aposporous apomictic Ranunculus. Future researches will be aimed at understanding if the failure of the meiotic program, the lack of a functional megaspore and/or the induction of apospory are sensed by the sporophytic cells embracing the gametophyte as internal stimuli, which in turn could directly or indirectly elicits wide transcriptional changes in the ovule.

Induction of Apospory Is Concomitant with the Differential Expression of Key Genes Involved In the Sexual Reproductive Pathway and RNA-Directed DNA Methylation
Multiple genes involved in female sporogenesis were significantly modulated throughout pistil development. Strikingly, the vast majority of sporogenesis-related DEGs were modulated in aposporous pistils during their development (cluster C4) and compared to their sexual counterparts (cluster C6) while only HPANT, the putative bipolar kinesin KRP-130 and HPRAD54 were down-regulated in sexual pistils (cluster C2).
To explain this finding, several scenarios can be taken into account. While it is possible that some genes not included in the array are involved in sporogenesis and gametogenesis, the present design of the array includes over 2500 unigenes that, based on their annotation, were predicted to be involved in ovule development and gamete formation (Galla et al., 2015). Otherwise, we could hypothesize a relatively low level of synchronism in sporogenesis progression alongside the proximal-distal axis of the pistil (Galla et al., 2011), leading to high variability within and among biological replicates and a concomitant decrease in our statistical power to detect significantly modulated genes. Interestingly, this hypothesis would imply that once established, transcriptional changes linked to the early events of aposporous apomixis (i.e., apomeiosis and differentiation of the aposporous initial) might be within view for a relatively long time.
Interestingly, the analysis of ontological terms significantly overrepresented in the cluster of genes down regulated in aposporous pistils at stage 14 (clusters C3 and C4) underlined a significant enrichment of terms associated with DNA binding and DNA-methyltransferase activity, which was not observed in sexual pistils collected at comparable developmental stages. Recently, Garcia-Aguilar et al. (2010) reported that the DNA methylation pathway active during reproduction is essential for gametophyte development in maize and likely plays a critical role in the differentiation between apomictic and sexual reproduction. In the same study, the authors observed that lossof-function mutants for two genes involved in the maintenance of DNA-methylation status result in phenotypes that are strikingly reminiscent of apomictic development, suggesting that, in addition to a crucial role in gametophyte development, DNA methylation in the maize ovule might regulate the transcriptional expression of genes involved in the differentiation between apomixis and sexual reproduction (Garcia-Aguilar et al., 2010). Similarly, a study examining the distribution of cytosine methylation at the apomixis locus in two aposporous Paspalum species revealed that artificial demethylation induced a significant depression of parthenogenesis, suggesting that factors controlling the repression of parthenogenesis might be inactivated in apomictic Paspalum by DNA methylation (Podio et al., 2014). In this frame, it is worth noting that the two genes involved in the maintenance of DNA methylation at CpG and non-CpG sites: MEE57 (Borges and Martienssen, 2013) and CMT3 (Henikoff and Comai, 1998), decreased their expression throughout pistil development in aposporous but not in sexual samples. It is worth noting that CMT3 is a key determinant for CpNpG methylation in A. thaliana (Lindroth et al., 2001) and its expression is critical for the maintenance of methylation at transposon-related sequences (Kato et al., 2003). Furthermore, its expression is also required for H3K9 dimethylation in the egg cell and for normal embryogenesis during the first few divisions of the zygote (Pillot et al., 2010). Intriguingly, CMT3 is also involved in the RNA directed DNA methylation pathway in Arabidopsis (Kawashima and Berger, 2014), suggesting that differentiation of the aposporous initial might be somehow related to miss regulation of the RNA-directed DNA methylation (RdDM) pathway or improper spatial/temporal expression of its components. Further experimental evidence that RNA directed DNA-methylation changes might be critical for the switch from sexual to aposporous developmental pathways derives from the finding that an apospory-like phenotype is associated with argonaute 9 (ago9) mutants in Arabidopsis thaliana (Olmedo-Monfil et al., 2010). AGO9 preferentially interacts with 24-nucleotide small interfering RNAs (siRNA) derived from transposable elements to direct homolog-based RNA-directed DNA methylation (RdDM). Moreover, mutations in SUPPRESSOR OF GENE SILENCING3 (SGS3) and RNA-DEPENDENT RNA POLYMERASE6 (RDR6), two genes required for siRNA biogenesis, also lead to a defect identical to that in ago9 mutants (Olmedo-Monfil et al., 2010). In a previous study, we demonstrated that the expression of HPAGO9 is significantly reduced in aposporous pistils at floral stage 11, while the expression of HPSGS3 and HPRDR6 has no or little variation between sexual and aposporous pistils (Galla et al., 2015). These findings were confirmed by the present expression analysis. Additionally, our microarray data and Real-Time qPCR assays indicated that the unigene isotig06679, which encodes for an orthologue of Arabidopsis INVOLVED IN DE NOVO 2 (IDN2), a doublestranded RNA-binding protein involved in de novo methylation and small interfering RNA (siRNA)-mediated maintenance of methylation (Finke et al., 2012), was significantly downregulated in aposporous pistils at stage 11 compared to sexual pistils (FC: −10.42). Previously, Zhu et al. (2013) demonstrated that IDN2 binds to Pol V-produced noncoding RNAs and physically interacts with MORC family proteins and SWI3B, a subunit of the SWI/SNF chromatin-remodeling complex, thus contributing to noncoding RNA-mediated transcriptional silencing. Noteworthy, our microarray and Real-Time qPCR assays showed that the unigene HHKO36P01A0D2H, sharing high similarity with a MORC-like CW-type Zinc finger protein, is significantly down regulated in aposporous pistils at both floral developmental stages. Additionally, the expression of the unigene isotig25064, which is predicted to encode a SWI/SNF ASSOCIATED PROTEINS 73 (HPCHC1), was also found to be reduced in aposporous pistils (Table 6, Figure 4). The two unigenes: isotig25064 and HHKO36P01A0D2H were both included in the expression correlation sub-network (Table 6, Figures 3, 5) whose modulation in aposporous pistils is concomitant with apomeiosis and aposporous initial induction and development. Interestingly, Jégu et al. (2015) demonstrated that CHC1 controls cytokinin biosynthesis and cell cycle progression in the Arabidopsis root, by directly modulating the deposition of H3K27 and H3K4me3 histone marks as well as RNA Pol II occupancy at the IPT3 and IPT7 loci. Whether the same actors regulate the expression of HPIDN2 (unigene: isotig06679), the MORC-like CW-type Zinc finger protein (unigene: HHKO36P01A0D2H) and HPCHC1 (isotig25064) and what is nature of these regulators is currently elusive. Nevertheless, the finding that expression of HPCHC1 and the MORC-like CW-type Zinc finger protein is highly correlated with that of several transcription and cell cycle regulators (Table 6, Figure 4) might lead us to speculate that possible regulators or targets of the RdDM machinery in aposporous pistils could be included among these co-expressed genes.

Unigenes Involved in Homeostasis, Signalling and Response to Hormones Are Differentially Expressed in Aposporous Pistils
Cytokinins and auxins are classes of plant hormones with a broad range of functions in regulating plant growth and development as well as physiological processes (Werner and Schmülling, 2009;Kieber and Schaller, 2014). It is known that activity of the two hormones in promoting stem cell maintenance or cell differentiation is defined based on the specific environmental context (Bencivenga et al., 2012). Key components of cytokinin biosynthesis and response modulators are expressed during ovule development (Bencivenga et al., 2012). Overall, the ontological term "cytokinin catabolic process" was found to be significantly enriched in B. gunnisoniana aposporous initial cells compared to the cell types comprising the mature gametophyte (Schmidt et al., 2014). The expression data collected in sexual and aposporous pistils indicates that several genes involved in the biosynthesis and perception of cytokinin stimuli are downregulated in aposporous pistils. Among these genes, IPT2, a tRNA isopentenyl transferase whose mis-expression has been associated with meiotic defects in Arabidopsis, was down regulated in aposporous pistils. This finding is coherent with the observed reduction in the expression of HPCHC1 (FC: −2.27), which is known to control cytokinin biosynthesis by regulating the expression of several IPT genes.
Interestingly, among genes responsive to cytokine stimuli, the receptors and/or signal transducer ACR4 (Hsieh and Goodman, 2002) was down regulated (FC: −2.27) in aposporous pistils at stage 11. It has long been known that an abundance of the ACR4 transcript is increased in response to cytokinin (Hsieh and Goodman, 2002). It is worth noting that the expression of the homeodomain transcription factor BEL1, one of the major players in chalaza development, which is modulated in response to cytokinin levels in Arabidopsis (Bencivenga et al., 2012), was increased in Hypericum aposporous pistils at developmental stages spanning late gametogenesis by a fold change of 1.88.
As cytokinin was proposed to be involved in ovule development by modulating auxin fluxes through the control of PIN1 expression, we were also interested in the expression patterns of genes involved in auxin biosynthesis, homeostasis and responses. Among the auxin efflux carriers involved in ovule intercellular auxin gradients and maxima, such as PIN1, no significant changes in expression were observed in sexual and/or aposporous pistils. The same was true for the unigenes encoding HPPIN3 and HPPIN5. Instead, Real-Time qPCR assays revealed that the endoplasmic reticulum-localized HPPIN8 is up regulated in aposporous pistils at flower stage 11 and 14. Following the publications of Ding et al. (2012) and Dal Bosco et al. (2012), we have indications that PIN8 is involved in the regulation of intracellular auxin homeostasis and affects auxin-regulated gene transcription by either preventing access of auxin to the nucleus or altering the biologically active auxin pool within the nucleus.
Taken together, the present data indicate that modulation in the expression of genes involved in cellular auxin homeostasis, either by intracellular transport or hydrolysis of IAA-conjugates, is concomitant with aposporous gametophyte development. The finding that sexual gametophyte development is not associated with modulation of the same genes might lead to speculate a possible role for these proteins in modulating IAA cellular levels in ovules undergoing the aposporous developmental program. Accordingly, several genes regulated by auxin, some of which are known to be involved in meiosis and/or gametogenesis, were differentially expressed in our samples. Among these genes, the expression of ARF2, ARF7, a SAUR-like auxin-responsive family protein (homologous to AT5G20820), TRN2, the RNA polymerase II C-terminal domain phosphatase-like 2 (CPL2) and the SKP1-like protein 1B (ASK2) was significantly lower in the aposporous than the sexual pistils. Conversely, the auxinresponsive protein IAA4 and the SAUR-like auxin-responsive protein SAUR20 were up regulated in aposporous pistils (stage 14 vs. stage 11).

CONCLUDING REMARKS
Expression analyses were aimed at investigating expression changes associated with the transition from sporogenesis to gametogenesis in sexual pistils, concomitant with apomeiosis and development of aposporous initials in aposporous pistils. Data collected are coherent with the asynchronous or heterochronic expression of multiple genes in aposporous pistils, a phenomenon which, for the most part, could be restricted to the stage corresponding to apomeiosis and development of the aposporous initials. Ontological annotation of DEGs led to the identification of several classes of annotation associated with the regulation of cell cycle and DNA metabolic processes along with multiple terms related to responses to biotic and abiotic stimuli (Table 3, Figure 5). The three genes involved in the RdDM pathway: HPCHC1, IND2, and CMT3 along with a MORC-like protein likely involved in the same biological process, were significantly down regulated in aposporous pistils. This finding provides further evidences supporting the already proposed idea that phenotypic expression of apomixis is dependent upon the chromatin state in sporophytic cells embracing the gametophyte (Garcia-Aguilar et al., 2010). Furthermore, the annotation of unigenes included in the expression correlation network reported in Figure 5 provides evidences for a possible link between genes involved in the RdDM pathway (i.e., HPCHC1), stress responses and flower development. Remarkably, a number of genes involved in hormones perception and homeostasis, including cytokinin, auxin, and brassinosterois, were differentially expressed in aposporous pistils. The connection between hormones and the switch from the sexual to the aposporous developmental programs is further strengthened by the finding that cytokinin biosynthesis in Arabidopsis is regulated by the activity of CHC1 (Jégu et al., 2015). These findings provide new evidences in support of the idea that chromatin state in sporophytic cells embracing the gametophyte might influence the cellular levels of plant hormones and ultimately, based on the specific environmental context, modulate the expression of a wide range of genes, including reproductive and cell fate related genes, in response to internal or external stimuli.

AUTHOR CONTRIBUTIONS
GB and GG conceived the study. SZ and LAV performed the Array design, cDNA synthesis and hybridization. GG carried out the computational analysis, validated the expression data and drafted the manuscript. LAL helped to annotate the Hypericum unigenes. All authors helped to draft the manuscript. All authors read and approved the final manuscript.