Ovule Transcriptome Analysis Discloses Deregulation of Genes and Pathways in Sexual and Apomictic Limonium Species (Plumbaginaceae)

The genus Limonium Mill. (sea lavenders) includes species with sexual and apomixis reproductive strategies, although the genes involved in these processes are unknown. To explore the mechanisms beyond these reproduction modes, transcriptome profiling of sexual, male sterile, and facultative apomictic species was carried out using ovules from different developmental stages. In total, 15,166 unigenes were found to be differentially expressed with apomictic vs. sexual reproduction, of which 4275 were uniquely annotated using an Arabidopsis thaliana database, with different regulations according to each stage and/or species compared. Gene ontology (GO) enrichment analysis indicated that genes related to tubulin, actin, the ubiquitin degradation process, reactive oxygen species scavenging, hormone signaling such as the ethylene signaling pathway and gibberellic acid-dependent signal, and transcription factors were found among differentially expressed genes (DEGs) between apomictic and sexual plants. We found that 24% of uniquely annotated DEGs were likely to be implicated in flower development, male sterility, pollen formation, pollen-stigma interactions, and pollen tube formation. The present study identifies candidate genes that are highly associated with distinct reproductive modes and sheds light on the molecular mechanisms of apomixis expression in Limonium sp.


Introduction
Apomixis (agamospermy), the asexual seed production found in less than 1% of angiosperms [1], can be either independent of or dependent on pollination [2,3]. Most natural apomicts produce meiotic-reduced pollen involved in the fertilization of the polar nuclei in the embryo sac (pseudogamy) [2,3], but others reproduce independently of pollination for the initiation of the embryo or endosperm formation (autonomous apomicts) [4,5].
Different from sexual reproduction, apomixis is characterized by alterations in the developmental program during the formation and development of the female germline [2,6].
Apomicts can reproduce via gametophytic apomixis, involving the formation of an unreduced embryo sac (apomeiosis) that gives rise to a parthenogenetic embryo and a functional endosperm without the 2 maternal: 1 parental genome ratio [2,6]. The unreduced gametophytes can develop via restitutional meiosis or mitotic division (diplospory) or by a somatic, unreduced cell of the nucellus, which develops into an embryonic sac (apospory) [2,6]. In sporophytic apomicts (adventitious embryonies), the embryos originate from somatic tissues of the nucellus and/or integument cells. The formation of these apomictic embryos usually occurs in parallel with the formation of sexual embryos [7].
The emergence of apomixis in natural systems has been a long-standing topic of debate. It was hypothesized that the different types of apomixis are caused by different mutations that destabilize meiosis (megasporogenesis), the gametophyte (embryo sac), and egg formation [2,6]. Loci genetically linked to components of apomixis have been identified in various species, and sequencing of these loci has revealed several genes with the potential to play critical roles in apomixis [17,18]. It was hypothesized that apomixis could be caused by asynchronously expressed germline genes in the ovules of certain hybrids [19]. Transcriptome comparisons show that the genetic control of apomixis in gametophytic apomicts has been related to a wide range of mechanisms regulating gene expression, including protein degradation, transcription, cell cycle control, stress response, hormonal pathways, cell-to-cell signaling, and epigenetic mechanisms [18,20,21]. Several common genes found to be differentially expressed in multiple stages of apomictic and sexual seed production support the view that sexual and apomictic reproduction are closely related developmental pathways [22]. Recent studies present substantial evidence in support of a polyphenic condition for meiosis and determine that polyphenic shifts from apomeiosis to meiosis and vice versa are regulated by metabolic states [23].
[Comparisons: uppercase refers to test samples, lowercase refers to control samples, and numbers refer to respective stages].

Ovule Extraction
Flower buds at distinct developmental stages were sampled prior to anthesis according to their size (between 2 and 5 mm) based on cytoembryological observations as in [13], at standardized times (between 9:00 a.m. and 12:00 p.m.). The ovules were selected with respect to the timing of apomeiosis (Stage 1-S1), megagametogenesis (embryo sac, Stage 2-S2), parthenogenesis, and endosperm formation (Stages 3/4-S3/S4), detailed in [13]. Dissection of ovules was performed in a sterile laminar air flow cabinet under a stereoscopic microscope (Stemi 2000-C, Zeiss) with the aid of tweezers and pencil-point spinal needles (Transmed). In total, 280 ovules were extracted from each ovary, containing one ovule each, and about ten to twenty ovules per stage were isolated and placed in a sterile Petri dish with B5 medium [39] to maintain hydration before RNA extraction. This procedure generated a total of 18 samples, including nine samples of ovules from apomictic plants (four biological replicates in stage S1 and five biological replicates in stage S2) and three samples from each of the remaining species (sexual: two species x three stages-S1, S2, and S3/S4; facultative apomictic: three biological replicates in stage S3/S4) ( Table 1).

Total RNA Extraction and Library Preparation
Total RNA was extracted from all samples using the Spectrum™ Plant Total RNA Kit (Sigma-Aldrich, St. Louis, MO, USA). Nevertheless, some modifications were required for obtaining high-quality RNA, as detailed. Samples were collected in 150 µL of lysis solution and grounded with the help of a micropestle in the microfuge tube. Then, another 300 µL of lysis solution supplemented with β-mercaptoethanol was added to the tube, which was vortexed vigorously and incubated at 56 • C for 5 min. The lysate was centrifuged at 14,000 rpm for 2 min, and then the supernatant was transferred into a filter column and centrifuged for 1 min at 14,000 rpm. The clarified flow-through lysate was transferred to a new tube, and 250 µL of Binding Solution was added. The mixture was applied to a binding column and centrifuged for 1 min at 14,000 rpm. The remaining steps followed the manufacturer's instructions. The quantity of RNA was determined using a BioDrop cuvette (BioDrop, Cambridge, UK) and electrophoresis on a 1% agarose gel. The RNA integrity number (RIN) was determined using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) and ranged from 9.19 to 9.45. The messenger RNA (mRNA) libraries were constructed with the Illumina "TruSeq Stranded mRNA Sample Preparation kit" (Illumina, San Diego, CA, USA) and sequenced on an Illumina NovaSeq6000 2× 100 bp at Macrogen facilities (Macrogen, Geumcheon-gu, Seoul, Korea).

Processing, Mapping, and Quantification of Illumina Reads
All raw reads have been deposited in the NCBI Sequence Read Archive (SRA), Bio-Project accession PRJNA752506. Quality control of the raw reads, including contaminants survey, was performed using FastQC version 0.11.9 [40] and FastQ Screen version 0.14.0 [41], which ran against the genomes of their default pre-indexed species and adaptors. Then, since all raw reads presented a quality base score over 36, Trimmomatic version 0.39 [42] was used to eliminate adaptors and filter reads of length below 36 base pairs (bp). A de novo transcriptome assembly was performed using Trinity version 2.11.0 [43], in which cleaned reads from all samples were combined to generate one global assembly since this software has shown consistent performance and has a high read alignment rate [44]. The assembly was assessed for completeness using BUSCO version 5 [45] through gVolante2 [46]. After alignment against the transcriptome using Bowtie2 aligner version 2.3.5 [47], sequences were quantified at gene-level expression with RSEM version 1.3.3 [48] through the Trinity pipeline. A principal component analysis (PCA) was performed to survey the relatedness of normalized gene counts using the function plotPCA in R Studio version 4.0.2 [49].

Differentially Expressed Genes Detection
To study significant differences between apomictic and sexual plants, differential expression analysis was performed with edgeR version 3.30.3 [50], which is a flexible empirical Bayes approach that uses weighted likelihood methods to estimate gene-specific variation even with very few or no replicates [51]. Overall, when studying differences between reproductive strategies, apomictic plants were set as the samples to test, while sexual and facultative apomictic plants were set as the controls, according to each comparison ( Table 1). As such, up-regulated DEGs are more expressed in apomictic than in sexual plants, while down-regulated DEGs are less expressed in apomictic and more expressed in sexual plants.
Genes with a normalized |log2 fold change (log2 FC)| > 2 were defined as differentially expressed and used in the downstream analysis. In the comparison between apomictic and facultative apomictic plants, in which all samples have at least 3 replicates, DEGs were previously filtered by p < 0.01. Venn diagrams were used to plot DEGs between different comparisons through matplotlib version 3.3.3 [52] in Python version 3.9.0 (Python Software Foundation 2020). Additionally, DEGs commonly triggered by more than one comparison were searched for opposite regulation.

Functional Annotation
Functional annotation of DEGs was performed with the Basic Local Alignment Search Tool (BLAST) version 2.10.1 command-line tool from the NCBI C++ Toolkit (National Center for Biotechnology Information 2020). Blastx was used to map DEGs to A. thaliana homologs against a local Swissprot database, filtering gene hits by a maximum E-value of 1.0E −3 and a minimum identity of 40%. [53]. Then, to avoid duplicated results, DEGs annotated to the same A. thaliana homolog were filtered by identities and sequence length, keeping the transcripts with the highest values.
Housekeeping genes (HKGs) are typically required for the maintenance of basal cellular functions that are essential for its existence, regardless of their specific role in the tissue or organism. Since these genes are usually highly conserved, genes stably transcribed in all comparisons were filtered. Additionally, DEGs were searched for the most common housekeeping homolog genes in A. thaliana to study their potential role in reproduction in Limonium plants.
In addition, DEGs were searched for sRNA biogenesis, which is known to play pivotal roles in reproductive development [54,55], and oxidative stress-related genes, which negatively affect reproductive development in plants [56]. DEGs associated with tryptophan and ethylene metabolism were investigated, which are associated with the biosynthesis and regulation of the phytohormone auxin, a vital component of plant reproduction since it regulates both male and female reproductive organs [57,58]. Moreover, DEGs related to aminoacyl-tRNA metabolism and lysine degradation, which are respectively essential to produce ribosomes and proteins, and epigenetic processes through DNA methylation [59,60] were searched.

Enrichment Analysis
Uniquely annotated DEGs were characterized with GO terms using the REST API on the UniprotKB website (The Uniprot Consortium, 2019). Finally, GO enrichment analyses were applied to log2FC-ordered lists of DEGs through an over-representation analysis (ORA) using the g:GOSt functional profiling tool from the gProfiler website [67], with the g:SCS tailored algorithm under FDR < 0.01, using a predefined A. thaliana custom background including only genes expressed by the samples in analysis. Enrichment results were summarized using REVIGO [68] through the removal of redundant GO terms with allowed similarity = 0.5 and then plotted with the R ggplot2 version 3.3.2 library [69]. To better understand the differences between the initial and final stages of both sexual and apomictic plants, ORA results were filtered by specificity to a particular stage.

Overall Differential Expression Analysis
The highest number of annotated DEGs (3544) was found in M2o4, followed by M2a4 (1174) ( Table 1). In intraspecies comparisons, L. ovalifolium showed the lowest number of DEGs in O2o1 (693), followed by O4o2 (1650), with the major differences being observed in O4o1 (2067; Supplementary Figures S2B and S3C,D). Conversely, L. auriculifolium showed the lowest differences in A4a1 (351), followed by A2A1 (796), thus presenting the major differences in A4a2 (1346; Supplementary Figures S2A and S3A,B). Notably, DEGs showed a proportion of 24% to 35% unannotated genes across samples (Table 1). On the other hand, M1o1 and M1o2 showed just slightly more DEGs than M1a2, while M2m1 presented more DEGs than M2d4. When comparing the two sexual plants, there was an increasing number of DEGs with the progression of stages, which varied from 616 to 1611 (Table 1; Supplementary Figure S2C).

Gene ID
Analyzing DEGs associated with sRNA biogenesis showed an up-regulation of various genes like AGO1, AGO4, AGO5, AGO7, AGO8, AGO9, DML2, and DNMT2 in M2o4. Nevertheless, in the remaining comparisons between apomictic and sexual ovules, while there was a down-regulation of AGO5 and AGO10 in apomictic ovules, DML2 was upregulated in the same plants (Supplementary Figure S5B).

General GO Enrichment
DEGs between apomictic and sexual ovules in early stages were found to be enriched in seven main ancestor terms ( Figure 4A   A.s thaliana, the most similar homolog of each differentially expressed gene (DEG) was mapped to the respective functional annotation, and enriched terms were summarized using REVIGO. Significantly (FDR < 0.01), gene ontology (GO) and biological processes (BP) terms are among DEGs from (A) apomictic in S1 relative to sexual ovules in S1 (M1a1 and M1o1) or S2 (M1a2 and M1o2), and from (B) apomictic in S2 relative to sexual ovules in S3/S4 (M2a4 and M2o4), and facultative apomictic in S4 (M2d4). The dot's size indicates the number of DEGs annotated with each term (counts), and the color shows the differential expression (red: up-regulated; blue: down-regulated). Enriched terms are grouped by their respective ancestors (ontology level 2).
was enriched among up-DEGs in M2a4. Metabolic process and response to stimulus terms were mainly enriched in up-DEGs from apomictic ovules relative to both sexual species ( Figure 4B).

Discussion
An increasing number of molecular studies have identified several candidate genes implicated in the shift from sexual to apomixis reproduction [18,21,22,70]. In different species, apomixis has been found to arise due to the action or deregulation of different genes associated with the normal sexual pathway [23,[71][72][73]. Nevertheless, it is still not fully understood how these genes alter reproductive pathways to establish apomixis.
In most apomictic wild species. Implementing omics approaches can be particularly challenging as complete genomic sequences are not available and, therefore, genome annotation information is not available. Additionally, obtaining plant material for transcriptomic studies can be an experimentally challenging task in Limonium since each plant presents a single ovary, enclosed in a calyx and inner, medium, and outer bracts that yield just a single basal ovule [13]. In the current study, we performed a comparative transcriptome analysis between sexual and asexual plants and identified candidate genes that are specifically or differentially expressed between reproductive modes and among stages of ovule development. This approach allowed us to disclose differential regulation of both HKGs as well as genes specifically involved in flower development, male sterility, and pollen recognition, besides major pathways potentially central to apomixis, including protein degradation, transcription, stress response, hormonal signaling, signal transduction, and epigenetic regulation.

Differential Regulation of HKG and Metabolic Pathways in Sexual and Apomictic Plants
In this study, the total number of expressed unigenes among samples was higher in ovules from apomictic plants than in those from sexual plants, together with a differential regulation of genes, particularly in the later stages of ovule development (Table 1). Previous studies between sexual and asexual plants provided support for the deregulation of reproductive pathways, including HKGs in, e.g., Boechera holboellii complex [72], Brachiaria [74], Cenchrus ciliaris [75], and Ranunculus [73], among others. In this study, although many homolog genes in Arabidopsis were stably expressed in Limonium, such as ACT domaincontaining proteins. Cytosolic Fe-S cluster assembly factors NBP35, NADH dehydrogenase [ubiquinone] flavoprotein 2, phosphoglycerate kinase 1, polyubiquitin 3, TATA-box binding protein 2, and other HKGs were found to be differentially expressed. These include genes related to the ubiquitin degradation process, such as ubiquitin-conjugating enzymes, tubulin, actin, and elongation factor-1 α as found in other sexual and apomictic plants' complexes above referred. Therefore, some of the HKGs identified (Table 3) in our study can be potentially used as reference genes to be validated in future quantitative gene expression studies using different developmental stages of specific tissue types or different reproductive modes.
A differential representation of DEGs in Limonium sexual and apomictic plants associated with the oxidative stress response was found. In apomictic plants, some of these DEGs (Supplementary Figure S5), e.g., ACO3, CDSP32, GSH1, etc., were up-regulated while others were down-regulated (GR1, GASA14, and GSTF11) in both the initial (apomeiosis) and later (parthenogenesis) stages of ovule development. Nonetheless, apomicts present more up-regulated DEGs regarding oxidative stress than sexual plants, supporting the involvement of redox reactions in this reproductive mode. Alteration of homeostasis-based processes of stress perception and attenuation in sexual species of several genera would induce apomeiotic spores and gametophyte formation [23]. Apomeiosis occurs when the redox balance is more toward H 2 O 2 catabolism, and the transition from meiosis to apomeiosis can be changed by a disturbance in this homeostasis [23].
Among DEGs between sexual and Limonium apomicts, most were up-regulated in the latter stages of ovule development (parthenogenesis), such as ATRX genes coding for chromatin remodelling proteins as well as multiple histone methylation genes concerning epigenetic developmental mechanisms in plants [76] (Tables 5 and 6). These DEGs were also previously found to be upregulated, for instance, in parthenogenetic eggs of Cenchrus ciliaris [77]. Moreover, other DEGs implicated in small RNA biogenesis and DNAmethylation pathways, such as the AGO9 and AGO4 homolog genes in Arabidopsis mutants found to be associated with phenotypes reminiscent of apospory or diplospory [78], were also detected in our study. In Boechera apomicts, AGO9 was found at low levels in the megaspore mother cell itself, becoming an apomictic initial cell [79]. However, in our study, AGO4 and AGO9 seem to have more specific roles in ovules at later stages of development (parthenogenesis), likely being involved in eggs assuming a parthenogenesis fate. MAPK signaling and aminoacyl-tRNA biosynthesis pathways that perform roles in translational regulation, RNA splicing, and tRNA proofreading [80] also showed transcriptional changes at this stage in apomictic Limonium plants (Tables 5 and 6).
DEGs were also remarkably enriched in genes implicated in hormonal signaling, such as the ethylene signaling pathway, in which the apomictic gametophytes overexpressed 26 ethylene-responsive transcription factors (Tables 5 and 6). For example, in Cenchrus ciliaris, EIN2 (ethylene insensitive 2) together with 14 ethylene responsive transcription factors were up-regulated in parthenogenetic eggs [77], although in our study both genes showed contrasting expression patterns in the same developmental stage (parthenogenesis). Moreover, in the parthenogenetic ovules, overexpressed genes were related to tryptophan metabolism, such as TAA1 (l-tryptophan pyruvate aminotransferase), which converts tryptophan to indole-pyruvic acid, a direct biosynthetic precursor of the auxin in Arabidopsis (IAA [81]). Crosstalk between ethylene signaling and auxin pathways is involved in the regulation of developmental processes [82].

Feminization of Apomicts Is Related to Down-Regulation of Floral Genes Specifying Stamens
Besides auxin, other hormones like gibberellic acid (GA) contribute to flower development, the development of male and female gametophytes, and seed germination [83.84]. The GASA genes as well as the GA biosynthesis genes in Arabidopsis, implicated in controlling floral induction, seed maturation, and germination [83,84], were differentially expressed between apomictic and sexual plants in our study ( Table 4). One of the targets of GA signaling are the floral homeotic genes encoding MADS-box transcription factors involved in floral development in accordance with the ABCDE model [85]. In our study, among the top DEGs between sexual and apomictic plants and between the different ovule stages, MADS-box transcription factors were identified, including floral homeotic genes with a MADS-box domain. In A. thaliana, the MADS-box from A-class genes (AP1; AP2) specifies the formation of sepals; the combination of A-and B-class genes (AP3; PI) determines petal' development; the B-, C-(AGL), and E-class (SEP) genes specify stamens; and the C-and E-class genes specify carpels. Only the expression of genes from class C specifies carpel formation. Class E genes (SEP3) are associated with the formation of all flower whorls. The gene classes A and C are expressed antagonistically; the A gene class is expressed in sepals and petals, and the C gene class is expressed in stamens and carpels [86,87].
GA promotes reproductive development by upregulating expression of the floral meristem identity gene LEAFY (LFY), which in turn upregulates expression of AP3 and AGL that, in conjunction with PI and SEP3, regulate floral organ identity [87]. In our study, we found changes in the expression of MADS-box genes in the different stages of ovule development between sexual and apomictic plants, particularly PI, SEP1, and AP3 genes, which were downregulated in apomicts. The PI/AP3 genes have a role in sexual dimorphism and have been identified as masculinizing factors in spinach [88]. In dioecious plants such as Populus, constitutive overexpression of PI/AP3 produces male flowers, but in female flowers the presence of a feminizing factor F downregulates PI/AP3, diverting development to a female developmental pathway, inhibiting stamens, and allowing carpels to form [89]. Therefore, it could be hypothesized that male sterile Limonium plants with homeotic changes in floral organs lack the gene function of the corresponding class B genes.

Male Sterility Appears to Be Linked with Downregulation of Genes Connected to Pollen Wall Formation and Assembly and Pollen Tube Growth
Various genes are involved in pollen wall development and assembly, which is a specialized extracellular cell wall matrix that encases the male gametophytes [90]. A specific cell wall polymer also known as β-glucan is synthesized by callose synthases in pollen mother cells and microspore tetrads that acts as a template for primexine, thus providing a structural basis for exine formation [90]. Callose synthase 5 (CSL5) is a key isoform of callose synthases responsible for the formation of the callose wall, which is essential for the accumulation of callose in the tube wall and the callose plug in growing pollen tubes [91]. In CSL5 mutants, the viability of pollen grains is greatly reduced in Arabidopsis [91,92] and rice [93]. In our study in Limonium, CSL5 is downregulated in the initial (apomeiosis) and later phases of development (parthenogenesis) in the apomicts. One of the characteristic features of plants in the L. binervosum group is the widespread existence of male sterility, i.e., a lack of pollen [12,13,26,36]. Electron microscopy studies showed that L. multiflorum plants had many flowers with empty anthers and sometimes flowers with no pollen at all; the few microspores that formed showed collapsed morphology and lacked the typical exine patterns [12]. In L. multiflorum apomicts, after anther dehiscence releases pollen, the plants never undergo their first mitosis, only attaining the "ring vacuolate" stage, and the male germ unit is not produced [12]. Interestingly, in this study, callose synthase isoforms such as CSL9 or CSL11 that were upregulated in the apomicts have a role in Arabidopsis pollen mitosis by disrupting pollen mitosis and producing pollen with only one or two nuclei, the generative cell being degenerated, undifferentiated, or mislocalized [94,95], as found here in Limonium.
Moreover, in our study, the gene PEX4 that codes for extracellular glycoproteins that belong to the hydroxyproline-rich glycoprotein family and have a role in pollen germination and pollen tube growth in A. thaliana [96,97]. was downregulated in Limonium apomicts (apomeiosis, M1a1). The PEX4 mutants have an excessive deposition of callose [96,97], leading to abnormal pollen tubes that develop bulges and burst [97]. While in L. multiflorum apomicts, pollen tubes are never observed since unicellular pollen never undergoes its first mitosis. In L. ovalifolium sexual plants, pollen grains follow a first asymmetric mitotic division, producing a generative cell within the vegetative pollen grain cell in the binucleate pollen stage [12]. Therefore, our results support a role for PEX4 in pollen tube growth. Perhaps this gene, along with other unknown genes, might create a terminal combination that does not allow the development of pollen tubes.

Pollen-Stigma Interactions
Pollen-pistil interactions can be viewed as a major prezygotic pollination reproductive barrier and are active systems of pollen rejection [98]. Some of these systems act at the level of the stigma, with a genetic control independent from embryo sac development involving S-alleles [99]. In the Brassicaceae family, the SI is controlled sporophytically by a single S locus that incorporates stigma-expressed and anther-expressed genes composed of multiple alleles or variants [100]. In various gametophytic apomicts, non-functional pollen can cause a weakening of SI and a breakdown of the sporophytic SI system (mentor effects [9]).
Limonium species show a polymorphic sexual system associated with flower polymorphisms and a sporophytic self-incompatibility that prevents self-and intramorph mating [29.30]. Limonium gametophytic apomicts that form diplosporic (apomictic) embryo sacs of Rudbeckia type like in L. multiflorum [13] show abnormal and non-functional pollen due to a sporophytic defect [12]. In our study, in the GO term "recognition of pollen" (GO:0048544; Table 6). All DEGs were lectin receptor kinases (LECRKs) [101], except for AT3G49500. These kinases belong to the class of G-LECRKs [101], particularly the S-locus Receptor Kinase (SRK), known for its role in self-incompatibility [101] and potentially of high interest in our studied species. In our study, DEGs from this LECRKs complex were detected under the GO term "recognition of pollen", namely AT1G61380, AT1G65800, AT5G24080, AT4G27300, AT4G21380, and AT4G27290. Four DEGs of the SRK were overexpressed, namely in the initial stages of ovule development (apomeiosis), such as AT1G61380 (M1a1 and M1a2) and AT1G65800 (M1a1), as well as in later ovule stages (M2a4 and M2o4), AT4G27300 (M2o4), and AT4G27290 (M2o4). These findings indicate that the genes implicated in pollen recognition in Limonium were already expressed at earlier ovule stages. This implies that the fate of the gametophytic apomict male spores is decided by the maternal genes in the initial ovule stages. Nonetheless, a genetic linkage between SI and apomixis cannot be easily assumed, since the breakdown of SI mostly affects the pollen and the stigma, while apomixis affects the development of the embryo sac.

Up-Regulation of Specific Genes Related with Embryo Formation in the Apomicts
In Arabidopsis, flowering can be promoted by repressing the transcription of the central flowering repressor and vernalization regulatory gene FLC (FLOWERING LOCUS C), which belongs to the MADS-box class of transcription factors [102]. The FLC seems to regulate several transcription factors involved in important biological processes such as reproductive and embryonic development [103]. FLC impedes the floral transition by inhibiting the expression of the floral primordium identity genes such as FT (FLOWERING LOCUS T), SOC1/AGL20 (SUPPRESSOR OF CONSTANS OVEREXPRESSION 1/AG20), LFY, AP1, and floral organ identity AG and AP3 genes [104]. Another gene involved in the genetic control of flowering time in Arabidopsis is FPF1 (FLOWERING PROMOTING FACTOR), which modulates the acquisition of competence to flower in the apical meristem and is expressed earlier than AP1 [105]. In our study, the FPF1 and SOC1 showed higher levels of regulation in apomicts than in sexual plants, whereas AP1 showed reduced levels in the first ones. These results indicate differences in the regulation of major genes controlling the transition from a vegetative to a reproductive mode in the apomict's apical meristem. Interestingly, in our study, both FLC, AGL6, and AGL15, as well as other MADSbox transcription factors, were specifically upregulated in the later stages of development. In A. thaliana AGL6 functions in the early stages of the flowering signal transduction pathway by inhibiting the transcription of FLC genes [106]. In Brachiaria brizantha, AGL6 is differentially expressed during embryo sac formation in apomictic and sexual plants [107]. In our study, AGL15, which plays an essential role during early zygotic embryogenesis in Arabidopsis [108], in Brassica napus, and in soyabean somatic embryos [109], was specifically upregulated in the later phase of ovule development in Limonium apomicts. The AGL15 is a component of the SERK protein complex [110] that is part of a molecular network linked to zygotic and somatic embryogenesis. However, in this study, we were unable to find any differential expression of SERK-like genes.
Some genes related to the bioactive gibberellins' deactivation reaction from the GA2OX family were differentially expressed. For example, GA2OX6 was up-regulated in both the initial and later stages of ovule development. While in A. thaliana, GA2OX6 expression is activated by AGL15 during embryogenesis [111], in our study, AGL15 was down-regulated in M2a4 and up-regulated in M2o4, suggesting differences among species. Moreover, GA2OX6 is found to be expressed in sepals, stigmas, and immature anthers. Regarding seed development, it was only expressed in the antipodal cells before the 8-cell stage, suggesting that this gene is a negative regulator of seed germination [111]. Remarkably, from the same family, GA2OX8, which is exclusively expressed in stomatal cells in A. thaliana but is not expressed or has distinct expression patterns in flower tissues or seed development [111], was down-regulated in Limonium apomicts in later stages of ovule development (parthenogenesis). This finding supports the idea that this gene can have other or different roles in Limonium.

Conclusions
This study sheds light on genes involved in Limonium sexual and apomictic reproduction. The findings substantiate the deregulation of gene expression in the regular sexual pathway. While several HKGs are found to be differentially expressed between sexual and asexual plants, other genes are found to be stably expressed. These can be potentially used as reference genes to be validated for specific tissue types (vegetative and reproductive) or reproductive modes (sexual and apomictic) in future quantitative gene expression studies. Our findings reveal that the latter stage of ovule development (parthenogenesis) was the most contrasting phase in terms of differential gene expression between asexual and sexual plants. Among them, the MADS-box domain TFs are central players in many developmental processes, including control of flowering time, homeotic regulation of floral organogenesis, fruit development, and seed pigmentation.
Since L. multiflorum male sterile plants form parthenogenetic ovule sacs, it could be interesting to analyze candidate genes such as PEX4 in pollen tube development and the function of AGL genes (e.g., AGL6) specifically modulated in the latter stages of development (parthenogenesis). Nonetheless, given the high number of 71% unannotated genes in Limonium, other studies are required to clarify the regulatory roles of these genes.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/genes14040901/s1. Figure S1: Principal component (PC) analysis of all rlog-transformed gene expression data; Figure S2: Weighted Venn diagrams of specific and overlapping differentially expressed genes (DEGs) found in the ovules of sexual L. auriculifolium and L. ovalifolium; Figure S3: Weighted Venn diagrams of specific and overlapping differentially expressed genes (DEGs) found in different stages of apomictic L. multiflorum, sexual L. auriculifolium, and L. ovalifolium ovules; Figure S4: Number of differentially expressed genes (DEGs) common to different comparisons with opposite regulations; Figure S5: Distribution of differentially expressed genes (DEGs) related to oxidative stress; Table S1: Quantification of basic quality and completeness metrics of Limonium de novo transcriptome assembly. Table S2: List of annotated differentially expressed genes (DEGs) shared between two comparisons with opposite regulation in ovules from Limonium plants; Table S3: Uniquely annotated top down-and up-regulated differentially expressed genes (DEGs) in different ovule stages in L. multiflorum relative to L. auriculifolium; Table S4: Uniquely annotated top down-and up-regulated differentially expressed genes (DEGs) in different ovule stages in L. multiflorum relative to L. ovalifolium; Table S5: Uniquely annotated top down-and up-regulate differentially expressed genes (DEGs) in different ovule stages in L. multiflorum relative to L. dodartii; Table S6: Uniquely annotated top down-and up-regulate differentially expressed genes (DEGs) in ovules from apomictic L. multiflorum; Table S7: Uniquely annotated top down-and up-regulated differentially expressed genes (DEGs) in ovules from sexual Limonium auriculifolium; Table S8: Uniquely annotated top down-and up-regulate differentially expressed genes (DEGs) in ovules from sexual L. ovalifolium; Table S9: Uniquely annotated top down-and up-regulate differentially expressed genes (DEGs) in ovules from sexual L ovalifolium (O) relative to the control L. auriculifolium; Table S10: Overrepresentation analysis (ORA) performed by gProfiler of transcription factors (TFs) and differentially expressed genes (DEGs).