Investigation of Differences in Fertility among Progenies from Self-Pollinated Chrysanthemum

Most chrysanthemum cultivars are self-incompatible, so it is very difficult to create pure lines that are important in chrysanthemum breeding and theoretical studies. In our previous study, we obtained a self-compatible chrysanthemum cultivar and its self-pollinated seed set was 56.50%. It was interesting that the seed set of its ten progenies ranged from 0% to 37.23%. Examination of the factors causing the differences in the seed set will lead to an improved understanding of chrysanthemum self-incompatibility, and provide valuable information for creating pure lines. Pollen morphology, pollen germination percentage, pistil receptivity and embryo development were investigated using the in vitro culture method, the paraffin section technique, scanning electron microscopy and transmission electron microscopy. Moreover, RNA sequencing and bioinformatics were applied to analyzing the transcriptomic profiles of mature stigmas and anthers. It was found that the self-pollinated seed set of “Q10-33-1①”,”Q10-33-1③”,”Q10-33-1④” and “Q10-33-1⑩” were 37.23%, 26.77%, 7.97% and 0%, respectively. The differences in fertility among four progenies were mainly attributable to differences in pollen germination percentage and pistil receptivity. Failure of the seed set in “Q10-33-1⑩” was possibly due to self-incompatibility. In the transcriptomic files, 22 potential stigma S genes and 8 potential pollen S genes were found out.


Introduction
The cultivated chrysanthemum is generally allohexaploid and aneuploidy [1,2] with a highly heterozygous genotype, and its trait inheritance and genetic background are extremely complex [3,4] so it is very hard to create chrysanthemum pure lines. However, chrysanthemum pure lines are very important for the utilization of heterosis and genome research. Faced with this, self-fertilization is an important means to develop chrysanthemum pure lines. Usually, pure lines can be created by the multigenerational selfing of self-compatible chrysanthemum cultivars. However, most cultivated chrysanthemum is self-incompatible [1,5]. Therefore, understanding the mechanism of self-incompatibility (SI) is crucial for the development of chrysanthemum pure lines.
Self-incompatibility, also called self-infertility (SI), means the inability of a fully fertile hermaphroditic plant to produce zygotes when self-pollinated [6]. It is distributed widely [7] throughout all the principal lineages of angiosperms, being present in approximately 19 orders, 71 families and 250 genera comprising approximately 60% of all angiosperm species. Therefore, the study of SI is always the hotspot of reproductive developmental biology, and scientists have made major achievements in this field.

Pistil Receptivity of "Q10-33-1"'s Four Progenies
In previous study, we found that average numbers of pollen grains germinating on each stigma among "Q10-33-1"'s ten progenies were nearly positively proportional to their seed set. Among these, the average numbers of pollen grains germinating on each stigma of "Q10-33-1①", "Q10-33-1③", "Q10-33-1④" and "Q10-33-1⑩" were 30.43, 8.57, 5.33 and 7.20, respectively. Thus it could be seen, except for "Q10-33-1⑩," the average numbers of pollen grains germinating on each stigma among other three progenies were positively proportional to their seed set. In this study, the pistil receptivity of "Q10-33-1"s four progenies was observed at the ultrastructural level ( Figure 4). It was found that pollen tube stopped growing and failed to penetrate the stigma's surface in "Q10-33-1⑩" ( Figure  4H). This indicated that abnormal pistil receptivity was the main reason leading to the failure of seed set in "Q10-33-1⑩."

Ovary Development of "Q10-33-1"'s Four Progenies
In all four progenies, the percentage of full ovaries were nearly negatively proportional to the days after self-pollination (Table 3), indicating that ovary abortion occurred gradually with the extension of the days after self-pollination. Meanwhile, the percentage of full ovaries among four progenies were nearly positively proportional to their seed set (Table 3), showing that ovary development of "Q10-33-1"'s four progenies was directly related to their seed set. Among the four progenies, no full ovary was found in "Q10-33-1⑩," and pollen tube growth of this progeny was abnormal, these might come to a conclusion that "Q10-33-1⑩" was self-incompatible.

Pistil Receptivity of "Q10-33-1"'s Four Progenies
In previous study, we found that average numbers of pollen grains germinating on each stigma among "Q10-33-1"'s ten progenies were nearly positively proportional to their seed set. Among these, the average numbers of pollen grains germinating on each stigma of "Q10-33-1", "Q10-33-1", "Q10-33-1" and "Q10-33-1" were 30.43, 8.57, 5.33 and 7.20, respectively. Thus it could be seen, except for "Q10-33-1," the average numbers of pollen grains germinating on each stigma among other three progenies were positively proportional to their seed set. In this study, the pistil receptivity of "Q10-33-1"s four progenies was observed at the ultrastructural level ( Figure 4). It was found that pollen tube stopped growing and failed to penetrate the stigma's surface in "Q10-33-1" ( Figure 4H). This indicated that abnormal pistil receptivity was the main reason leading to the failure of seed set in "Q10-33-1."

Ovary Development of "Q10-33-1"'s Four Progenies
In all four progenies, the percentage of full ovaries were nearly negatively proportional to the days after self-pollination (Table 3), indicating that ovary abortion occurred gradually with the extension of the days after self-pollination. Meanwhile, the percentage of full ovaries among four progenies were nearly positively proportional to their seed set ( Table 3), showing that ovary development of "Q10-33-1"'s four progenies was directly related to their seed set. Among the four progenies, no full ovary was found in "Q10-33-1," and pollen tube growth of this progeny was abnormal, these might come to a conclusion that "Q10-33-1" was self-incompatible.  Seed set were quoted from our previous study [23].
In addition, the differences between full and unful ovaries on different days after self-pollination were observed under a stereoscope. The ovule of a full ovary was also full and gradually became dark in color until seed formation, while the ovule of an unful ovary was smaller, the color was white, and it had no tendency to form seed ( Figure 5).  In addition, the differences between full and unful ovaries on different days after self-pollination were observed under a stereoscope. The ovule of a full ovary was also full and gradually became dark

Embryo Development of "Q10-33-1"'s Three Progenies at Microscopic Level
If the process of embryo development was normal, heart embryos could be seen in ovaries at 12 days after self-pollination, and torpedo or cotyledon embryos could be seen in ovaries at 18 days after self-pollination ( Figure 6A-C). However, embryo abortion occurred in some ovaries with the extension of the days after self-pollination ( Figure 6D,E). Meanwhile, no embryos could be seen in "Q10-33-1⑩", we could only find tetranucleate embryo sacs, showing one central cell, one egg cell, two synergids near the micropylar end and two degenerated megaspores ( Figure 6F,G).

Embryo Development of "Q10-33-1"'s Three Progenies at Microscopic Level
If the process of embryo development was normal, heart embryos could be seen in ovaries at 12 days after self-pollination, and torpedo or cotyledon embryos could be seen in ovaries at 18 days after self-pollination ( Figure 6A-C). However, embryo abortion occurred in some ovaries with the extension of the days after self-pollination ( Figure 6D,E). Meanwhile, no embryos could be seen in "Q10-33-1", we could only find tetranucleate embryo sacs, showing one central cell, one egg cell, two synergids near the micropylar end and two degenerated megaspores ( Figure 6F,G).

Embryo Development of "Q10-33-1"'s Three Progenies at Microscopic Level
If the process of embryo development was normal, heart embryos could be seen in ovaries at 12 days after self-pollination, and torpedo or cotyledon embryos could be seen in ovaries at 18 days after self-pollination ( Figure 6A-C). However, embryo abortion occurred in some ovaries with the extension of the days after self-pollination ( Figure 6D,E). Meanwhile, no embryos could be seen in "Q10-33-1⑩", we could only find tetranucleate embryo sacs, showing one central cell, one egg cell, two synergids near the micropylar end and two degenerated megaspores ( Figure 6F,G).  (D, E) Embryo abortion of "Q10-33-1"; (D) Embryo abortion at 12 days after self-pollination, the embryo was degenerating; (E) Embryo abortion at 18 days after self-pollination, the embryo was degenerating and the integument was thickening; (F, G) Unfertilized ovules of "Q10-33-1", tetranucleate embryo sacs, two black arrows indicated two degenerated megaspores. Abbreviations: Cc: central cell, Ec: egg cell, Sy: synergid.

Embryo Development of "Q10-33-1" at Ultrastructural Level
The results of TEM showed that there were significant differences between normal and abortive embryonic cells (Figure 7). In normal embryonic cells, the nucleus and nucleolus could be seen clearly and some typical cell organelles such as mitochondria, Golgi complex and vacuoles were also found ( Figure 7A,B). Meanwhile, at the early stage of embryo abortion, cytoplasmic vacoulation, condensed cytoplasm and degradative organelles were observed in abortive embryonic cells ( Figure 7C,D). At late stage of embryo abortion, the nucleus and organelles were fully degraded. No clear cellular tissues could be seen in abortive embryonic cells ( Figure 7E,F).

Embryo Development of "Q10-33-1①" at Ultrastructural Level
The results of TEM showed that there were significant differences between normal and abortive embryonic cells (Figure 7). In normal embryonic cells, the nucleus and nucleolus could be seen clearly and some typical cell organelles such as mitochondria, Golgi complex and vacuoles were also found ( Figure 7A,B). Meanwhile, at the early stage of embryo abortion, cytoplasmic vacoulation, condensed cytoplasm and degradative organelles were observed in abortive embryonic cells ( Figure 7C,D). At late stage of embryo abortion, the nucleus and organelles were fully degraded. No clear cellular tissues could be seen in abortive embryonic cells ( Figure 7E,F).

Transcriptome Sequencing and Read Assembly
To study transcriptomic gene expression differences between mature stigmas of self-compatible and self-incompatible progenies, two cDNA libraries from 3S (mature stigmas of "Q10-33-1③") and 10S (mature stigmas of "Q10-33-1⑩") were subjected to Illumina sequencing. After reads filtering, we obtained 44.52 Mb and 44.19 Mb clean reads from two cDNA libraries, containing 6.68 Gb and 6.63 Gb clean bases, respectively. The clean reads ratios of the two samples were about 80%, the Q20 percentages were more than 97% and the Q30 percentages were more than 91% (Table 4). After clustering the high-quality reads, we finally obtained 113,800 unigenes with a total length of 94,376,715 bp, a mean length of 829 bp and a GC percentage of 39.40% (Table 5).
To study transcriptomic gene expression differences between mature anthers of self-compatible and self-incompatible progenies, two cDNA libraries from 3A (mature anthers of "Q10-33-1③") and 10A (mature anthers of "Q10-33-1⑩") were subjected to Illumina sequencing. After reads filtering, we obtained 44.25 Mb and 44.37 Mb clean reads from two cDNA libraries, containing 6.64 Gb and 6.65 Gb clean bases, respectively. The clean reads ratios of the two samples were about 80%, the Q20 percentages were more than 97% and the Q30 percentages were more than 91% (Table 6). After

Transcriptome Sequencing and Read Assembly
To study transcriptomic gene expression differences between mature stigmas of self-compatible and self-incompatible progenies, two cDNA libraries from 3S (mature stigmas of "Q10-33-1") and 10S (mature stigmas of "Q10-33-1") were subjected to Illumina sequencing. After reads filtering, we obtained 44.52 Mb and 44.19 Mb clean reads from two cDNA libraries, containing 6.68 Gb and 6.63 Gb clean bases, respectively. The clean reads ratios of the two samples were about 80%, the Q20 percentages were more than 97% and the Q30 percentages were more than 91% (Table 4). After clustering the high-quality reads, we finally obtained 113,800 unigenes with a total length of 94,376,715 bp, a mean length of 829 bp and a GC percentage of 39.40% (Table 5).
To study transcriptomic gene expression differences between mature anthers of self-compatible and self-incompatible progenies, two cDNA libraries from 3A (mature anthers of "Q10-33-1") and 10A (mature anthers of "Q10-33-1") were subjected to Illumina sequencing. After reads filtering, we obtained 44.25 Mb and 44.37 Mb clean reads from two cDNA libraries, containing 6.64 Gb and 6.65 Gb clean bases, respectively. The clean reads ratios of the two samples were about 80%, the Q20 percentages were more than 97% and the Q30 percentages were more than 91% (Table 6). After clustering the high-quality reads, we finally obtained 113,638 unigenes with a total length of 91,284,692 bp, a mean length of 803 bp and a GC percentage of 39.66% (Table 7).    N50: a weighted median statistic that 50% of the total length is contained in unigenes great than or equal to this value. GC (%): the percentage of G and C bases in all unigenes.

Unigene Functional Annotation
To predict the potential functions of unigenes in chrysanthemum stigmas, all assembled unigenes were subjected to functional annotation using 7 functional databases, including NR, NT, Swiss-Prot, KEGG, COG, Interpro and GO. Consequently, 61,731 unigenes were annotated with the databases, and the number of unigenes annotated with each database was 57,127,42,115,39,577,42,586,19,476,38,563 and 23,423, respectively (Table 8).
According to the NR annotation, we used GO assignments to classify the functions of unigenes. As a result, 23,423 unigenes were classified into 55 functional subcategories, which belonged to 3 main categories: biological process, cellular component and molecular function. There were 23, 17 and 15 functional groups in each category, respectively (Figure 9). In the biological process category, "metabolic process" and "cellular process" were the main functional subcategories. In terms of cellular component category, two major subcategories were "cell" and "cell part". For the molecular function category, "catalytic activity" and "binding" were remarkable. Finally, KEGG annotation was performed to identify the biological pathways activated in chrysanthemum's mature stigmas. A total of 42,586 unigenes were assigned to 135 KEGG pathways (Table S3). The majority of these pathways were "metabolic pathways (ko01100)" (21.41%), "biosynthesis of secondary metabolites (ko01110)" (12.27%), "plant-pathogen interaction (ko04626)" (4.02%) and "endocytosis (ko04144)" (3.69%).
To predict potential functions of unigenes in chrysanthemum anthers, all assembled unigenes were performed functional annotation with 7 functional databases, including NR, NT, Swiss-Prot, KEGG, COG, Interpro and GO. Consequently, 63,517 unigenes were annotated with the databases, According to the NR annotation, we used GO assignments to classify the functions of unigenes. As a result, 23,423 unigenes were classified into 55 functional subcategories, which belonged to 3 main categories: biological process, cellular component and molecular function. There were 23, 17 and 15 functional groups in each category, respectively (Figure 9). In the biological process category, "metabolic process" and "cellular process" were the main functional subcategories. In terms of cellular component category, two major subcategories were "cell" and "cell part". For the molecular function category, "catalytic activity" and "binding" were remarkable. and the number of unigenes annotated with each database was 59,346, 42,785, 40,774, 43,970, 20,580, 40,653 and 24,313, respectively (Table 9).  Overall: the number of unigenes which are annotated with at least one functional database. Finally, KEGG annotation was performed to identify the biological pathways activated in chrysanthemum's mature stigmas. A total of 42,586 unigenes were assigned to 135 KEGG pathways (Table S3). The majority of these pathways were "metabolic pathways (ko01100)" (21.41%), "biosynthesis of secondary metabolites (ko01110)" (12.27%), "plant-pathogen interaction (ko04626)" (4.02%) and "endocytosis (ko04144)" (3.69%).
According to the NR annotation, we used GO assignments to classify the functions of unigenes. As a result, 24,313 unigenes were classified into 54 functional subcategories, which belonged to 3 main categories: biological process, cellular component and molecular function. There were 22, 17 and 15 functional groups in each category, respectively ( Figure 11). In the biological process category, "metabolic process" and "cellular process" were the main functional subcategories. In terms of cellular component category, two major subcategories were "cell" and "cell part". For the molecular function category, "catalytic activity" and "binding" were remarkable.
Finally, KEGG annotation was performed to identify the biological pathways activated in chrysanthemum's mature anthers. A total of 43,970 unigenes were assigned to 135 KEGG pathways (Table S4). The majority of these pathways were "metabolic pathways (ko01100)" (21.36%), "biosynthesis of secondary metabolites (ko01110)" (12.16%), "endocytosis (ko04144)" (4.14%) and "plant-pathogen interaction (ko04626)" (3.95%).  According to the NR annotation, we used GO assignments to classify the functions of unigenes. As a result, 24,313 unigenes were classified into 54 functional subcategories, which belonged to 3 main categories: biological process, cellular component and molecular function. There were 22, 17 and 15 functional groups in each category, respectively ( Figure 11). In the biological process category, "metabolic process" and "cellular process" were the main functional subcategories. In terms of cellular component category, two major subcategories were "cell" and "cell part". For the molecular function category, "catalytic activity" and "binding" were remarkable.

DEGs (Differently Expressed Genes) Related to SI in Mature Stigmas and Anthers of Self-Compatible and Self-Incompatible Progenies
Based on the FPKM (fragments per kb per million fragments) value, we explored gene expression level in 3S and 10S. Compared with 3S, in 10S, 8130 and 3766 genes were up-and down-regulated, respectively. The details about these DEGs were shown in Table S5. To determine putative functions of the DEGs, GO and KEGG analysis were applied. As a consequence, 7483 DEGs were associated with GO categories, 7194 DEGs were mapped to 134 KEGG pathways (Table S6). The high number of DEGs provided an abundant list of candidate SI-related genes in stigmas.
Based on the FPKM value, we explored gene expression level in 3A and 10A. Compared with 3A, in 10A, 7287 and 4538 genes were up-and down-regulated, respectively. The details about these DEGs were shown in Table S7. To determine putative functions of the DEGs, GO and KEGG analysis were applied. As a consequence, 7657 DEGs were associated with GO categories, 7020 DEGs were mapped to 133 KEGG pathways (Table S8). The high number of DEGs provided an abundant list of candidate SI-related genes in anthers. Finally, KEGG annotation was performed to identify the biological pathways activated in chrysanthemum's mature anthers. A total of 43,970 unigenes were assigned to 135 KEGG pathways (Table S4). The majority of these pathways were "metabolic pathways (ko01100)" (21.36%), "biosynthesis of secondary metabolites (ko01110)" (12.16%), "endocytosis (ko04144)" (4.14%) and "plant-pathogen interaction (ko04626)" (3.95%).

DEGs (Differently Expressed Genes) Related to SI in Mature Stigmas and Anthers of Self-Compatible and Self-Incompatible Progenies
Based on the FPKM (fragments per kb per million fragments) value, we explored gene expression level in 3S and 10S. Compared with 3S, in 10S, 8130 and 3766 genes were up-and down-regulated, respectively. The details about these DEGs were shown in Table S5. To determine putative functions of the DEGs, GO and KEGG analysis were applied. As a consequence, 7483 DEGs were associated with GO categories, 7194 DEGs were mapped to 134 KEGG pathways (Table S6). The high number of DEGs provided an abundant list of candidate SI-related genes in stigmas.
Based on the FPKM value, we explored gene expression level in 3A and 10A. Compared with 3A, in 10A, 7287 and 4538 genes were up-and down-regulated, respectively. The details about these DEGs were shown in Table S7. To determine putative functions of the DEGs, GO and KEGG analysis were applied. As a consequence, 7657 DEGs were associated with GO categories, 7020 DEGs were mapped to 133 KEGG pathways (Table S8). The high number of DEGs provided an abundant list of candidate SI-related genes in anthers.
After deep analysis and screening of these DEGs, some potential stigma S genes and pollen S genes were found out, they were listed in Tables 10 and 11, respectively.

Validation of Gene Expression Profiles by qRT-PCR
In the library of stigmas, ten differentially expressed genes were selected for qRT-PCR to test the reliability of RNA-seq data. The result showed that the expression trend of almost all genes was consistent with sequencing data (Figure 12). Most of DEGs selected were related to SI female determinants in other plant families, such as S-receptor-like serine/threonine-protein kinase, protein in ribonuclease T2 family and epidermis-specific secreted glycoprotein.
In the library of anthers, ten differentially expressed genes were also selected for qRT-PCR to test the reliability of RNA-seq data. The result showed that the expression trend of almost all genes was consistent with sequencing data (Figure 13). Most of DEGs selected were related to pollen germination, pollen tube reception and growth, and SI male determinants in other plant families, such as pollen coat protein.

Cellular Mechanism of Differences in Fertility among Progenies
This study shows that the self-pollinated seed set of chrysanthemum is mainly affected by pollen germination percentage, the germination behavior of pollen grains on stigmas, pollen tube growth and embryo development.
Pollen germination percentage has some effects on seed set [24,25], while abnormal pollen morphology can cause low pollen germination percentage. In "Q10-33-1"s four progenies, the number of abnormal pollen grains in "Q10-33-1" was the largest, so its pollen germination percentage was the lowest; while the number of abnormal pollen grains in "Q10-33-1" was the least, so its pollen germination percentage was the highest. In addition, except for "Q10-33-1", pollen germination percentages of other three progenies were positively proportional to their seed set. This indicated that pollen germination percentage indeed had an effect on seed set. However, seed set of "Q10-33-1" was not connected with its pollen germination percentage. Its pollen germination percentage was the highest in four progenies, but its seed set was 0. This indicated that "Q10-33-1"s pollen germination percentage was not the main reason for its self-pollinated failure.
In previous study, we found that the interaction between some chemical products secreting from stigmas and those secreting from pollen walls was incompatible, and consequently prevented pollen grains from germinating normally on stigmas, resulting in low seed set [26]. Therefore, the germination behavior of pollen grains on stigmas also has an effect on seed set [27][28][29]. In "Q10-33-1"s four progenies, except for "Q10-33-1", the average numbers of pollen grains germinating on each stigma among other three progenies were positively proportional to their seed set. This result just confirmed the above view.
Lastly, embryo development is another important factor related to seed set [32][33][34]. Except for "Q10-33-1", embryo abortion was observed in other three progenies, this indicated that embryo development also really had an effect on seed set.
In conclusion, seed set of: "Q10-33-1", "Q10-33-1", and "Q10-33-1" were related to pollen germination percentage, the germination behavior of pollen grains on stigmas and embryo development. Meanwhile, seed set of "Q10-33-1" was mainly affected by pollen tube growth. We also came to a conclusion that "Q10-33-1" was self-incompatibility. In addition, as shown in the Figure 4H, the inhibition of pollen tube growth in "Q10-33-1" was taken place on the surface of stigma, rather than in the style. This indicated that SI in chrysanthemum was more likely to SSI.

Molecular Mechanism of SI
In molecular study, SI of angiosperm is mainly controlled by allele-specific interactions between two highly polymorphic S-determinants. So far, however, little research has been done on chrysanthemum SI. In our transcriptome data, some DEGs related to S-determinants in other plant families were found, and they might also play an important role in chrysanthemum SI.
According to a previous study, GSI was found typically in Solanaceae, Rosaceae, Scrophulariaceae and Papaveraceae, whereas SSI was found in Brassicaceae, Asteraceae and Convolvulaceae [8]. The most studied SSI system is Brassicaceae SSI. In this family, the female S-determinant is SRK (S-receptor serine/threonine kinase), moreover, SLG (S-locus glycoprotein) also plays a supporting role in SI. The male S-determinant is SCR (S-locus cysteine-rich protein)/SP11 (S-locus protein 11), which is a kind of pollen coat protein [35]. This may provide a key clue for analyzing SI of chrysanthemum. In the analysis of stigma transcriptome, some unigenes encoding S-receptor-like serine/threonine-protein kinase or S-locus glycoprotein were up-regulated in 10S, such as Unigene8775 (4.3 in 10S, 0.84 in 3S), Unigene12923 (8.86 in 10S, 3.51 in 3S) and so on. Meanwhile, in the analysis of anther transcriptome, some unigenes encoding pollen coat-like protein were up-regulated in 10A, such as CL10523.Contig1 (42.99 in 10A, 0.96 in 3A), Unigene12281 (1008.08 in 10A, 217.32 in 3A) and so on. That is, the unigenes related to S-determinants in Brassicaceae were up-regulated in self-incompatible chrysanthemum progeny. The results raise a hypothesis that SI system of chrysanthemum is similar to Brassicaceae SSI system.
Although it was reported that SSI was found in Asteraceae, there was no precise evidence to confirm it. Maybe GSI system of Rosaceae also plays an important role in chrysanthemum's SI. In this GSI system, the female S-determinant is S-RNase (S ribonuclease), belonging to ribonuclease T2 family, and the male S-determinant is SFB/SLF (S-locus F-box protein) [36][37][38]. In the comparative analysis of stigma transcriptome, some unigenes encoding ribonuclease belonging to ribonuclease T2 family were up-regulated in 10S, such as CL9204.Contig2 . The three kinds of proteins, F-box, Skp1 and Cullin1, can form a SCF complex (Skp, Cullin, F-box containing complex), this complex plays a pivotal role in S-RNase-based SI system [36]. Likewise, these results raise a hypothesis, that is, S-RNase-based GSI system may also function in chrysanthemum SI.
Furthermore, in stigma transcriptome analysis, some unigenes participating in recognition of pollen were found, such as CL1851.Contig3 and CL2175.Contig6. We also found some unigenes encoding specific proteins which had been reported in stigma transcriptome files of other plants [39,40]. For instance, CL9691.Contig3 and Unigene23027 encoded stigma-specific peroxidase, CL831.Contig1 and Unigene47365 encoded mitogen-activated protein kinase (MAPK), and CL916.Contig6 and Unigene590 encoded pistil-specific extensin-like protein. Meanwhile, on the basis of previous studies [41,42], in another transcriptome analysis, some unigenes participating in pollen germination, pollen exine formation, pollen tube reception and pollen tube growth were found. For example, CL500.Contig1 and CL3143.Contig2 participated in pollen germination, CL12767.Contig2 took part in pollen exine formation, CL6315.Contig9 participated in pollen tube reception, and CL9992.Contig1 took part in pollen tube growth. The above unigenes may also play an important role in chrysanthemum SI.
From the above, we surmise that in chrysanthemum SI, the SSI system-similar to Brassiceae SI-plays a main role, and S-RNase-based SI system may play a supporting role.

Determination of Flowering Traits
On the basis of relative studies of other researchers [43][44][45], we selected some well-grown inflorescences and measured their flowering traits, including flower color, petal shape, anthotaxy diameter, flower disc diameter, length and width of ray floret, and the number of ray florets and tubiform florets. Each trait was determined five times. In addition, we observed and photographed the inflorescence, tubiform floret and ray floret of it under a SLR (single lens reflex) camera and a stereoscope.

Ultrastructural Observation of Pollen Morphology
Abnormal pollen morphology may lead to low pollen germination percentage, thus affecting seed set. In order to support this view, fresh pollen grains were collected between 10:00 and 14:00 on a sunny day. The pollen grains were dried at 42 • C for 10-15 days, until they were not together [46][47][48]. Then the pollen grains were coated with gold, and observed under a scanning electron microscope.

Determination of Pollen Germination Percentage
We adopted in vitro culture method [26,49] to determine pollen germination percentage. A few drops of culture medium ME3 (0.1 g·L −1 H 3 BO 3 , 0.3 g·L −1 CaCl 2 ·2H 2 O, 0.001 g·L −1 CoCl 2 ·6H 2 O) + 200 g·L −1 PEG4000 were added on a glass slide, then fresh pollen grains were distributed on the glass slide. After that, the glass slide was incubated at 20 • C for about two hours, because it was found from the past research that pollen grains would not germinate after incubating two hours. Germinated pollen grains were observed under an Olympus BX41 microscope (Olympus Corporation, Tokyo, Japan) and germination percentage was counted from ten optical fields. When pollen tube length was longer than the diameter of the pollen grain, the pollen grain was identified as germinating. Each experiment was conducted three times. Notably, pollen germination percentage must be determined immediately after pollen grains collection, because as time goes on, the germination percentage will reduce rapidly.

Ultrastructural Observation of Pistil Receptivity
Five self-pollinated inflorescences were cut off from each progeny at 24 h after pollination, immediately fixed in 2.5% glutaraldehyde solution, and put under vacuum until submerged. Then they were stored at 4 • C until use. The pistils were dissected from the ray florets, dehydrated in an ethanol series, subjected to critical-point drying, and coated with gold. Samples were observed with a scanning electron microscope.

Examination of Embryo Development
Five inflorescences at 12, 18 and 30 days after self-pollination were respectively cut off from each progeny. Then full and unful ovaries of these inflorescences were separately counted and photographed under a stereoscope. Each inflorescence was a repetition. Meanwhile, about five inflorescences of "Q10-33-1", "Q10-33-1" and "Q10-33-1"-at 12 and 18 days after self-pollination-were respectively collected, immediately fixed in FAA solution (95% ethanol:distilled water: formalin:glacial acetic acid = 63:27:5:5), put under vacuum until submerged and then stored at 4 • C until use. Full and unful ovaries were separately dissected from florets, dehydrated in an ethanol series, and then embedded in paraffin wax. Sections were cut to a thickness of 8-10 µm and stained with Heidenhain's hematoxylin (Heidenhain, Berlin, Germany) [50,51]. The sections then were observed and photographed under an Olympus BX41 microscope. In addition, five inflorescences of "Q10-33-1" at 12 and 18 days after self-pollination were collected, fixed in 2.5% glutaraldehyde solution, put under vacuum until submerged, and stored at 4 • C until use. Full and unful ovules were separately dissected from ovaries, dehydrated in an ethanol series, subjected to critical-point drying, and coated with gold. Samples were examined with a transmission electron microscope, and all images were processed digitally.
Total RNA was extracted using the Trizol reagent according to the manufacturer's protocol (Takara Bio Inc., Otsu, Japan) [53]. The quantity and quality of total RNA were confirmed by Agilent 2100 RNA 6000 Kit (Agilent Technologies Inc., Santa Clara, CA, USA) and electrophoresis on 1% agarose gel [53,54].

cDNA Library Construction and Illumina Sequencing
Illumina sequencing was performed at Beijing Genomics Institute (BGI)-Shenzhen, China according to the manufacturer's instructions (Illumina, San Diego, CA, USA). Each cDNA sequencing library was constructed with the RNA mixture from three biological replicates.

Sequencing Reads Filtering and De Novo Assembly
The sequencing reads which containing low-quality, adaptor-polluted and high content of unknown base reads, should be processed to remove these reads before downstream analysis [55].
After filtering, we obtained clean reads from raw data, then Trinity [56] was applied to perform de novo assembly with clean reads. Next, Tgicl [57] was used to cluster transcripts to unigenes. The transcriptome datasets are available in the NCBI repository, (https://submit.ncbi.nlm.nih.gov/ subs/sra/?from) Accession No. for library SRP131835.

Unigene Functional Annotation
Unigene sequences were aligned by Blastx to protein databases such as NR (NCBI nonredundant protein), Swiss-Prot, KEGG (Kyoto Encyclopedia of Genes and Genomes) and COG (Cluster of Orthologous Groups) and aligned by Blastn to the nucleotide database NT (E-value < 0.00001). In addition, we used Blast2 GO program [58] to obtain GO (Gene Ontology) [41] annotation of unigenes. GO has three ontologies describing molecular function, cellular component and biological process. After that, WEGO software [59] was used to GO functional classification for all unigenes. Finally, with the KEGG database, we could further study the complicated biological behaviors of unigenes and get pathway annotation of them.

Analysis of Differences in Unigene Expression
To predict unigene expression levels in different samples, we calculated the unigene expression using the FPKM (fragments per kb per million fragments) method [60]. This method can eliminate the influence of different gene lengths and sequencing levels, so the calculated unigene expression can be directly used for comparing unigene expression differences between samples. After calculation, differentially expressed genes (DEGs) were detected between two samples by using the method of Audic and Claverie [61]. In our analysis, DEGs should meet the criteria that FDR (false discovery rates) ≤ 0.001 and fold change ≥ 2, then GO functional analysis and KEGG pathway analysis were carried out for DEGs.
In GO functional analysis, all DEGs were mapped to each term of Gene Ontology database (http://www.geneontology.org/) and the gene number of each GO term was counted. After that, hypergeometric test was used to find significantly enriched GO terms in DEGs. In KEGG pathway analysis, significantly enriched metabolic pathways or signal transduction pathways were identified in DEGs.

Quantitative Real-Time PCR (qRT-PCR) Validation
To ensure the libraries' reliability, we randomly selected 10 DEGs from each library and validated the data by qRT-PCR with three biological replicates. Gene specific primers were designed using Primer Premier 5.0 (Premier Biosoft International, Palo Alto, CA, USA) and the reference gene was the chrysanthemum EF1α (Elongation Factor 1a) gene (Genbank accession number KF305681). The sequences of all primers including those of EF1α were listed in Tables S1 and S2. The amplification with these primers was specific and the efficiency of qPCR was about 100%. qRT-PCR was performed on the Mastercycler ep realplex device (Eppendorf, Hamburg, Germany). Each 20 µL qPCR mixture contained 10 µL SYBR Premix Ex Taq II (Takara), 1 µL forward primer, 1 µL reverse primer, 3 µL ddH 2 O and 5 µL cDNA template. The cDNA was synthesized from 1 µg total RNA with PrimeScript ® Reverse Transcriptase (Takara), following the manufacturer's instructions. The quality of RNA met the criterion that the values of A260/A280 were around 2.0 and the values of A260/A230 were all larger than 2.0. The PCR cycling regime was initially denatured (95 • C/2 min), followed by 40 cycles of 95 • C/10 s, 55 • C/15 s and 72 • C/20 s. A melting curve was followed each assay to confirm the amplicons' specificity. Relative expression levels of the DEGs were calculated by 2 −∆∆Ct method [41,62].

Conclusions
We performed a series of cellular experiments and RNA-seq to investigate the differences in fertility among progenies from self-pollinated chrysanthemum. Finally, we concluded that seed set of "Q10-33-1", "Q10-33-1" and "Q10-33-1" were mainly determined by pollen germination percentage, while self-pollinated failure of "Q10-33-1" was determined by pistil receptivity. That is to say, "Q10-33-1" was self-incompatible. Moreover, we guessed both SSI system similar to Brassicaceae SI and S-RNased-based GSI system functioned in chrysanthemum SI. According to abnormal pollen tube in Figure 4H and previous studies [7,63], we came to a conclusion that SSI system played a crucial role. These results will help us to produce self-compatible chrysanthemum cultivars in future, which can further help us to create chrysanthemum pure lines.