The Embryo Characteristics and Genes Expression Analysis Involved in Apomixis of Wild Germplasm Materials of Kentucky Bluegrass in Gansu Province of China

Background: Apomixis is mainly used to maintain the heterosis, stability, and consistency of crops. Its main advantage is to reduce the costs of seed production and shorten the breeding process. In the eld of hybrid breeding, apomixis has been referred to as the "asexual revolution" and has realized a new green revolution.Apomixis is mainly used to maintain the heterosis, stability, and consistency of crops. Its main advantage is to reduce the costs of seed production and shorten the breeding process. In the eld of hybrid breeding, apomixis has been referred to as the "asexual revolution" and has realized a new green revolution. Kentucky bluegrass (Poa pratensis) is a natural apomictic species that mainly exhibits facultative apospory. Its main feature is that the somatic nucellar cells bypass meiosis and double fertilization to form unreduced embryos, and the development of the endosperm requires pseudogamy. Although apomixis is of great signicance in breeding, the genetic control of apomixis remains unclear. Therefore, we report the results of a global gene expression analysis of wild germplasm materials of Kentucky bluegrass spikelets in Gansu province of China, exhibiting signicant differences in apomictic rates to identify the genes, biological processes, and molecular functions related to apomixis. Result: At four reproductive periods, there were 5400 differentially expressed genes (DEGs) between the two genotypes, including 2476 downregulated and 2823 upregulated genes. Further analysis of the gene functions, pathways, expression patterns, networks, and transcription factors (TFs) showed that the occurrence of apomixis in Kentucky bluegrass was related to changes in the time- and space-related expression of genes associated with sexual reproduction, which led to disordered sexual reproduction and thus the production of offspring by apomixis. Conclusion: At the transcriptional level, the genesis and development of apomixis was regulated by TFs. It also involved the coexpression of many genes associated with disordered meiosis, hormone signal transduction, embryonic development, stress response pathways, and epigenetics. We also veried these 16 DEGs by real-time quantitative PCR. The expression results were basically consistent with the transcriptome expression prole, indicating that the transcriptome data were reliable. The results of this study may provide a theoretical basis for revealing the mechanism of occurrence and development of apomixis in Kentucky bluegrass. genes indicates a large-scale change in the regulation of apomixis ovules[44, 45]. Sharbel et al.[45] suggested that TFs are important motifs for the differential expression between sexual and apomictic ovules, and their mediation of the inhibition of gene expression is a characteristic of the development of sexual ovules to apomictic ovules. These studies suggest that TFs may regulate the transformation of plants from sexual reproduction to apomixis[46], and thus we further analyzed the differential expression of TFs in the two genotypes of Kentucky bluegrass. The results showed that a total of 84 TF families were differentially expressed in the two genotypes of Kentucky bluegrass, of which 52 were common across all four stages. nucellar cells embryogenesis ability. Our results showed that the stress response plays an important role in the development of apomixis

aestivum), rice (Oryza sativa), and maize (Zea mays), its apomictic trait cannot easily be transformed into cereal crops by hybridization, and there is a lack of genomic information available for Kentucky bluegrass. In this case, it is crucial to analyze its gene expression patterns during embryonic development.
Therefore, in this study, we selected wild Kentucky bluegrass germplasm resources in Gansu exhibiting signi cant differences in apomictic rates and used high-throughput sequencing to obtain spikelet reference transcriptomes of two genotypes of Kentucky bluegrass, screen and identify relevant DEGs, and clarify the biological processes that these genes may participate in. Our ndings should help elucidate the molecular basis of apomixis in Kentucky bluegrass and provide basic data for the cultivation of excellent new varieties of Kentucky bluegrass through apomixis.

Experiment material
We collected wild germplasm materials of Kentucky bluegrass from various areas in Gansu, China, in the early stage, and the apomictic rate of these wild germplasm materials was calculated using cytoembryology, as described in our previous study [80]. In this study, wild germplasm materials of Kentucky bluegrass collected from Longnan (LN) and Gannan (GN), which exhibited signi cant differences in apomictic rates, were selected as research materials.
Information on the collection sites is shown in Table 4. In late April 2019, GN and LN were planted in the lawn training base of Gansu Agricultural University (longitude 103°34′ east, latitude 36°5' north). The climatic conditions and soil information of the lawn training base are shown in Table 5. The spikelets of Kentucky bluegrass in the premeiosis (PreM), meiosis (M), postmeiosis (PostM), and anthesis (A) stages were collected according to the panicle type reported by Albertini et al. [25]. Some of the spikelets were immediately xed in FAA (formalin: glacial acetic acid: 70% ethanol = 1:1:18) and then used to make para n sections. Some of the spikelets were rapidly frozen in liquid nitrogen and stored at -80°C for endogenous hormone determination and RNA extraction. Three biological replicates were set at each time point, and each sample was taken from a mixed sample of at least ve plants.

Cytoembryological observation of the ovaries of Kentucky bluegrass
After the orets of the spikelets stored in 70% ethanol were separated, the samples were subjected to an ethanol dehydration series. Xylene was applied as a transparency agent, the samples were dipped in para n wax, and a split para n embedding machine (EG1150H+EG1150C) was then used for embedding. A Leica RM2265 electric rotary automatic microtome was used to produce sections of 10 μm thickness. Ehrlich (Ehrlich's) hematoxylin (1 min)-eosin (2 min) double staining was conducted in a Leica RM2265 automatic dyeing machine, following which a Revolve RVL-100-G inverted microscope was used to observe and photograph the samples.
RNA extraction, library construction, and Illumina sequencing Total RNA was extracted from all samples using an RNasy Plant Mini Kit according to the manufacturer's instructions (Qiagen, Hilden, Germany). After accurate detection of RNA by 1% agarose gel electrophoresis as well as a NanoPhotometer spectrophotometer (IMPLEN, Westlake Village, CA, USA) and Agilent 2100 biological analyzer (Agilent Technologies, Palo Alto, CA, USA), the cDNA library was constructed using Illumina's NEBNext ®UltraTM RNA Library Prep Kit (San Diego, CA). First, polyA mRNA was screened and enriched by Oligo (DT) magnetic beads. The mRNA was then randomly interrupted in NEB Fragmentation Buffer and used as a template, and random oligonucleotides were used as primers to synthesize the rst-strand cDNA in a M-MuLV reverse transcriptase system. Third, the RNA chain was degraded by RNaseH, and the second chain of cDNA was synthesized from dNTPs in the system of DNA polymerase I. Finally, cDNA of 250-300 bp was obtained and subjected to end-repair and polyA tail addition and screened using AMPure XP beads. PCR ampli cation and puri cation of the PCR products were carried out, and the nal library was obtained.
The constructed library was sequenced on an Illumina HiSeq 4000 platform (2×150 bp reading length) after passing a quality control test. Raw reads were obtained by sequencing and were ltered to obtain clean reads. The transcripts of the clean reads were assembled by Trinity [81] software, and the longest among similar classes was selected as the unique gene by CD-HIT [82], and a collection of all gene sequences was obtained as a reference sequence for analyzing gene expression. Finally, RSEM software [83] was used for reference sequence mapping of the clean reads.
Functional annotation of unigenes BLASTX was used to search unigenes in the Nr (NCBI non-redundant protein sequences), Nt (NCBI nucleotide sequences), Swiss-Prot, and KOG (euKaryotic Ortholog Groups) databases, and the Pfam (Protein family) database was used to predict protein domains.
The expression level of each transcript was determined by FPKM (the number of fragments per thousand exons per million mapped fragments). The DESeq R Bioconductor package [84] was used to determine the differential expression of expressed genes based on the negative binomial distribution model. Benjamini and Hochberg [85] was used to control the false discovery rate (FDR). A |log2 (FoldChange)| ≥1 and FDR <0.05 were used as the criteria for selecting DEGs.
When the Bonferroni-corrected P-value was <0.05, a DEG was considered to be signi cantly enriched. Using the method of GOseq [86], unigenes were annotated according to biological process, cellular composition, and molecular function in the GO (Gene Ontology) database, so as to obtain the biological function of signi cantly enriched DEGs. KEGG (Kyoto Encyclopedia of Genes and Genomes) and related software KAAS (KEGG Automatic Annotation Server, Quantitative real-time (qRT) PCR veri cation Sixteen genes were randomly selected and speci c PCR primers were designed by using the NCBI online website Primer-BLAST. The total RNA was extracted using an RNeasy Plant Mini Kit according to the manufacturer's instructions (Qiagen, Hilden, Germany). The total RNA was reverse-transcribed into cDNA using a Prime Script RT reagent Kit with a gDNA Eraser (Perfect Real Time) (TaKaRa, Japan) and diluted 20 times as a template. The qRT-PCR was performed with a StepOnePlus™ Real-Time PCR System (ABI) using SYBR premix Ex Taq (Takara, Japan). The reaction volume was 20 μL and included 2 × SuperReal PreMix Plus 10 μL, ddH 2 O 3 μL, cDNA template 5 μL, and forward and reverse primer 1 μL. The reaction conditions included pre-denaturation at 95°C for 15 min and 40 cycles of PCR ampli cation, including denaturation at 95°C for 10 s and annealing at 58°C for 30 s. Using Actin as the internal reference gene, the relative expression level of the selected gene was calculated by 2 −∆∆Ct [88].

Results
Histological study on the ovaries of Kentucky bluegrass The sexual embryo sac of Kentucky bluegrass is a typical polygonum-type embryo sac, and the normal sexual reproductive embryo sac develops from the MMC through meiosis (Fig. 1A-G). The MMC develops from a cell located at the top of the central axis cell mass under the nucellus epidermis. Compared with other nucellar cells, this cell is large in size and has a small and dense cytoplasm, and the nucleus accounts for the main part of the whole cell, which is called an archesporial cell (Fig. 1A). The archesporial cell grows and develops into MMC (Fig. 1B). Then the MMC undergoes meiosis to form a dyad and tetrad ( Fig. 1C and D). With growth and development, the tetrad forms a megasporangium (Fig. 1E). The megasporangium undergoes mitosis gradually to form a two-nucleus, four-nucleus embryo sac and a "7-cell-8-nucleus" structure ( Fig. 1F and G). At this point, the development of the female sexual gametophyte of Kentucky bluegrass is complete.
Developmental variety in the AICs was observed ( Fig. 1H-P). First, the multiple nucellar cells in the ovary were specialized, and the cell volume was signi cantly larger than that of the surrounding nucellar cells (Fig. 1H). During the stage of the archesporial cell, the specialization of the nucellar cells also occurred ( Fig. 1I). Subsequently, the specialized nucellar cell developed and expanded to form an AIC, which was large in size and had an obvious nucleus. The AIC could coexist with MMC, and they could be arranged horizontally (Fig. 1J) or longitudinally (Fig. 1K). In addition, multiple specialized nucellar cells developed together to form multiple AICs. At this time, the MMC may disintegrate and form nutrients that are absorbed by the AICs (Fig. 1L). In some ovaries, the MMC developed normally and coexisted with two AICs (Fig. 1M). The AIC could form and coexist after the MMC formed a dyad (Fig. 1N). The AIC formed a longitudinally arranged two-nucleus embryo sac that was clearly distinguished from the development of sexual reproduction after mitosis (Fig. 1O). In addition, the existence of apomixis led to the formation of a three-nucleus embryo sac (Fig. 1P).

Changes in endogenous hormone contents in the spikelets of Kentucky bluegrass
The change trend in ZT content of the spikelets of the two genotypes was basically the same, increasing at rst and then decreasing, and reaching the maximum value at meiosis. One difference was that the increase and decrease in GN materials were higher than that in LN ( Fig. 2A). The variation trend of GA 3 content in both genotypes rst increased and then decreased, being highest at meiosis in LN and highest at the postmeiosis stage in GN (Fig. 2B). The content of IAA in LN exhibited no difference between the rst two periods and was signi cantly lower than that in the latter two periods. In GN, the content of IAA decreased rst, then increased, and nally decreased. Comparing the two genotypes, the IAA content of LN was signi cantly lower than that of GN in the rst two periods and was higher in the latter two periods (Fig. 2C). The ABA content increased gradually with the development of the spikelets in both genotypes, and there was no difference between the two genotypes in the rst three periods. Only at anthesis was the ABA content in LN signi cantly higher than that in GN (Fig. 2D).

RNA sequencing analysis
Each sample had more than 21.03 million raw reads, and a total of 152.35 Gb clean reads were obtained after ltering. With the exception of LM, the error rate of the other samples was 0.02%. The Q20 and Q30 percentages were greater than 98.12% and 94.46%, respectively, and the GC content percentage ranged from 52.17% to 54.06% (Table 1). Furthermore, the clean reads of each sample were mapped to obtain the total mapped reads after ltering. The percentage of total mapped reads in all samples to total reads was more than 63.05% (Table S1). In addition, we found there were 833 (57.85%) complete sequences, 366 (25.42%) fragment sequences, and 241 (16.74%) homologous genes (Fig. S1).

Transcriptome de novo assembly and annotation
A total of 892774 transcripts and 288678 unigenes were obtained. The length and quantity statistics of the transcripts and unigenes are shown in Fig. 3A.
According to the statistics, there were 300685 and 114763 transcripts and unigene sequences with a length of 300 to 500 bp, respectively. The number of transcripts between 501 and 1000 bp was the largest, accounting for 40.90% of the total. The number of transcripts and unigenes larger than 2000 bp was the lowest, accounting for only 4.70% and 4.92% of the total, respectively (Fig. 3B, Table S2, Fig. S2). Furthermore, functional annotation information of the unigenes of the Kentucky bluegrass spikelet transcriptome was obtained (Fig. 3C, Table S3). The number of annotations in the Nr database was the largest, accounting for 47.51%. The results showed that the spikelet gene sequence of Kentucky bluegrass had the highest homology with Aegilops tauschii, followed by Brachypodium distachyon and Hordeum vulgare. The homology rates of Triticum aestivum and Triticum urartu accounted for 7.3% and 7.2%, respectively (Fig. 3D). In addition, analyzing the distribution of its e-value, it was found that only 1.2% of the unigenes annotated to the Nr database were compared to an e-value equal to 0 (Fig. S3A). The distribution of sequence similarity showed that 33.8% of the Nr database-annotated unigenes reached 80%-95% similarity, and only 15.2% of the sequence similarity was less than 60%, which indicated that the sequence similarity of the species was relatively high (Fig. S3B).

Identi cation and analysis of DEGs
In LN, the number of DEGs was the highest in LA vs. LPreM, which had 11849 upregulated and 19019 downregulated genes. The number of DEGs was lowest in LPostM vs. LM, which included 176 and 147 upregulated and downregulated genes, respectively (Fig. 4, Table 2). The change in GN was the same as that of LN, in which the number of DEGs was the largest in GA vs. GpreM and the lowest in GPostM vs. GM (Table 2, Table S4). To analyze the speci c expression of genes of the four stages and visually illustrate the common and unique DEGs, we constructed a Venn diagram of the DEGs (Fig. 4). We further analyzed the DEGs between the two genotypes. It was found that 23234 DEGs were identi ed at PreM, of which 13583 were upregulated and 9651 were downregulated. A total of 37773 DEGs were identi ed in LM vs. GM, including 19507 upregulated and 18266 downregulated genes, respectively. The number of upregulated and downregulated genes identi ed was 22722 and 23446, respectively, at PostM. A total of 15691 and 9766 upregulated and downregulated DEGs were identi ed in LA vs. GA, respectively (Fig. 4C, Table S4). There were 5400 DEGs expressed at all four stages. The number in PreM was the lowest, only accounting for 8.76% of the total DEGs, while the percentage of DEGs in PostM was the highest (23.05%).

Functional annotation of DEGs
We selected GO terms that were signi cantly enriched in at least one stage and that may be related to apomixis to construct a heatmap showing expression characteristics (Fig. 5). We found that the DEGs were generally related to stress response (regulation of response to stress, response to oxidative stress, and oxidoreductase activity), hormone signal (regulation of hormone levels, hormone metabolic process, cellular hormone metabolic process), meiotic (reciprocal meiotic recombination), female gamete generation (female gamete generation), and cell wall (cell wall modi cation) (Fig. 5, Table S5). Among them, only oxidoreductase activity was signi cantly enriched at all four stages, cell wall was signi cantly enriched at three stages, and the other 10 terms and 6 terms were only signi cantly enriched at two and one stages, respectively. To further explore the biological pathways of the DEGs, we performed KEGG analysis (Table S6). The results indicated that the pathways with the highest scores were photosynthesis−antenna proteins, thiamine metabolism, avone and avonol biosynthesis, and cutin, suberine, and wax biosynthesis. According to existing research reports, we listed the pathways that may be related to apomixis in this study (Table 3). Although the scores of these pathways are not the highest, they were still found.

Analysis of Transcription Factors (TFs)
TFs are key to plant growth, development, and the regulation of gene expression in response to abiotic and biotic stress. The top 20 TF families are shown at each stage in Fig. 6 and Table S7. Compared with GN, the TF families with the most DEGs in ProM were NAC, FAR1, and C2H2 (Fig. 6A). The top three TF families in the meiotic stage and the PostM stage were the same and included NAC, FAR1, and TARF ( Fig. 6B and C). This was quite different from the previous three periods to anthesis, in which the TF families with the most DEGs were C2H2, FAR1, and bHLH (Fig. 6D). In addition, there was a signi cant difference in the number of DEGs corresponding to the upregulated and downregulated adjustments at each stage. There was little difference in the number of differentially expressed TF families in the four stages, ranging from 64 to 75. The stage PostM had the highest numbers of DEGs, and the PreM stage had the lowest numbers of DEGs (Table S7).
Analysis of gene expression patterns, DEG clustering, and functional enrichment According to the clustering heatmap analysis of the DEGs, the gene expression pro les underwent a series of changes between the different genotypes and/or at different developmental periods (Fig. 7A). According to the gene expression pro le, the DEGs of two different genotypes at different developmental stages were clustered into 10 pro les ( Fig. 7B and C). It was found that with time, the two genotypes exhibited differences in gene expression in response to spikelet growth and development. Pro le 0, 7, 8, and 9 were signi cantly overexpressed in LN (P<0.05, Fig. 7B, Table S8). The difference was that Pro le 0, 8, and 9 were obviously overexpressed in GN (P<0.05), and Pro le 7 was not overexpressed (Fig. 7C, Table S9). To determine the functional signi cance of the transcriptional changes of each genotype, GO and KEGG function enrichment analysis were performed on genes belonging to the overexpression pro le ( Fig.  7D-G). In LN, with the exception that response to hormone, gene silencing by RNA, and gene silencing were not enriched in Pro le 0, the remaining terms were all signi cantly enriched. There were signi cantly less terms enriched in Pro le 7 and 9 than in Pro le 0. Only protein binding in Profile 8 was signi cantly enriched (Fig. 7D). The KEGG enrichment results showed that RNA transport, RNA degradation, and phenylpropanoid biosynthesis were only signi cantly enriched in Pro le 0. Zeatin synthesis, ribosome, and oxidative phosphorylation were only signi cantly enriched in Pro le 9. No pathway was enriched in Pro le 0 and Pro le 7 (Fig. 7E). In GN, RNA binding was signi cantly enriched in Pro le 0 and Pro le 8. Oxidoreductase activity, oxidation-reduction process, cell wall modi cation, and cell wall were signi cantly enriched in Pro le 0 and Pro le 9, and other terms were signi cantly enriched only in a single pro le (Fig.   7F). Pathway enrichment analysis showed that RNA transport was signi cantly enriched in Pro le 0 and Pro le 8. Phenylpropanoid biosynthesis and avone and avonol biosynthesis were signi cantly enriched in Pro le 0. RNA degradation and DNA replication were signi cantly enriched in Pro le 8. The remainder of the pathways were only signi cantly enriched in Pro le 9 (Fig. 7G).

Gene coexpression network analysis
Weighted gene correlation network analysis (WGCNA) has a scale-free topological structure and is a systems biology method that is used to describe gene association patterns among different samples. It can be used to identify highly synergistic gene sets and candidate biomarker genes based on the interconnectivity of the gene sets and the association between the gene set and the phenotype. According to the pairwise correlation and gene expression trend of all samples, the unigene expression data from all 24 samples (two genotypes, three biological repeats, and four time points) were constructed into a coexpression network. The hierarchical clustering tree displays the coexpression module identi ed by WGCNA, and each leaf on the tree represents a gene.
The main branch consisted of 18 modules (from a low serial number to a high serial number, M1 to M28), which are marked with different colors (Fig. 8A, Table S10). The 18 modules contained a total of 19981 unigenes, and the module with the least number of genes was midnightblue, which contained only 81 genes. The largest number of genes was in the turquoise module, which contained 4966 genes. The average number of genes per module was 1110, and the genetic correlation between the modules was low (Fig. 8B, Table S10), which indicated that the analysis results of the module were reliable. In addition, to clarify the correlation between different modules, we also further conducted clustering and heatmap analysis. The two analyses exhibited the same result, both showing a high degree of correlation between the MEblue and MEred, MElightcyan and MEmidnightblue, MEgreen and MEyellow, MEcyan and MEgrey60, MEblack and MEpink, and MEgreenyellow and MEpurple modules (Fig. 8C, Table S11, S12). Finally, the correlation between the module and the sample was analyzed, and the result indicated that MEcyan and MEgrey60 had a high correlation with sample LA_2. In addition, MEsalmon, MEpink, and MEgreenyellow had a higher correlation with LPostM_1, GM_3, and GPostM_2, respectively (Fig. 8D, Table S13 and S14).

qRT-PCR veri cation
Sixteen genes possibly related to apomixis in Kentucky bluegrass were selected for qRT-PCR to verify the authenticity and reliability of the transcriptome results (Fig. 9, Table S15). The results indicated that the genes basically showed the same expression patterns and the same trend of change between the qRT-PCR and RNA-seq results. However, the expression of CAT2 (catalase isozyme 2), which was downregulated at PostM in our results compared to GN, was not consistent with the RNA-seq results. Excluding the postmeiosis stage, the change trends in the premeiosis, meiosis, and anthesis stages were consistent with the RNA-seq results. The expression trend of the other 15 DEGs at four time points was the same as that of the RNA-seq results. Therefore, we conclude that the qRT-PCR results validated the transcriptional map data obtained from our RNA-seq analysis, con rming that the transcriptome sequencing results were accurate and reliable.

Discussion
Embryological study on apomixis in Kentucky bluegrass Apomixis can x any genotype across generations, and Kentucky bluegrass is a natural apomictic species closely related to major crops. The existing research on apomixis in Kentucky bluegrass has focused on various aspects, particularly cell embryology. In the late 1970s, Yong et al. [32] used the techniques of pistil cleaning and thick sectioning to observe the apomixis of Kentucky bluegrass. Subsequently, Marshall et al. [33] noted that the apomixis of Kentucky bluegrass was mainly aposporic. In addition, Kentucky bluegrass possesses a small number of parthenogenetic embryos, cocellular embryos, and antipodal embryos [11]. The occurrence of apospory in Kentucky bluegrass is caused by the expansion of one or more somatic cells in the ovules, which restores reproductive ability, and the unreduced embryo is then produced by several mitotic events [34,35]. The number of AICs is mostly one or two, though three are occasionally observed [36]. In this process, sexual embryos and apomictic embryos can coexist in the same ovule [37] and ultimately produce seeds with two or three embryos [11]. These previous reports support our research results, which demonstrated that during the embryonic development of Kentucky bluegrass, one or more nucellar cells located near the MMC under the nucellus epidermis became enlarged and exhibited an obvious nucleus, and their meristematic ability was restored. They then further developed to form AICs. The AICs developed into unreduced and unfertilized apomictic embryos, which could coexist with the sexual embryos.

Transcriptome analysis of apomixis in Kentucky bluegrass
Plant apomixis was rst discovered in Alchornea ilicifolia in 1841 [48]. Since this discovery, there has been much research on the classi cation, genetic evolution, embryonic development, morphogenesis, physiology, biochemistry, and breeding of apomixis [38,39]. Researchers have long tried to collect natural apomixis resources and analyze their biological characteristics and genetic mechanisms, as well as introduce key genes into crops to make use of apomixis xed heterosis in crops. However, due to the complexity of species, occurrence processes, and formation mechanisms, as well as the diversity in apomixis types, it is extremely di cult and time-consuming to mine and clone related genes using traditional breeding methods. Thus, a reverse genetic strategy based on transcriptomics, which is used to screen and identify candidate genes and then verify them by experimental mapping and functional analysis, has become the core of apomixis studies [23]. Furthermore, transcriptome methods are applicable to all types of plants and traits [40].
In this study, a reference transcriptome of wild germplasm materials of Kentucky bluegrass with signi cant differences in apomictic rates in Gansu, Chain, was constructed using RNA-seq, and the gene expression changes and differences between the two genotypes were studied. However, the ovary of Kentucky bluegrass is very small, and RNA is also very easily degraded, and thus the biological sample we used included the entire spikelet of Kentucky bluegrass, which is composed of 3-4 orets, including the glumes, lemma, palea, lodicules, pistil, and stamen. Therefore, our transcriptome originates from a variety of cell types, including somatic cells and male and female germ cells from the premeiosis to anthesis stages. Based on this, we believe that this database will be very helpful in identifying any transcripts expressed in the owers of Kentucky bluegrass that can be identi ed at a detectable level. However, the spatial and temporal speci city of gene expression will need to be evaluated by other experiments such as qRT-PCR [23].
Regulation of TFs in the ovule development of apomixes: During plant growth and development, the expression of genes in cells is time-speci c and spacespeci c, which leads to differentiation and differences among different cells. One of the main reasons for this phenomenon is the regulation of TFs at the transcriptional level [41]. TFs are protein molecules that can speci cally bind to cis-acting elements in the promoter region of eukaryotic genes to regulate the expression of target genes. They play an important role in plant growth and development, regulation, hormone signal transduction, and stress responses [42]. By comparing the transcripts of obligate apomixis and sexually reproductive Boechera plants, Shah et al. [43] showed that the dysregulation of TFs is one of the main components of apomixis-speci c transcripts. Before this, scientists noted that no obvious developmental pathways or time changes could simply explain the shift from sexual reproduction to apomixis. The overexpression of TFs in apomixis-speci c genes indicates a large-scale change in the regulation of apomixis ovules [44,45]. Sharbel et al. [45] suggested that TFs are important motifs for the differential expression between sexual and apomictic ovules, and their mediation of the inhibition of gene expression is a characteristic of the development of sexual ovules to apomictic ovules. These studies suggest that TFs may regulate the transformation of plants from sexual reproduction to apomixis [46], and thus we further analyzed the differential expression of TFs in the two genotypes of Kentucky bluegrass. The results showed that a total of 84 TF families were differentially expressed in the two genotypes of Kentucky bluegrass, of which 52 were common across all four stages.
Present studies on plant apomixis have reported that many TFs are involved in the regulation of apomixis. Kumar et al. [40] reported a transcriptome analysis of ovule development during nucellar embryogenesis in citrus (Citrus reticulata), and the results showed that MYB and WRKY were key TFs and were differentially expressed in single and polyembryonic ovules as well as between the pre-anthesis and post-anthesis stages of polyembryonic varieties. In addition, there are many reports that TFs such as B3 [47], NAC and bZIP [43], MADS [45], OFP [46], TCP and Al n-like [40] and bHLH, ARF and AUX/IAA [44] can regulate apomixis. These TFs have diverse functions and regulate plant responses to biotic and abiotic stresses in addition to many aspects of plant growth and development. Similarly, in our results, these TFs were differentially expressed in different genotypes and were time-speci c, indicating that they may play an important role in regulating the apomixis development of Kentucky bluegrass. In addition, we also found that TFs such as FAR1, C2H2, TRAF, SWI/SNF-BAF60b, and CPP were differentially expressed in different genotypes. They can regulate different genes and perform a variety of regulatory functions. For example, FAR1 can act as a Cdk inhibitor (CKI) to resist the transformation of G1/S, block cells in the G1 phase and then regulate the cell cycle [48]. It is speculated that they may be involved in the regulation of apomixis in Kentucky bluegrass. However, since our sampling included the whole spikelet, we suspect that many of the TFs we have screened may be involved in ower development, cell division, cell evolution, reproductive organ development, and so forth.
Disorder of meiosis leads to apomixes: The loss or failure of meiosis is one of the three basic events involved in apomixis [15]. Apomixis in Kentucky bluegrass is mainly aposporic, which is usually characterized by the failure of meiosis, followed by differentiation into one or more unreduced embryos from nucellus cells [49]. In Arabidopsis thaliana, apomixis also occurs through important meiotic genes, such as DYAD/SWITCH1 (SWI1), which controls meiotic chromosome tissue. Additionally, the MIME (mitosis instead of meiosis) of a triple mutant, which could be achieved in A. thaliana by combining different combinations of meiotic mutants spo11-1, osd1, and rec8, exhibits characteristics similar to diplospory [50,51]. Although the role of these genes or potential mutations in asexual reproduction in natural apomixis, such as in Boechera, has not been documented, a comparative MMC study between A. thaliana and Boechera indicated that meiotic genes in the MMC of Boechera such, as PARTING DANCERS and SPO11-2 (meiotic recombination protein), may be downregulated compared with those in the MMC of sexual A. thaliana [52]. Interestingly, Okada et al. [53] found that meiotic genes such as DMC1 (DNA meiotic recombinase 1), RAD50 (RAD50 double strand break repair protein), and other genes required for meiosis in A. thaliana, including SWITCH1/DYAD, MULTIPOLAR SPINDLE1, SPOROCYTELESS/NOZZLE, REC8 (meiotic recombination protein) or similar transcripts, were not found in AIC, EAE, or SO cell types. It was previously reported that the meiosis control gene DMC1 was expressed in both sexual and apomixis MMCs in Hieracium spp., but could not be detected in AI cells by situ hybridization [54]. However, the results of this study indicated that these meiosis-related genes were differentially expressed between two genotypes, which may be due to the fact that our samples contained a variety of cellular tissues and originated from multiple time points. In addition, this discussion also implies the difference between diplospory and apospory. The former can be realized through loss of meiosis on the premise of avoiding fertilization, while the ability of nucellus cells to develop into mature embryo sacs through mitosis is more important for apospory.
Regulation of plant hormone signal transduction in apomixis development: Plant hormones directly or indirectly regulate the determination of cell fate, including embryogenesis and postembryonic development, and thus it is speculated that plant hormones may also affect the development of apomixis ovules [55]. In this study, many plant hormone-related genes were differentially expressed between the two genotypes, such as the auxin-related genes SAUR (SAUR-like auxin-responsive protein family), ARF (auxin response factor), IAA (auxin-responsive protein), AUX1 (AUX transcriptional regulator family protein1), cytokinin synthesis-related genes IPT (tRNA dimethylallyltransferase 2), and CKX (cytokinin dehydrogenase 7), indicating that the expression of apospory in Kentucky bluegrass may be accompanied by a response to hormone stimulation. Giulio et al. [16] found that the DEGs between sexual reproduction and apomixis in Hypericum perforatum included genes involved in cytokinin biosynthesis, several auxin response factors, and SAUR-like proteins. After further analysis, it was noted that some genes related to hormone perception and dynamic balance, including cytokinin, auxin, and brassinol, were differentially expressed in the apomictic pistils. Schmidt et al. [52] reported that compared with mature gametophytes, the upregulation of cytokinin degradation genes was detected in apomictic genotypes, while the egg cell was marked by the activity of genes that lead to cytokinin modi cation. At the same time, it was found that the genes related to auxin signal transduction were enriched in AIC, but not expressed in sexual MMC [13]. We also determined the endogenous hormones of two genotypes and found that ZT and IAA may have signi cant effects on reproductive patterns. The change in the regularity of GA 3 and ABA was not obvious. We speculate that the differential expression of genes related to hormone signal transduction disrupts the balance of endogenous hormones in Kentucky bluegrass. The unbalanced hormone content affects plant morphogenesis, and so the two Kentucky bluegrass genotypes exhibit different reproductive patterns. In summary, our analysis of the results shows that the developmental pathway of apospory in Kentucky bluegrass is accompanied by the expression and regulation of multiple genes involved in hormone homeostasis.
Regulation of apomixis in Kentucky bluegrass by embryogenesis-related genes: The successful development of embryos is necessary in both sexual reproduction and apomixis. At present, it has been reported that SERK (somatic embryogenesis receptor-like kinase), LEC (laryngotracheo esophageal cleft), and WUS (WUSCHEL-related homeobox) are the main genes that regulate the development of apomixis. The expression of SERK is closely related to apomixis. Tucker et al. [56] studied the expression of AtSERK1 in sexual and apomictic Hieracium umbellatum, and the differential activity of SERK was also observed between sexual and asexual genotypes of Kentucky bluegrass [25]. In our study, the expression of SERK1 was upregulated in LN during the meiosis stage, while SERK2 was upregulated at the anthesis stage. WUS can promote the transformation from vegetative to embryogenic cells. In the absence of exogenous auxin, it binds with other genes that can induce embryo formation, such as LEC1, LEC2, and SERK1, to regulate the formation and development of embryos [57]. LEC regulates embryo formation by establishing an intercellular environment suitable for embryo development and is an important regulator for inducing embryo morphogenesis and controlling embryo development [58]. Compared with GN, we found that WUS1 and WDR23 (WD40 repeat-containing protein), as the homologous genes of LEC14B, were upregulated and time-speci c. Comprehensive analysis showed that the upregulated expression of genes related to somatic and zygotic embryogenesis may regulate the occurrence of apomixis in Kentucky bluegrass.
Regulation of apomixis by the stress response pathway: One study showed that stress can induce the occurrence of somatic cell embryogenesis [59]. Many studies have observed a close correlation between stress-induced signal pathways and somatic embryogenesis induction [60][61][62]. One of the key steps in apospory is to confer nucellar cells embryogenesis ability. Our results showed that the stress response plays an important role in the development of apomixis in Kentucky bluegrass, and a large number of TFs related to stress responses, such as MYB, WRKY, and NAC, were differentially expressed between the two genotypes. In fact, many studies have found that genes involved in reproduction are associated with responses to stimulation and defense processes [63,64].
Kumar et al. [40] found that many genes related to abiotic stress were signi cantly expressed in citrus polyembryos, including chalcone synthase, heat-shock factor protein, and chaperone protein CLPB1, indicating that the overexpression of stress-related signals was the main characteristic of the pre-anthesis ovules of citrus polyembryonic varieties. Shah et al. [43] also found that apomixis exhibited better tolerance to osmotic stress than sexual reproduction in vitro and was accompanied by the signi cant upregulation of a subset of NAC genes. These results indicate that stress response pathways are involved in early preparatory events before apomixis transition [23]. Previous studies hypothesized that the stress response pathway regulates apomixis due to prolonged UV exposure, which may be responsible for activating meiosis-speci c proteins and initiating double-strand break formation and DNA repair, thus increasing the frequency of recombination [65,66].
Regulation of apomixis by epigenetic pathways and other factors: The genetics of apomixis is quite complex, often with large scale segregation aberrations [67]. Increasing studies have shown that mutations in the apparent pathway lead to phenotypes similar to apomixis phenotypes [68]. Studies on Paspalum notatum and Paspalumsimplex showed that the methylation level of genomic regions controlling apomixis was very high, and it was suggested that the relevant reproductive control factors were inactivated by methylation in Paspalum [69]. APOLLO (DNA cross-link repair 1B protein) is both an apomictic allele and a sexual allele [70], and so this gene has the potential to produce apomixis or sexual reproduction by regulation. In this study, APOLLO was upregulated compared with GN during the meiosis and postmeiosis stages, indicating that this gene regulates the development of apomixis in Kentucky bluegrass. In European rape (Brassica napus), the existence, deletion, and distribution of BnMET1a-like DNA methyltransferase genes are related to methylation, while methylation genes can transform pollen growth and development to embryogenesis, and the resulting embryo and zygotic embryo exhibit similar DNA methylation and MET1a-like expression patterns [71]. Our results also demonstrated that MET1 (methyltransferase 1) was differentially expressed between the two genotypes. Thus it can be seen that the regulation of epigenetics may be a key focus for the continued study of natural apomixis [68]. AGO9 can combine with different types of small RNAs to play a role in the regulation pf apomixis [72], and other genes involved in RNA-guided DNA methylation pathways also play important roles [73]. In parallel, according to a study that reported that small RNAs can promote DNA methylation in plants [74], including in maize, mutations in the effectors of the small RNA pathway have been shown to be important for the reproduction of events similar to diplospory [75]. In addition, some studies have shown that the mechanisms dependent on small RNAs play an important role in the cell speci cation within the ovule [76,77].
Therefore, the regulatory mechanisms of epigenetics on apomixis in Kentucky bluegrass from a small RNA perspective should be explored in the future.
In addition, SPL (squamosa-promoter binding protein-like protein)/NOZZLE, which is involved in the initiation of primary spore differentiation, also plays an important role in the regulation of apomixis [78,79]. In this study, SPL11 and SPL15 were differentially expressed in the two genotypes, indicating that the expression of SPL/NOZZLE regulates the apomixis development of Kentucky bluegrass.

Molecular regulation model of the apomixis mechanism of Kentucky bluegrass
Based on the observation of the histological process of apomixis embryos and the study of gene regulation at the transcriptome level of wild germplasm materials of Kentucky bluegrass, we predicted the molecular regulatory model of apomixis embryogenesis in Kentucky bluegrass (Fig. 10). Compared with sexual reproduction, the apospory of Kentucky bluegrass avoids meiosis and double fertilization, which may be caused by temporal and spatial changes in related gene expression. Therefore, we speculate that the temporal and spatial changes in gene expression lead to the inhibition or imbalance of normal sexual reproduction in Kentucky bluegrass. As the plant needs to ensure the survival of offspring and produce seeds, it promotes the occurrence of apomixis.
In essence, in facultative apomixis species such as Kentucky bluegrass, the occurrence of apomixis is mainly due to the disorder of the sexual reproductive pathway. The differential expression of stress response-related genes is extremely important for the acquisition of embryogenic cells and the change in hormone content. We speculate that the occurrence of events related to the stress response may be a preliminary preparation for apomixis in Kentucky bluegrass,which promotes the embryonic ability of the nucellus cells, causes a hormone imbalance, and leads to a disordered sexual reproductive pathway, which promotes the occurrence and development of apomixis. In summary, the regulatory mechanism of apomixis in Kentucky bluegrass is extremely complex, and its occurrence and development involve the coexpression of genes related to many aspects, such as the imbalance of meiosis, hormone signal transduction, embryonic development, stress response pathway, and epigenetic regulation pathway. In addition, we also found a number of new DEGs that have not been previously reported to be related to apomixis regulation, and so their apomixis regulatory mechanisms require further research.

Conclusion
The molecular basis of apomixis in Kentucky bluegrass was evaluated in detail by transcriptome analysis. In this study, the global gene expression of wild Kentucky bluegrass spikelets exhibiting signi cant differences in apomictic rates was compared, and many genes related to apomixis development were identi ed. The results showed that apomixis in Kentucky bluegrass was caused by changes in the expression of genes related to the sexual reproduction pathway, which led to disordered normal sexual reproduction, resulting in the production of offspring by apomixis. Through the functional enrichment analysis of DEGs, it was found that apomixis was regulated by TFs at the transcriptional level. In addition, it also involves the co-expression and co-regulation of genes related to meiosis, hormone signal transduction, embryonic development, the stress response pathway, and epigenetic pathway. Although this provides a basis for the molecular pathway of apomixis development, further studies on the genomic and functional characteristics are needed to elucidate the nature of its genetic determinants.