Cytological and transcriptional insights of late-acting self-incompatibility in tea plants (Camellia sinensis)

Self-incompatibility (SI) is an important genetically controlled mechanism that inhibits inbreeding and facilitates outcrossing for many hermaphrodite species^{[1, 2]}. According to different modes of genetic control of pollen incompatibility phenotype, SI was classified into two types, gametophytic self-incompatibility (GSI) and sporophytic self-incompatibility (SSI)^{[1, 3]}. The GSI mechanism in Rosaceae, Solanaceae and Scrophulariaceae, the SI reaction is dominated by the interaction between the female determinant S-RNase and the male component S-locus F-box protein (SLF/SFB)^{[2, 4]}. The GSI system of the Papaveraceae, where the Ca²⁺ signal cascade leads to programmed cell death (PCD) in the pollen tube (PT)^[5]. During SSI in Brassicaceae, the interaction between a stigma-specific S-locus receptor kinase (SRK) and a pollen-specific S-locus cysteine-rich protein (SCR)^{[6, 7]}. In some plants including Ipomopsis aggregata, Citrus reticulate and Cyrtanthus breviflorus, pollen germination and pollen tube growth are not inhibited, and the pollen tube can pass through the style to reach the ovary after self-pollination, while it still appears as self-incompatible, which is called late-acting self-incompatibility (LSI) or ovarian self-incompatibility (OSI)^[8−10].

As a barrier to fertilization, SI can maintain the stability and genetic variability of different plant species. However, SI makes it difficult to cross-pollinate, self-pollinate and cultivate inbred lines in the plant breeding process for breeders, thereby various methods have been conducted to overcome SI. For instance, the application of flower organ extracts to stigmas and the temperature treatment of pollen to overcome SI in Lilium longiflorum^[11], modification of bud pollination to breakdown SI during pistil development in Lycopersicon peruvianum^[12], overcome SI through S-Rnase gene knockout using CRISP-Cas9 in diploid potato ^[13]. Among them, breeding of plants having SI through bud pollination or delayed pollination is one of the most effective methods. Bud pollination is to pollinate fertile pollen on the stigma of flower buds 2−7 d before flowering, and delayed pollination is to pollinate fertile pollen on the stigma of old pistil 3−9 d after flowering^[14−16].

The tea plant (Camellia sinensis) has been cultivated for thousands of years, it is one of the three popular nonalcoholic beverage crops consumed worldwide^{[17, 18]}. Although the tea plant has hermaphrodite flowers, it is a self-incompatible species where pollination with pollen of the same genotype will cause fertilization failure^[19−21]. Evidence has shown that LSI is caused by the suppression of pollen tubes in the ovaries after self-pollination in tea plants^[22], and LSI was initiated in style and sustained up to ovary with the high expression of SRKs, CsRNS, and SKIPs during self-pollination^[23]. In styles after self- and cross-pollination, several SCF complexes and a putative S-RNase gene were identified based on transcriptome data, indicate that SI maybe under gametophytic control^{[24, 25]}. In SI system of tea plant, Ca²⁺ and K⁺ are involved in signal transduction, and LecRLK and UGT74B1 may play a role in controlling SI^[26]. Nevertheless, the molecular mechanism of tea plant SI system remains unclear, severely limiting the technology development for overcoming the SI to obtain inbred germplasm.

In the present study, we found that self-pollination at bud stage and full bloom stage showed lower SI, while at initial bloom stage showed high SI in tea plants. To understand the mechanism of SI, paraffin-embedded sections and aniline blue staining were performed to observe the pollen tube elongation in styles and ovaries after self-pollination at different flowering stages. Subsequently, the transcriptomes of self-pollinated pistils at the three stages identified many DEGs that may be related to the SI system. This study provides a theoretical basis for breeding of inbred lines and directional selection of excellent varieties in tea plant.

Six-year-old clonal tea cultivars (Camellia sinensis 'Shuchazao') were selected from the Tea Plant Cultivar and Germplasm Resource Garden, Anhui Agricultural University (China) (31°25' N, 117°09' E). Pollen were collected from balloon stage flowers, and then placed in a desiccator for 12 h at room temperature. In order to ensure good viability of the pollen, a viability test on the pollen in vitro was conducted. The different tissues were collected from unpollinated pistils: filaments, styles, ovaries, calyxes, anthers, petals and leaves. The samples were frozen immediately in liquid nitrogen and stored at −80 °C in a refrigerator.

Flowers in balloon stage were bagged after emasculation and self-pollinated on sunny days in November. The pistils were harvested at different time points (12, 24, 36, 48, and 60 h) after self-pollination and fixed in carnoy fixative solution (Phygene Life Sciences Company, China). Furthermore, flowers at different flowering stages were bagged and self-pollinated. Six pollination treatments were designed with three biological replicates, including unpollinated pistil at bud stage (CKB), unpollinated pistil at initial bloom stage (CKIB), unpollinated pistil at full bloom stage (CKFB), pistil at 48 h after self-pollination (HASP) at bud stage (BP), pistil at 48 HASP at initial bloom stage (IBP), and pistil at 48 HASP at full bloom stage (FBP) after emasculation. After harvesting, some of the pistils of each sample were fixed in carnoy fixative solution and the remaining pistils were stored at −80 °C.

Cytological observation

The pistil fixed in carnoy fixative solution for 3 h was cut along the base of the style to separate the style and ovary. After dissecting the style, it was placed carefully in 70% ethanol and stored at 4 °C for more than 2 h. The style was washed by 50% ethanol, 20% ethanol and distilled water for 5 min. Then, the style was stained in 0.5% aniline blue solution that dissolved in 0.3 M K₃PO₄ for 5 min. The style was placed in 0.5% aniline blue solution on a glass slide, and squashed under a cover slip to ensure the style spread out evenly. Olympus IX73 fluorescence microscope (Olympus Corporation, Japan) was used for microscopy observation and taking photographs.

The fixed ovary was embedded in paraffin, sectioned, stained, and finally stained with aniline blue. LEICA ASP200S tissue dehydrator (LEICA Corporation, Germany) and LEICA EG1150C embedding machine (LEICA Corporation, Germany) were used to perform gradient dehydration and paraffin embedding of the ovary. Fifteen µm thick sections were used LEICA RM2255 microtome (LEICA Corporation, Germany). Finally, the sections were stained with 0.5% aniline blue solution, scanned and captured using Olympus IX73 fluorescence microscope. At least 20 styles and ovaries were observed for each sample.

RNA extraction, library construction and sequencing

The total RNA was obtained using the RNAprep pure Plant Kit (Tiangen Biochemical Technology Corporation, China). The quality and quantity of the total RNA was determined using 1% agarose gel electrophoresis and the Nanodrop 2500 (Thermo Fisher Scientific, USA). Total RNAs from six samples (CKB, CKIB, CKFB, BP, IBP, and FBP) were subsequently used for RNA library construction and sequencing with the BGISEQ-500 platform (Beijing Genomic Institution, China)^[27]. To obtained the clean sequencing reads for subsequent analysis, all the generated raw reads were filtered to remove low-quality, joint contamination or unknown base N content reads.

Genome alignment, correlation analysis and DEGs identification

The clean reads were mapped to the reference genome of tea plant (Camellia sinensis 'Shuchazao') using HISAT2 (v.2.1.0)^[28]. StringTie software (v.1.3.5) was applied to assembly transcripts with default parameters, and the fragments per kilobase of exon per million reads mapped (FPKM) were calculated based on sequencing depth and the length of gene^[29]. The correlation of biological repeats was evaluated by R (v.3.3.2) software with 'corrplot' package (v.0.84). Differential expression analysis was performed using the Cuffdiff software. The transcripts with false discovery rate (FDR) ≤ 0.05 and log₂(FC) > 1 were extracted to be differentially expressed.

Clustering and functional enrichment analysis

DEGs cluster analysis was performed by R with 'pheatmap' and 'mfuzz' packages. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were conducted using TBtools (v.1.0). GO terms and KEGG analysis with P-values ≤ 0.01 were regarded assignificantly enriched by DEGs, visualizing the enrichment result by R with 'ggplot2' package.

Identification of putative SI genes

The Ribonuclease_T2.hmm (PF00445) and S_locus_glycop.hmm (PF00954) were downloaded from the Pfam database (www.pfam.xfam.org/), which were used to identify the putative S-RNase and SRK genes from the tea plant reference genome using HMMER (v.3.3.1) with an E-value 1.0E−5. The SMART website (http://smart.embl-heidelberg.de/) was further used to ensure these proteinscontained conserved domains.

Validation by qRT-PCR

The specific primer pairs of qRT-PCR were designed in non-conserved region with the software Primer 5.0. TB Green Premix Ex TaqII (Takara, China) were used for qRT-PCR on a CFX96 real-time PCR machine (Bio-Rad, USA). The 2^−ΔΔCᴛ and 2^−ΔCᴛ method was used to calculate the relative expression levels based on expression of the CsGAPDH gene, with three biological replicates and three technical replicates^[30].

The mature pollen of 'Shuchazao' tea variety was collected and its viability was tested by pollen germination in vitro. It was shown that the pollen grains have excellent viability and are suitable for pollination (Fig. 1a & b). Tissue staining was performed on the pistils at 12, 24, 36, 48 and 60 HASP. Microscopic observation showed that pollen germinated on the stigma at 12 h, the pollen tubes passed through the stigma and reached the upper part of the style at 24 h, the pollen tubes have extended to the lower part of the style at 36 h, the pollen tubes have passed through the style and entered the ovary cavity at 48 h, and the number of pollen tubes entering the ovary cavity increased at 60 h after self-pollination (Fig. 1c & d). It can be seen that self-pollination at 48 h has met the requirements of the sampling time at different flowering stages.

Figure 1. 

Growth of pollen tube in pistil after self-pollination. (a) 1−2 Picture and internal structure of balloon stage flower used to collect pollen; 3 Pollen. (b) Pollen germination experiment in vitro. (c) 1−5 Pollen tube growth in style after self-pollination. (d) 1−5 Longitudinal section of the ovary after self-pollination. HASP: hours after self-pollination; PT: pollen tube.

Comparison of SI at different flowering stages

The flowers at bud stage, initial bloom stage and full bloom stage were used for self-pollination, and pistils of 48 HASP were collected for microscopic observation. The pollen tubes at all the three different flowering stages can pass through the base of the style and enter the ovary cavity after self-pollination (Fig. 2a). However, the proportion of pollen tubes entering the ovule at bud and full bloom stages were higher than at initial bloom stage (Fig. 2b). Although the number of pollen tubes in styles at bud stage and full bloom stage is relatively small, the proportions of pollen tubes accessing the ovules at bud and full bloom stages are high, accounting for 70.0% and 50.0%, respectively. In comparison, only 10.0% of the samples have pollen tubes in the ovules at initial bloom stage. The results indicate that self-pollination at bud and full bloom stages may increase the self-fertilization rate of tea plants.

Figure 2. 

Growth of pollen tube in different flowering stage at 48 HASP. (a) 1−3 Pictures of flowers at different flowering stages; 4−6 Growth of pollen tubes in style at 48 HASP at different flowering stages; 7−9 Cross section of ovary at 48 HASP. (b) Fluorescence microscope was used to observe and count the PTS and PTO at 48 HASP. HASP: hours after self-pollination; BP: pistil at 48 HASP at bud stage; IBP: pistil at 48 HASP at initial bloom stage; FBP: pistil at 48 HASP at full bloom stage; PT: pollen tube; PTS: proportion of pollen tube through the base of the style; PTO: proportion of pollen tube into ovule.

Identification of DEGs among samples at different flowering stages

To study the molecular mechanism of SI at different flowering stages, six RNA libraries (CKB, CKIB, CKFB, BP, IBP, and FBP) were sequenced. The detailed information including the number of clean reads, the proportion of clean reads aligned to the reference genome, the proportion of clean reads with Q20 quality, and the average GC content was shown in Supplemental Table S1. The correlation analysis of FPKM values of 18 samples showed that the correlation coefficients among the three biological repeats of each treatment were all above 0.84, representing that there was a high correlation among the biological replicates of these samples (Supplemental Fig. S1).

Based on transcriptome analysis, the identified DEGs were divided into seven groups for comparison (CKB vs BP, CKIB vs IBP, CKFB vs FBP, CKB vs CKIB, CKIB vs CKFB, BP vs IBP, and IBP vs FBP) (Fig. 3a). Subsequently, DEGs that may be associated with self-pollination were identified from five groups (CKB vs BP, CKIB vs IBP, CKFB vs FBP, BP vs IBP, and IBP vs FBP). Venn diagram analysis showed that a total of 38 DEGs were identified in all these comparison groups (Fig. 3b). These genes are mainly annotated as peroxidase, serine/threonine protein kinase, ubiquitin protein ligase, synaptic fusion protein, calcium dependent protein kinase and polygalacturonase (Supplemental Table S2). Cluster analysis showed that these DEGs can be divided into five clusters in pistils after pollination at different flowering stages, and the genes in each cluster had similar expression patterns (Fig. 3c). Among them, cluster 2 contained the most genes, cluster 3 had the highest expression at initial bloom stage, cluster 1 and cluster 4 had the highest expression at full bloom stage, and cluster 2 and cluster 5 had the highest expression at bud stage.

Figure 3. 

Analysis of differentially expressed genes (DEGs) after self-pollination. (a) Numbers of DEGs in different comparisons. (b) Venn-flower diagram of DEGs. CKB: unpollinated pistil at bud stage; CKIB: unpollinated pistil at initial bloom stage; CKFB: unpollinated pistil at full bloom stage; BP: pistil at 48 HASP at bud stage; IBP: pistil at 48 HASP at initial bloom stage; FBP: pistil at 48 HASP at full bloom stage. (c) Cluster analysis of all DEGs.

Analysis of DEGs related to SI at different flowering stages

To screen DEGs that may be related to SI between the bud stage and initial bloom stage, we analyzed the two comparison groups (CKB vs CKIB, BP vs IBP) by Venn diagram (Fig. 4a). The results showed that 1,463 DEGs may be related to low SI at bud stage. Cluster analysis showed that these genes were divided into three groups at different flowering stage pistils after pollination (Fig. 4b). The expression of genes in cluster1 (823 genes) was the highest at bud stage pistils after pollination, which may be the key genes that led to low SI at bud stage. The expression of genes in cluster2 (369 genes) was the lowest at bud stage pistils after pollination, and it might be the related genes of SI that had not been expressed. KEGG enrichment analysis showed that 'metabolism', 'carbohydrate metabolism' and 'amino acid metabolism' included the most genes, followed by 'transcription factor', 'environmental adaptation' and 'lipid metabolism'. In addition, 'plant-pathogen interaction', 'cysteine and methionine metabolism', 'lycine, serine and threonine metabolism' and 'cytochrome P450' were also annotated (Fig. 4c).

Figure 4. 

Analysis of SI related genes in BP. (a) Venn diagram displays the number of specific genes in BP. (b) Cluster analysis of DEGs in BP. (c) KEGG pathway enrichment analysis of genes in cluster 1 and cluster 2. The rich factor represents the ratio of the number of DEGs and the number of all the unigenes in the pathway.

Two comparative groups (CKIB vs CKFB, IBP vs FBP) were analyzed by Venn diagram, demonstrating that 1409 genes might be related to the low SI at full bloom stage (Fig. 5a). Cluster analysis showed that these genes could be divided into four patterns (Fig. 5b), and the genes in cluster 1 (274 genes) and cluster 4 (679 genes) might be related to low SI at full bloom stage. KEGG enrichment showed that the pathways containing the most genes were the same as those at bud stage, and also enriched to 'transporter', 'signal transduction', 'MAPK signal pathway-plant' and 'cysteine and methionine metabolism' pathways (Fig. 5c).

Figure 5. 

Analysis of SI related genes in FBP. (a) Venn diagram displays the number of specific genes in FBP. (b) Cluster analysis of DEGs in FBP. (c) KEGG pathway enrichment analysis of genes in cluster 1 and cluster 4.

It is noteworthy that high SI at initial bloom stage was observed, thereby we further analyzed the DEGs of CKIB vs IBP. Cluster analysis showed that 507 DEGs could be divided into five patterns, which cluster 4 contained the most genes, the DEGs in cluster 2 (85 genes) and cluster 5 (41 genes) were probably related to SI (Supplemental Fig. S2a). KEGG enrichment showed that 'signaling and cellular processes', 'transporters', 'membrane trafficking', 'transcription factors', 'cytoskeletal proteins', 'endocytosis', 'plant hormone signal transduction', 'MAPK signaling pathway-plant' and 'ubiquitin system' were mainly enriched (Supplemental Fig. S2b). To further understand the function of DEGs related to the difference of SI at different flowering stages, we obtained annotation information from six databases (GO, KEGG, KOG, Pfam, TrEMBL and NR) (Supplemental Table S3).

Validation of DEGs at different flowering stages

To verify the reliability of DEGs based on transcriptome data, we randomly selected 20 DEGs and verified their expression patterns by qRT-PCR (Fig. 6). The sequences of primers were listed in Supplemental Table S4. Among them, six genes (CsSHMT, CsGCSH CsBCAA2, CsCYP1, CsCYP2and CsMYB4) were related to low SI at bud stage, seven genes (CsMYB24, CsC4Hc, CsMYB102, CsERF1, CsACC1, CsACC2, and CsACC3) were related to low SI at full bloom stage, and the remaining seven genes (CsBST5, CsMGL, CsCDPK1, CsLIMD2c, CsAVT3B, CsGF14, and CsPRK4) were related to high SI at initial bloom stage. As a consequence, the results of qRT-PCR validation of these genes were highly consistent with the expression patterns of transcriptome data. The results indicated that the transcriptome data had high reliability and accuracy, providing many candidate genes that maybe related to SI system in tea plant.

Figure 6. 

qRT-PCR validation of 20 DEGs related toSI. Gray pillars represents the results of qRT-PCR, and the black line represents the results of RNA-seq.

Expression patterns of S-RNase and SRK genes in different flower organs

In various plant species, S-RNase and SRK genes were reported to function in SI system. By alignment, we screened three S-RNase genes (CsRNSL1～CsRNSL3), one known S-RNase gene (CsRNS), and eight SRK candidate genes (CsSRKL1～CsSRKL8). To further identify candidate genes that involved in SI, tissue specific expression patterns were performed in different flower organs. The sequences of the 12 primer pairs are listed in Supplemental Table S5.

As a result, only CsRNS was highly expressed in style among the four S-RNase candidate genes, CsRNSL1 had the highest expression in calyx, CsRNSL2 and CsRNSL3 had the highest expression in leaf (Fig. 7a). Among the eight SRK candidate genes,CsSRKL5 and CsSRKL8 had extremely higher expression level in style than in other organs, while CsSRKL1 was highly expressed in calyx and leaves (Fig. 7b). Interestingly, CsSRKL2, CsSRKL3, CsSRKL4, CsSRKL6 and CsSRKL7 had extremely high expression level in calyx. Therefore, CsRNS, CsSRKL5 and CsSRKL8 are probably critical candidate genes related to SI system in tea plant.

Figure 7. 

Expression patterns of the putative 12 SI genes in different flower organs. (a) The relative expression levels of four S-RNase genes in seven different flower organs. (b) The expression patterns of eight SRK genes in seven different flower organs.

Expression levels of three candidate SI genes at different flowering stages

To further understanding the potential roles of CsRNS, CsSRKL5 and CsSRKL8 in SI system, we analyzed their expression levels in pistils after pollination at different flowering stages. The results showed that the expression levels of CsRNS at initial bloom and full bloom stages were significantly higher than at bud stage (Fig. 8). Similarly, both CsSRKL5 and CsSRKL8 had extremely high expression level at initial bloom and full bloom stages. It is speculated that the low expression levels of CsRNS, CsSRKL5 and CsSRKL8 in pistils may be one of the reasons for the low SI at bud stage.

Figure 8. 

The relative expression levels of putative three SI genes at different flowering stage pistils.

Availability of data and materials

The RNA-seq data have been uploaded to the SRA database of NCBI (accession number: PRJNA985274).

Aniline blue staining solution can be used to stain callose in pollen tube, which shows bluish green under UV excitation^[31]. Herein, we observed the growth state of pollen tube after self-pollination. The results showed that the pollen tubes could pass through the base of style and enter the ovary cavity at 48 h after self-pollination (Fig. 1c), which was basically consistent with previous observations^{[23, 26]}. However, a previous study showed that the pollen tube could enter the ovary cavity at 72 h after self-pollination^[24]. It is probably due to some biological differences between tea cultivars, such as different style length or different pollen germination efficiency, or it may also be caused by geographical location and flowering climate. Previously, it was reported that the pollen tube would grow twisted and callose would accumulate at the base of style after self-pollination, while this phenomenon was not found in our study and a previous study^[22]. Interestingly, we found that the pollen tube in the style at bud stage was significantly less than that at initial bloom and full bloom stage after pollination (Fig. 2a), which may be that the mastoid cells of the stigma were not fully developed, the mucus secreted was less, and the pollen on the stigma was less when pollinated.

In plants with SSI and GSI systems, self-pollination on pistils at bud and full bloom stages has been reported to overcome SI^{[32, 33]}. For instance, bud pollination can significantly reduce their SI in broccoli^[34], self-pollination at bud and full bloom stage could increase the seed setting rate in pear^[35], self-pollination at bud stage could avoid the SI barrier while the SI did not improve significantly with the aging of pistil after flowering in Lycopersicon peruvianum^[12]. In SSI and GSI systems, although there are some differences in the results of bud pollination and full bloom pollination to overcome SI, they all affect the seed setting rate of selfing. In tea plants, we found that there were differences in SI during self-pollination at different flowering stages.

During self-pollination at bud stage in pistils, the incompatible substances had not been formed or the content was lower than the threshold of toxicity, thereby their own pollen could not be recognized, and the pollen tube grew naturally and entered the ovule to complete fertilization. Based on qRT-PCR validation, CsRNS had extremely high expression in style and the expression level of CsRNS at bud stage pistils was significantly lower than initial bloom and full bloom stage, which was consistent with the results that S-RNase in citrus was specifically expressed in style and was lower at bud stage^[2]. CsSRKL5 and CsSRKL8 were highly expressed in the style (Fig. 7b), and the expression level at bud stage was significantly lower than initial and full bloom stages (Fig. 8), which had the same rule as some plants in SSI system^[34]. Studies have shown that S-RNase and SRK homologous genes played an important role in SI system in many plants, such as Citrus grandis^[36], Drosera adelae^[37], Senecio squalidus^[38] and Fragaria viridis ^[39]. Although some S-RNase and SRK homologous genes were not specifically expressed in pistils, they may played an indirect functions in tea plant SI system. Many DEGs of self-pollination at bud stage were annotated into the 'plant pathogen interaction' pathway, which is correlated with the origin hypothesis of SI system^[33], and is consistent with the previous reports in tea plant^[24] and citrus^[40].

During self-pollination at full bloom stage, the pistil began to wilt but no insect pollinated, and the pistil began to accept the pollen from the flower to continue the offspring. Based on transcriptome data, we found that the expression levels of three ACC genes in pistils at full bloom stage were increased significantly (Fig. 6). ACC is usually used to induce ethylene reaction, and ACC signal is involved in promoting the secretion of chemical inducer in pollen tube and promoting pollen tube entering ovule^[41]. In Petunia, pretreatment of stigma with an ACC synthesis inhibitor aminoxyacetic acid (AOA), before pollination self-incompatible lines can promote the growth of pollen tube, and the pollen tube will not appear any characteristics of programed cell death (PCD)^[42]. It is suggested that ACC may play an important role in self-pollination at full bloom stage in tea plants.

During self-pollination at initial bloom stage, the expression level of CsSRK genes in pistil increased significantly (Fig. 8). SRK might recognize some proteins in pollen, resulting in phosphorylation of SRK, ubiquitination of downstream proteins and PCD of pollen tube in ovary^{[1, 43]}. Moreover, the expression level of CsRNS gene in pistils was high, S-RNase enter the pollen tube through endocytosis and plays a cytotoxic role, which degrade mRNA and rRNA in pollen tubes, resulting in PCD of pollen tubes. Ca²⁺ is a secondary messenger in the signaling network, which play a pivotal role in pollen tube elongation and reorientation^[44]. Therefore, the increase of Ca²⁺ will transfer the downstream SI signal, the actin of pollen tube will be depolymerized, and finally the pollen tube will undergo PCD.

Transcription factors (TFs), such as MYB, WRKY, ERF and NF-Y, played critical roles in regulating plant growth and development, secondary metabolism and signal transduction^[45−48]. MYB have diverse roles in regulating leaf and flower development, secondary metabolite biosynthesis, and environmental stress responses in tea plants^[49]. WRKY and ERF usually activate salicylic acid, jasmonic acid and ethylene resistance signaling pathways to further regulate downstream gene expression in response to biological stress^{[45, 46, 50]}. Moreover, WRKY participates in MAPK signaling pathway, which is very important in defense response^[51], and is similar to the identification of pollen and pistil in SI system. NF-Ys are sequence-specific TFs with histone-like subunits that specific bind to CCAAT boxes on promoters and enhancers, and they are played critical roles in embryo development, flowering and stress tolerance^{[48, 52]}. Therefore, these TFs may be involved in pollen pistil interaction, and may also cause SI by regulating cell apoptosis. It can be speculated that complex defense response and signal transduction occurred in the process of self-pollination at initial bloom stage, which may lead to the occurrence of PCD in the ovary cavity of pollen tube and ultimately lead to fertilization failure. Finally, based on the results from our study and previous studies, a proposed molecular mechanism model of SI difference at different flowering stage pistils was provided (Fig. 9).

Figure 9. 

Proposed molecular mechanism model of SI difference at different flowering stage pistils.

In this study, cytological observation and transcriptome sequencing were performed to study the self-pollination of tea plants at different flowering stages. The results showed that the SI intensity at bud stage and full bloom stage were lower than initial bloom stage, and many DEGs related to SI were screened. The low expression levels of CsRNS, CsSRKL5 and CsSRKL8 at bud stage may be related to the low SI, and three CsACC genes may be related to the low SI at full bloom stage. Our results provide a theoretical basis for the study of SI and are of great significance for the breeding efficiency and genetic improvement of SI plants.