A study of the heterochronic sense/antisense RNA representation in florets of sexual and apomictic Paspalum notatum

Apomixis, an asexual mode of plant reproduction, is a genetically heritable trait evolutionarily related to sexuality, which enables the fixation of heterozygous genetic combinations through the development of maternal seeds. Recently, reference floral transcriptomes were generated from sexual and apomictic biotypes of Paspalum notatum, one of the most well-known plant models for the study of apomixis. However, the transcriptome dynamics, the occurrence of apomixis vs. sexual expression heterochronicity across consecutive developmental steps and the orientation of transcription (sense/antisense) remain unexplored. We produced 24 Illumina TruSeq®/ Hiseq 1500 sense/antisense floral transcriptome libraries covering four developmental stages (premeiosis, meiosis, postmeiosis, and anthesis) in biological triplicates, from an obligate apomictic and a full sexual genotype. De novo assemblies with Trinity yielded 103,699 and 100,114 transcripts for the apomictic and sexual samples respectively. A global comparative analysis involving reads from all developmental stages revealed 19,352 differentially expressed sense transcripts, of which 13,205 (68%) and 6147 (32%) were up- and down-regulated in apomictic samples with respect to the sexual ones. Interestingly, 100 differentially expressed antisense transcripts were detected, 55 (55%) of them up- and 45 (45%) down-regulated in apomictic libraries. A stage-by-stage comparative analysis showed a higher number of differentially expressed candidates due to heterochronicity discrimination: the highest number of differential sense transcripts was detected at premeiosis (23,651), followed by meiosis (22,830), postmeiosis (19,100), and anthesis (17,962), while the highest number of differential antisense transcripts were detected at anthesis (495), followed by postmeiosis (164), meiosis (120) and premeiosis (115). Members of the AP2, ARF, MYB and WRKY transcription factor families, as well as the auxin, jasmonate and cytokinin plant hormone families appeared broadly deregulated. Moreover, the chronological expression profile of several well-characterized apomixis controllers was examined in detail. This work provides a quantitative sense/antisense gene expression catalogue covering several subsequent reproductive developmental stages from premeiosis to anthesis for apomictic and sexual P. notatum, with potential to reveal heterochronic expression between reproductive types and discover sense/antisense mediated regulation. We detected a contrasting transcriptional and hormonal control in apomixis and sexuality as well as specific sense/antisense modulation occurring at the onset of parthenogenesis.

Background Apomixis (i. e., agamospermy) is an asexual mode of plant reproduction via seeds, which generates progenies consisting of exact genetic replicas of the mother plant [1]. This puzzling trait occurs in at least 326 angiosperm genera, with no clear tendency to any specific group [2]. The apomictic and sexual developmental pathways are strongly related, since both take place within the ovule and involve common developmental features [3]. Traditionally, the asexual route was considered to be a deviation from the sexual one, repeatedly emerging during evolution from genetic or epigenetic mutations derived from polyploidization and/or hybridization events [4][5][6]. More recently, apomixis and sexuality were hypothetically classified as ancient polyphenic traits [7]. According to the latter view, the switching from one phenism to the other would be environmentally triggered by epigenetic mechanisms, with full sexual genera/species having lost the capacity to carry on this transition due to deleterious (epi)mutations affecting the molecular switch that connects both pathways [7].
The potential of apomixis for fixing heterosis has long been recognized [8,9]. Therefore, the use of this trait in combination with sexuality can dramatically accelerate the development of new hybrid varieties and reduce the costs associated with seed-production [10]. Among other advantages, the transference of apomixis to species of agricultural interest will allow the perpetuation of genetic resources including wide-cross hybrids, the rapid generation of adapted crops, the avoidance of monocultures and the development of maternal seeds from vegetatively propagated cultivars, like potatoes or strawberries [10][11][12]. The apomixis breeding technology impact in global agriculture could be comparable to that produced by the green revolution, initiated in USA, Mexico, India and further spread to other countries in the middle '70s [13].
While plant sexuality takes off with a specialized cell division process (meiosis) preceding the formation of haploid megaspores, apomictic mechanisms share the common characteristic of lacking any reductive division before female sporogenesis (apomeiosis) [1]. Moreover, while sexuality start the sporophytic life cycle by restoring the species-specific ploidy level through fertilization of the genetically reduced egg cell with an equally reduced male gamete, apomixis does it by inducing parthenogenesis, i. e., the spontaneous formation of an embryo from a reproductively-committed cell [1]. Finally, for the formation of the seed endosperm, the sexual route requires the fertilization of two reduced female polar nuclei with a reduced male gamete under a strict 2:1 maternal:paternal genomic contribution. Instead, apomixis may alternatively proceed with the spontaneous proliferation of maternal polar nuclei (autonomy) or the fertilization of one/two female unreduced polar nuclei with a male reduced gamete (pseudogamy), a path that often deviates from the expected genome contribution ratio [1]. While in some species the three components of apomixis (apomeiosis, parthenogenesis and endosperm development) seem to be controlled by a non-recombinant superlocus, in others these factors can be readily uncoupled [2,3,[14][15][16]. However, in nearly all cases there seem to be a consistent association between the expressivity of trait and the ploidy level increments [2,3,[14][15][16].
Paspalum notatum Flüggé is a warm-season perennial grass widely distributed in the Western Hemisphere [46], where it occurs as a primary constituent of natural grasslands, particularly in southern Brazil, Paraguay, Uruguay, and north-east Argentina [47]. The species form a multiploid complex in which the diploid cytotype (2n = 2x = 20) is self-sterile and sexual, while the tetraploid one (the common race) (2n = 4x = 40) is pseudogamous, aposporous apomictic and self-fertile [48]. Other infrequent polyploid cytotypes (3x and 5x) are also apomicts [49]. Moreover, numerous sexual tetraploid individuals were artificially synthetized from diploids by colchicine treatment, or obtained from experimental crosses involving facultative apomicts [50]. Interestingly, although sexual seeds form the endosperm under a strict 2 maternal: 1 paternal genomic contribution, apomictic ones are more permissive and can develop under different maternal:paternal genomic ratios (e.g., 4:1; 8:1) [49]. P. notatum has become a model system for apomixis research and breeding, mainly due to the existence of freely-crossable races of the same ploidy and different reproductive mode, a thoroughly-characterized living germplasm collection and advanced breeding programs exploiting apomixis for cultivar generation [48,[50][51][52][53].
Several transcriptome surveys were conducted to identify P. notatum apomixis-associated genes. Firstly, transcripts expressed in sexual and apomictic florets were compared by using differential display, allowing the identification of a pioneer apomixis candidate gene homologous to the maize kinesin-like motor protein KIN-14P [28]. Then, Laspina et al. [29] identified 65 transcripts differentially expressed in spikelets of apomictic and sexual genotypes at premeiosis/meiosis, several of which mapped in silico to a rice chromosome 2 region that had previously been associated with apospory by comparative mapping [54,55]. Moreover, endosperm RNA representation comparisons 3-24 h after pollination revealed more than 100 differentially expressed transcripts (DETs) in seeds that differed from the expected 2 m:1p genome contribution ratio, formed when apomictic plants were used as female parents in crosses (the endosperm involved 4 m:1p, 8 m:1p or 8 m: 3p contribution ratios, depending on the cross) [56]. Besides, transcripts related to endosperm development were identified in apomictic and sexual ovaries of Paspalum notatum 48 h after pollination, a stage prior to post-zygotic collapse [57]. These DETs were mainly associated with genes related to transcription, signal transduction, growth/division, protein destination and storage, as well as regulation of gene expression. Interestingly, several differential sequences identified at the onset of endosperm development showed high similarity with proteins expressed in response to changes in the levels of extracellular ATP [56,57]. Besides, the Roche 454/FLX + long-read technology was used to produce apomictic and sexual reference floral transcriptomes on an equitable mix of RNA extracted from spikelets at different developmental states from premeiosis to anthesis [31]. Recently, De Oliveira et al. (2020) [58] reported a global gene expression analysis using Illumina Hi-Seq technology on RNA isolated from leaf and floral tissues of 2x sexual, 4x sexual and 4x apomictic genotypes. Interestingly, 89 DETs expressed in apomictic or sexual plants mapped at the chromosome regions of rice syntenic to the Paspalum apomixis controlling locus (ACL) [58].
All the above-cited contributions have partially disclosed the nature of numerous apomixis candidates and evidenced the effects of polyploidy on gene expression. However, our knowledge on the chronological modulation of transcript levels along the sexual and the apomictic reproductive routes remained limited, since previous studies were conducted on samples collected at a particular timeframe or on mixed pools representing a group of developmental stages. Under these experimental limitations, either part of the DETs go undetected or those experimenting expression increments at different developmental stages for contrasting reproductive modes are mistakenly classified as non-differential. Moreover, the plus/minus orientation of the expressed transcripts remains globally unexplored, even when some apomixis candidates (ORC3, PsACR/H5, PsACR/H.13) display antisense differential expression [39,59]. The lack of comprehensive data prevents researchers from grasping the true dimensions of heterochronic expression and antisense regulation affecting apomixis development as well as from inspecting their biological consequences. The objective of this work was to gain a quantitative, statistically significant, massive characterization of sense/ antisense transcripts expressed across four crucial reproductive steps in sexual and apomictic P. notatum and, after reciprocal comparisons, produce a detailed picture of the main molecular pathways operative during apomixis.

Sequencing and de novo assembly
To get a compressive characterization of sense/antisense transcripts expressed during the P. notatum reproductive development, Illumina TruSeq floral cDNA libraries representing two reproductive modes (apomixis and sexuality) and four developmental stages (premeiosis, meiosis, postmeiosis, anthesis), each one including three biological replicas, were sequenced with Illumina HiSeq 1500 technology. The procedure involved 24 libraries and generated a total of 60.94 Gb, of which 97% had a Phred value > 30. After demultiplexing, 292,647,558 pass filter (PF) reads were selected, which, after cleaning and trimming, yielded 234,957,559 high-quality reads (Q score > 30) (available under the NCBI SRA accession PRJNA511813). De novo assemblies were carried out with the Trinity software [60,61], considering two groups of samples (apomictic and sexual), each of them containing 12 libraries. In a first trial, the available Roche 454/FLX + P. notatum reference floral transcriptome [31] was used as a guide reference. Then, a second de novo assembly, without any reference, was constructed to detect novel transcripts. Table 1 shows the output basic statistics derived from both procedures. As expected, the de novo assembly without a guide reference generated a higher number of contigs for both samples, probably reflecting the inclusion of low expressed transcripts and/or allelic variants absent in the Roche 454/FLX + transcriptome. Subsequently, the four assemblies were concatenated in one file and filtered to obtain 199,074 non-redundant transcripts corresponding to a global transcriptome assembly (GTA), representing both reproductive modes (apomixis and sexuality), and the four developmental steps (premeiosis, meiosis, postmeiosis, anthesis) ( Table 1). This transcriptome shotgun assembly was deposited at DDBJ/EMBL/Gen-Bank under the accession GIUR00000000. With an average length of 1182.31 bp and an N50 of 1508 bp, the assembled sequences resulted in good quality for annotation ( Table 1). Mapping of raw reads onto the GTA with Bowtie2 and TopHat showed 99.09 and 98.3% of match with the reference, respectively. These results indicated the GTA covered almost the complete set of sequence reads. To validate the assembly identity, a BLASTN assay was performed onto the Roche 454/FLX + reference transcriptome [31]. Ninety-four % of the transcripts showed homology with the reference and 92% of them exhibited more than 95% identity (Fig. 1). A survey of the assembled contigs with TransDecoder identified 108,011 (54.24%) entities with protein-coding capacity, while BUSCO estimated a gene coverage of 93.1% (35.6% single copy and 57.5% duplicated genes). This value was higher than the coverage achieved with the Roche 454/FLX + reference transcriptome [31], which rounds about 82.2% (Additional file 1).

Global transcriptome assembly (GTA) annotation
A total of 101,079 transcripts produced robust top BLASTN hits against the NCBI NT database (E-value: 1e − 15 ; % query coverage > 30), and more than 98% of them matched to monocot sequences, mainly corresponding to S. bicolor, S. italica, P. hallii and Z. mays (Fig. 2). Besides, more than 96,350 transcripts showed significant hits against the UniProt database (taxon identifier: 58,024, BLASTP, E-value: 1 x e − 5 ). In total, 66,617 transcripts were grouped into 755 Cellular Component-, 88,776 into 558 Biological Process-, and 15,620 into 147 Molecular Function-GO terms. Additional file 2 (A-C) shows the 30 most representative GO terms for each category. Besides, 43,275 transcripts were annotated into 168 Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, with carbon metabolism, biosynthesis of amino acids, spliceosome, endocytosis, mRNA surveillance and glycolysis/gluconeogenesis as the most represented routes (Additional file 2, D).

Apomixis vs. sexuality differential expression analyses
Based on the GTA, we analyzed the expression level of the complete set of sense and antisense transcripts by using the DESeq2 package. Comparisons included: i) a global differential expression analysis (GDEA), in which pairwise comparisons concerned the reproductive modes (apomictic vs. sexual) considering all developmental stages together (i.e., involving reads generated from all developmental stages), and ii) a stage-specific differential expression analysis (SSDEA) between reproductive modes at each developmental phase (i.e., involving reads generated from each particular developmental stage).

Gene ontology and KEGG pathway classification
Next, we established a Gene Ontology (GO) classification for all sense transcripts expressed at each developmental stage (Additional file 8). Based on this catalogue, we carried out a KEGG pathways prediction to identify specifically-modulated molecular routes (i.e., upregulated or downregulated in apomictic plants) (Additional file 10). Homologous recombination, endocytosis, thiamine metabolism, monobactam biosynthesis, lysine biosynthesis, protein export, RNA polymerase, photosynthesisantenna proteins, mismatch repair, mRNA surveillance pathway, fructose and mannose metabolism, various types of N − glycan biosynthesis, terpenoid backbone biosynthesis and phenylalanine metabolism pathways are regulated (by upregulation or repression) only at premeiosis, a clear turning point in which the onset of apospory initials occurs (Additional file 10). In the rest of the developmental stages these routes are similarly expressed in flowers of apomictic and sexual plants. Besides, mitogen-activated protein kinase (MAPK) signaling pathways, glycosphingolipid biosynthesis and galactose, pyrimidine, sphingolipid and biotin metabolisms are modulated exclusively during meiosis. Moreover, valine, leucine and isoleucine degradation, ABC transporters, circadian rhythms, autophagy, plantpathogen interaction, fatty acid metabolism, aliphatic and aromatic aminoacid biosynthesis and selenocompound metabolism are regulated exclusively at postmeiosis. Finally, protein processing in the endoplasmic reticulum, ribosome, citrate cycle, oxidative phosphorylation, phagosome, nitrogen metabolism, alanine, aspartate and glutamate metabolism, carbon fixation in photosynthetic organisms and fatty acid elongation are modulated only at anthesis, when parthenogenesis starts. Protein interaction predictions among members of these stage-specific routes analyzed by comparisons with the STRING database [62], revealed they integrate tight networks worth to get functionally explored (Fig. 4).
A similar GO analysis was conducted for DEATs (Additional files 9 and 10), but this time up-and downregulated entities were not distinguished due to the low number of candidates involved. We were able to identify sequences involved in mRNA surveillance and ubiquitinmediated proteolysis modulated only at premeiosis; others related to cyanoaminoacid metabolism modulated during postmeiosis; and pentose/glucuronate interconversions, vitamin B6 metabolism, flavonoid biosynthesis, ether lipid metabolism and fatty acid elongation occurring exclusively at anthesis. Interestingly, DEATs related to endocytosis and spliceosome resulted exclusively regulated at two developmental stages: meiosis and anthesis (Additional files 9 and 10).
Then, Gene Ontology (GO) Enrichment analyses were performed on a subset of 5268 DETs which were common to the four developmental stages analyzed (i.e., transcripts that were found differentially expressed at all stages of sexual and apomictic developments) (Additional file 11). The most represented GO terms regarding Cellular Component (CC) were ribonucleoprotein complex, nuclear lumen, vacuolar membrane, vacuolar part and chloroplast thylakoid. The most represented Biological Process (BP) terms were protein localization to organelle, the establishment of protein localization to organelle, nucleotide biosynthetic process, nucleoside phosphate biosynthetic process and pyruvate metabolic processes. The Molecular Function (MF) main classes were structural molecule activity, structural constituent of ribosome, adenyl nucleotide binding, adenyl ribonucleotide binding and ATPase activity. Moreover, the main KEGG-predicted pathways were ribosome, carbon metabolism and spliceosome (Additional file 11).
The same study was applied to common DEATs. In this case, out of 32 total transcripts, only a few could be assigned to GO terms corresponding to CC (nucleus), MF (zinc ion binding, stearoyl-CoA9-desaturase activity, ATP binding and amylase activity) and BP (mature ribosome assembly) (not shown). Although the limited number of DEATs did not allow a proper KEGG pathway evaluation, we could identify that fatty acid metabolism, MAPK signaling pathway, starch/sucrose metabolism and plant hormone signal transduction pathways were changed (not shown).

Transcriptome dynamics
Matrixes representing the normalized raw read counts for each transcript in each one of the libraries are provided in Additional file 12 (S1 and S2 for sense and antisense transcripts, respectively). The advantage of counting with samples representing different developmental stages allowed us to perform a cluster analysis to identify groups of transcripts with similar expression patterns (Fig. 5). Normalized counts were used to execute a hierarchical clustering using a simple euclidean distance metric and a complete linkage method, to find some structure in our transcript expression trends and consequently partition our transcripts into different groups. To enable the analysis, subsets of transcripts were used, corresponding to: 1) all DETs and DEATs; 2) transcripts that were differentially expressed at all the developmental stages (common DETs and DEATs). 3) stage-specific DETs and DEATs (premeiosis, meiosis, postmeiosis and anthesis). The number of transcripts within each cluster is displayed in Table 2.
Clustering analysis split all DETs into six clusters displaying distinct expression patterns (Fig. 5a). Cluster 1 included transcripts up-regulated in the sexual samples at all stages (Fig. 5a). Clusters 2 and 3 showed transcripts up-regulated in apomictic samples at all stages, while clusters 4, 5 and 6 displayed variable divergent expression patterns (Fig. 5a). Particularly, Cluster 6 includes a group of transcripts showing expression heterochronicity, with entities overexpressed at postmeiosis/anthesis in apomictic plants but at premeiosis/meiosis in sexual ones (Fig.  5a). In the case of all DEATs, two major clusters were established (Fig. 5c). In Cluster 1, members show variable expression in sexual samples in the course of development, while in apomictic ones the expression seems lower or null (Fig. 5c). The opposite behaviour was observed in Cluster 2 (Fig. 5c).
Interestingly, the analysis of the common DETs (transcripts that were differentially expressed at all developmental stages) revealed a striking contrast between sample types (Fig. 5b). Clusters 1, 3 and 4, showed upregulation in apomictic samples across the four developmental stages, while the same transcripts were consistently repressed in sexual ones (Fig. 5b). The opposite occurred in Clusters 2, 5 and 6 ( Fig. 5b). In the case of common DEATs, Cluster 1 showed high expression increasing in the course of development in apomictic samples, but low modulated expression in sexual ones (Fig. 5d). The opposite occurred in Cluster 2 (Fig. 5d).
A similar study was conducted with the stage-specific DETs and DEATs. Clustering graphs are presented in Additional file 12 (S3 and S4, respectively). In all cases (premeiosis, meiosis, postmeiosis and anthesis), transcripts were organized into 6 different clusters with contrasting behaviour. Expression heterochronicity is clearly visible in clusters like DETs Postmeiosis 5 and DETs Anthesis 5.

Differential expression of transcription factors
To investigate the nature, the level and the expression timing of transcription factors (TFs) detected in the apomictic and sexual samples, we contrasted our Illumina transcriptomes against the Plant Transcription Factor Database (http://planttfdb.cbi.pku.edu.cn). Based on the BLASTx top hits, we identified 63,076 (31.67%) transcripts highly similar to TFs, corresponding to 60 different families

Differential expression of transcripts associated with plant hormones
A search for hormone-related transcripts expressed during the apomictic and sexual reproductive developments was performed by testing our transcriptomes against the Arabidopsis hormone-related protein database (http://hormones. psc.riken.jp/pathway_ja.html). A BLASTx analysis showed 3781 top hits associated with plant hormones and related compounds, including 714 related with auxin, 595 with jasmonate, 592 with cytokinin, 520 with abscisic acid, 517 with gibberellin, 487 with brassinosteroids, 245 with salicylic acid and 111 with ethylene (Additional file 15). Of them, 56 (related to auxin), 70 (related to jasmonate), 39 (related to cytokinin), 42 (related to abscisic acid), 46 (related to gibberellin), 36 (related to brassinosteroids), 21 (related to salicylic acid) and 5 (related to ethylene) transcripts resulted differentially expressed between the apomictic and sexual samples (DE Hormones) (Additional file 15, Fig. 6).

Apomixis candidates survey
In previous work, several genetically-linked, differentiallyexpressed or functionally-associated apomixis candidate genes were reported in P. notatum and other related species [36,50,58,64]. Additional file 16 shows some of these candidate genes (with their respective identifiers) and 49 DETs displaying significant similarity them, as well as the stages at which differential expression is detected, the Log 2 FC and padj values for each developmental stage and the annotation of the sequence. Ocassionally, different DETs showed the same annotation but displayed contrasting expression profiles, pointing to the existence of transcripts variants of the same gene expressing at distinct reproductive stages. Ten (10) DETs showed different expression levels at all developmental stages, three were specific for premeiosis, five for meiosis, three for postmeiosis and nine for anthesis (Additional file 16). The rest revealed differential expression levels at more than one stage. The identities of DETs include, among others, KIN-14P [28], ENHANCED DISEASE RESISTANCE 2 APOS-TART1/2 [65], FIE [66], LORELEI-like (N20) [67], SERK1/ [68], TGS1 [69], ORC3 [39], MAP3K (N46) [41], GID1 [70] (Additional file 16). Moreover, an examination of the candidate expression levels revealed some variation between stages and/or reproductive modes. For example, transcript TRpn_185717 which codifies for an ARGO-NAUTE 104 protein [71], showed a constant higher level of expression in the apomictic genotype during all the developmental stages, while its expression decreased from premeiosis to anthesis in the sexual one (Fig. 7a). Transcript TRpn_57024, with high similarity to a gibberellin receptor (GID1), exhibited an insignificant expression in the apomictic genotype at all developmental stages with a peak at anthesis, while it was steadily expressed at low levels in sexual genotypes (Fig. 7b). Two other genes, ORC3b (TRpn_33086/TRpn_96407) and APOSTART2 (TRpn_175767 /TRpn_52614) showed the same expression patterns that exhibited in the in vivo experiments reported in P. simplex and P. pratensis (Fig. 7c-f) [39,65].

Discussion
Apomixis in Paspalum notatum is a genetically heritable trait controlled by a single complex-locus (ACL), which shows non-mendelian segregation against apomixis and a strong restriction of recombination involving around 36 M[bp [72][73][74]. The ACL is localized in a chromosome fragment syntenic to subtelomeric segments of rice chromosomes 2 and 12 and maize chromosomes 1 and 3 [55,75,76]. This genomic region displays typical heterochromatin features, like the presence of repetitive elements, gene degeneration and cytosine hypermethylation [75,77]. The complex genomic topography together with the unavailability of a reference genome have seriously compromised the identification of genes controlling the trait by positional map-based approaches. As a consequence, transcriptomic surveys became essential tools for identifying apomixis-related genes in the species [48,50].
In the early 2000s, PCR-based methods were used to predict several P. notatum transcripts associated with the occurrence of agamospermy [28,29]. These preliminary amplification analyses, although rather limited, served well to the selection of candidates and the establishment of further functional characterizations that led, finally, to the identification of major components of the apomixis pathway in this species. In particular, the MAP3K-encoding QGJ (QUI-GON JINN) transcript, initially identified by Laspina et al. (2008) as clone N46, was validated in sexual and aposporous plants through RNA in situ hybridization experiments [41] and functionally classified as an inducer of aposporous embryo sacs (AESs) formation in Paspalum RNAi lines [41]. Moreover, TGS1 (TRIMETHYLGUANOSINE SYNTHA SE 1) [29] was also confirmed to express differentially in sexual and aposporous ovules by RNA in situ hybridization experiments [69] and classified as an AESs repressor in Paspalum antisense lines [78]. The functional association between the apospory induction/repression pathways directed by QGJ and TGS1 is currently under investigation.
The development of the first P. notatum floral transcriptomes of reference for apomictic and sexual genotypes widened the opportunities for identifying reproductive genes and provided a whole catalogue of full-length molecules related to traits of agronomical interest [31]. However, the scarcity of reliable statistical information due to the absence of replicates was a serious challenge to network interaction predictions and functional analysis preparation [31]. Moreover, since the reference transcriptomes were produced from bulked samples representing different developmental stages, the temporal variation of expression during development remained unexplored at a wide-genome level. Indeed, stage-specific qRT-PCR analyses were carried out for some of the candidates in several genotypes of contrasting reproductive modes, in order to evaluate the chronological evolution of the expression in the course of development. Finally, no information on the antisense/ sense orientation of the DETs was available, a fact that partially constrained the elucidation of the mechanisms involved in the sexuality/apomixis transition. In particular, since the differential expression of several retrotransposons [79], long non-coding RNAs (lncRNAs) [40] and pseudogenes [39] was associated with the occurrence of apomixis in Paspalum, determining the sense/ antisense nature of these groups of transcripts is central to disclose their possible role on molecular mechanisms regulating reproduction.
The set of sequences presented here considerably increases the number of transcripts available from the Roche 454/FLX + reference library [31] and offers information on the expression dynamics for each sequence in the course of development. The GDEA analysis showed 13,205 (68%) and 6147 (32%) transcripts up-and down-regulated during apomixis, providing full-length sequences for most of them and pointing to a greater complexity for the molecular control of agamospermy with respect to sexuality. Meanwhile, the SSDEA study revealed a higher number of DETs than GDEA, a clear indication of Fig. 7 Expression patterns of apomixis-related transcripts across reproductive developmental stages from premeiosis to anthesis. a Argonaute 104 (TRpn_185717). b Gibberellin receptor GID1 (TRpn_57024), c, d ORC3b (TRpn_33086, TRpn_96407). e, f APOSTART2 (TRpn_175767, TRpn_52614). Repd: reproductive mode. Red lines: apomictic libraries. Blue lines: sexual libraries expression heterochronicity for a considerable number of transcripts (that is because those genes showing upregulation at different developmental stages in apomictic and sexual plants are classified as 'non-differential' in the GDEA). Particularly, a large number of DETs and DEATs was detected at anthesis, which reveals an abundance of exclusive pathways at this stage. Regarding this, chronological differences in the onset of the sexual and apomictic embryonic periods should be taken into account: in apomictic plants, parthenogenesis of the 2n (non-reduced) egg cell starts at the end of megagametogenesis, while in sexual plants, the legitimate n (reduced) egg cell remains quiescent until fertilization [78]. Then, at least a part of the differentially expressed transcripts detected at anthesis might be related alternatively to the induction of parthenogenesis or to embryo development, since, at anthesis, embryos are forming in apomictic plants only. Besides, other genes like those related to pollen-stigma interactions or pollen discharge into the egg cell could also be integrating this particular group of DETs.
We found 11,417 antisense transcripts expressed in flowers, and hundreds of them were differentially represented in sexual and apomictic genotypes. These results confirmed the evidence reported in previous articles, in which both sense and antisense transcripts were detected in apomictic and sexual genotypes at variable representation levels or altered localization [39,41,68]. Several of these transcripts, like those homologous to SERK2, MAPK3 and ORC3 [39,41,68] where found here represented by antisense strands and thus, a regulatory mechanism based on complementary hybridization might be associated with them. The confirmation of such a modulation process will require further functional analysis. Interestingly, a methylation-mediated silencing mechanism was reported to control parthenogenesis in other species of the Paspalum genus [77]. We detected a higher number of DEATs occurring at anthesis, in comparison with other developmental stages. Moreover, the anthesis DEATs have a tendency (59%) to appear downregulated in apomictic plants. Our results indicate that: 1) antisense-mediated regulatory mechanisms might be particularly active at anthesis in sexual plants; and 2) silencing of a considerable number of DEATs occurs in apomictic plants at this particular stage. It will be interesting to investigate the influence of genome methylation in the down-regulation of DEATs during apomixis, and its consequences on the representation of DETs and the emergence of autonomous embryos in the absence of fertilization.
Two of the main ontology classes previously associated with apomixis in Paspalum corresponded to transcription factors [31,50] and hormones [31,50,80]. Here, we identified a number of transcriptional regulators showing differential expression in reproductive organs of apomictic and sexual plants. Some of them, like MYB family members, show heterochronic upregulation in the apomictic and sexual libraries, with expression reaching peaks at different developmental stages. MYB proteins share the conserved MYB DNA-binding domain that is crucial to the control of proliferation and differentiation in several cell types [81]. Besides, upregulation in the apomictic libraries was observed for members of the bHLH, ERF, WRKY, B3, ARF and AP2 family, all of them previously related to plant development through cell division, proliferation and differentiation control [82]. We also identified numerous members of different hormonal pathways, especially the auxin and jasmonate routes, that show differential regulation during apomixis, a fact that had been anticipated when analyzing the representation of miRNAs in reproductive organs of sexual and apomictic plants [80]. Currently several members of these pathways are under functional characterization to determine its role in apomixis development. Besides, the expression characterization for several transcripts previously associated with apomixis provided information on the specific stages these transcripts are modulated and a proof of concept that the use of the catalogue presented here might contribute to a thorough comparison of both reproductive pathways.
The Illumina sequence database reported here represents a detailed chronological characterization of the sense/antisense gene expression landscape in reproductive organs of sexual and apomictic counterparts of the same species. The derived information could be of use in Paspalum research programs dealing with gene expression during sexual and asexual seed formation, as well as the molecular breeding of apomixis. Moreover, it will allow the identification of apomixis candidate genes, which could be further characterized in expression and function in other apomictic species. The reported information could not be completely systematized in a single scientific article, since a number of different analyses can be conducted, depending on the need of each specific research project. However, the main use for this tool no doubt will be the comprehensive identification of candidate genes that can be used as molecular markers in apomixis-based breeding programs or to induce asexuality from sexual genetic backgrounds through genetic engineering. In particular, it will allow for rapid discrimination of some of the sequences controlling apomeiosis and parthenogenesis, due to its potential to expose differential expression at specific stages. Such discrimination was impracticable when using the formerly-produced apomictic vs. sexual transcriptomic databases, since their construction invariably involved only one reproductive stage or a bulk of stages. Currently, the Paspalum genome is in process of sequencing and assembly in our laboratory and the apospory controlling locus (ACL) is being identified by positioning markers that were fully linked to the trait in former genetic mapping experiments. In this context, mapping the stage-specific candidate sequences exposed here onto the Paspalum ACL will help to identify the genomic controllers of apomixis, while RNA in situ hybridization will reveal the precise site of expression in reproductive tissues. Eventually, functional analysis will disclose the reproductive phenotypes that can be induced after up-or down-regulation in precise cell types pointed by the in situ analysis. Moreover, genetic engineering will allow the harnessing of these candidates to reproduce the desired reproductive phenotypes in species of interest, using the appropriate promoters. On another note, an additional application of this database, among many others, will be aimed at clarifying the functionality of the expression originated from the heterochromatic nonrecombinant ACL, since mapping the antisense transcripts identified here onto the Paspalum genome will reveal which of these transcripts emerge from this particular region.
Difficulties involved in the elucidation of the molecular control of reproduction in P. notatum range from the usual complexities of all reproductive systems (the molecular intricacy of the routes involved, the high temporal variation rate, the involvement of cell-specific expression patterns, among others) to those derived from the particular nature of apomictic species, like the evolutionary and physical characteristics of the ACL and the involvement of poorly characterized polyploid heterozygous genomes. Despite all these drawbacks, during the last years, an unprecedented research effort led to detailed characterization of the operational molecular routes for both P. notatum reproductive modes, through the establishment of publicly available solid reference transcriptomes [31] and replicate sRNA libraries [80]. Here, we are deepening the characterization of the molecular transcriptional landscape operational in sexual and apomictic plants, by providing a chronological, high yield, orientation-sensitive transcript database covering all reproductive stages from premeiosis to anthesis. We hope that this contribution will serve as a basis to promote future research on the functional mechanisms controlling agamospermy in plants and as a valuable resource for those plant breeders who are focused on the introduction of apomixis technology into their cultivar improvement programs.

Conclusions
Here we introduce a complete sense/antisense gene expression catalogue from florets of apomictic and sexual P. notatum, involving four subsequent reproductive developmental stages, from premeiosis to anthesis. This comprehensive sequence collection quantitatively reveals apomixis vs. sexual heterochronic expression and sense/ antisense mediated regulation. In particular, contrasting transcriptional and hormonal control was detected and thoroughly characterized. Our analysis exposed a considerable alteration of sense/antisense expression occurring at the onset of parthenogenesis. The experimental approach used in this work established a set of differentially expressed sequences well beyond the former group of candidates detected in Paspalum, which even discriminated the sequence orientation, giving important clues on antisense-mediated transcriptional and posttranscriptional regulation. This dataset will be applied to a more efficient selection of apomixis candidate genes, contributing to the future development of molecular tools for harnessing the trait.

RNA isolation
Inflorescences at different stages of the reproductive development were collected from both the Q4117 (apomictic) and C4-4x (sexual) plants. The classification of the reproductive stages was carried out following the protocol described by Laspina et al. [29], by analyzing in parallel both the mega-and microsporogenesis as well as the mega-and microgametogenesis processes. Four stages were considered: (I/II) premeiosis: megaspore mother cells (MMCs) and apospory initials (AIs) are visible in ovules, while tetrads start to appear in anthers; (III) meiosis: uninucleate pollen and female meiosis occurs in the sexual genotype; (IV-VI) postmeiosis: uninucleated/binucleate pollen, and first division of megagametogenesis; and (IV) anthesis: binucleate pollen and mature embryo sacs [29]. Total floral RNA was extracted with the SV RNA Total Isolation Kit (Promega) and quantified using the Quant-iT RiboGreen RNA Reagent and Kit (Invitrogen). Three replicates were established, corresponding to different floral RNA extractions from the same genotypes (Q4117 and C4-4x). The RNA quality was evaluated with RNA 6000 Pico-Chip (Agilent Bioanalyzer 2100).

Library preparation and Illumina HiSeq sequencing
Library preparation and sequencing experiments were carried out at Instituto de Agrobiotecnología de Rosario (INDEAR), Rosario, Argentina. Libraries were built using the TruSeq® Stranded mRNA kit (Illumina) starting from 1 μg of total RNA. Procedures for mRNA purification (using oligo-dT hybridization), RNA fragmentation, double-stranded cDNA synthesis, end-adenylation, ligation of adapters and enrichment (amplification) of library fragments were performed following the protocol described in TruSeq® Stranded mRNA Illumina (October 2017). The library quality was checked with the DNA 1000 Kit (Agilent Technologies), using 1 ul of each preparation in a 2100 Bioanalyzer. Libraries resulted in double-stranded DNA fragments with an average size of 260 bp. Three biological samples were processed for each one of the developmental stages (premeiosis, meiosis, postmeiosis, anthesis), for both reproductive modes. Thus, a total of 24 (3 × 4 × 2) libraries were constructed. Before sequencing, the individual libraries were quantified by qPCR (Light Cycler 480 Roche) and normalized to a concentration of 3 nM. One equimolar pool of all libraries was prepared and quantified by qPCR (Light Cycler 480 Roche) using the Qiagen Library Quantification Kit. The pool was used for the generation of clusters in the sequencing cell. A sequencing run was performed by generating paired ends (PE) 2 × 100 bp reads in a HiSeq-Illumina device.

Bioinformatics methods
Raw reads were de-multiplexed and checked for quality using the FastQC software (http://www.bioinformatics. babraham.ac.uk/projects/fastqc/). Adaptors, duplicated sequences, ambiguous reads, and low-quality reads were r e m o v e d b y u s i n g T r i m G a l o r e ( h t t p : / / w w w . bioinformatics.babraham.ac.uk/projects/trim_galore/). The high-quality reads (QC > 30) were used for assembling the transcriptomes with Trinity v2.0.2 [60,61]. The available Roche 454/FLX + P. notatum reproductive transcriptome generated by Ortiz et al. [31] was initially used for a reference-guided assembly with the parameters:"--SS_lib_type RF --normalize_by_read_set --min_ contig_length 400", and --genome_guided_bam. Then, a second de novo assembly (without a reference) was carried out using the default parameters. The four assemblies (apo_over reference, sex_over reference, apo_ without reference and sex_without reference) were combined in one file and the non-redundant transcripts were selected using CD-HIT [85]. The quality of the assemblies was measured using QUAST [86]. The raw reads were mapped to the global assembly using Bowtie2 (v 2.3.2.0) [87] and TopHat (v2.1.1) [88]. The transcriptome coverage was evaluated using Benchmarking Universal Single-Copy Orthologs (BUSCO, version 3.0.2) [89,90] with the following commands: "Python run_ BUSCO.py -i sequence_file -o output_name -l lineage -m tran", "Python generate_plot.py -wd working directory" and the "liliopsida_odb10" dataset. The coding competence for all expressed transcript sequences was predicted using the TransDecoder software (https:// github.com/TransDecoder/TransDecoder.wiki.git) with the default parameters (−m 100). For differential expression analysis, RNA-seq reads were analyzed with the Kallisto v.0.44.0 software to determine transcript counts and abundances. The libraries were normalized by size and low-count transcripts were filtered out (< 3 in the three replicas of each library). Differential expression analysis and the corresponding p-values were estimated using the Bioconductor software package DESeq2 [91]. The p-values attained by the Wald test were corrected for multiple testing using the Benjamini-Hochberg method. The adjusted p-values (DESeq2 padj or FDR) thresholds for considering transcripts as DETs/DEATs were < 0.001 and < 0.05 for sense and antisense transcripts, respectively. Moreover, analyses of DETs were restricted to those showing an absolute value of Log 2 FC > |3|. Comparisons of gene expression between modes of reproduction were carried out considering all stages of development (global comparison) or each developmental stage (stage-specific comparison). Venn diagrams were created by using the jvenn online tool/ software (http://jvenn.toulouse.inra.fr/app/example.html) [92]. The transcriptome dynamics was analyzed with R: normalized counts were used to execute a hierarchical clustering using a simple euclidean distance metric and a complete linkage method.

Availability of data and materials
The datasets generated and/or analyzed during the current study are available in the NCBI SRA repository at the following SRA accession: PRJNA511813. This Transcriptome Shotgun Assembly project has been deposited at DDBJ/EMBL/GenBank under the accession GIUR00000000. The version described in this paper is the first version, GIUR01000000. Other dataset(s) supporting the conclusions of this article are included within the article and its additional files. The plant materials used in this study belong to the living germplasm collection of Instituto de Botánica del Nordeste Ethics approval and consent to participate Not applicable.

Consent for publication
Not applicable.