Global transcriptome dissection of pollen–pistil interactions induced self-incompatibility in dragon fruit (Selenicereus spp.)

Plant materials

In this work, we used the Hongshuijing and Dahong dragon fruit varieties, which are well-known red flesh varieties with a high market value in mainland China. Hongshuijing is a self-incompatible variety that needs hand-pollination when blooming. However, it produces higher quality fruits. On the other hand, the Dahong variety is self-compatible and does not require hand-pollination. As a result, it is more popular on farms due to easier production. However, its fruit is of relatively low quality. All plants used in this study were grown and preserved in a greenhouse at the Institute of Fruit Tree Research, Guangdong Academy of Agricultural Science, China, ready to be shared with any researcher.

Hand pollination

In the self-pollination experiment, 120 well-developed Hongshuijing flowers were randomly selected and bagged at 8:00 am on a clear day. Then, the bagged flowers were manually shaken at about 8:30 pm to promote self-pollination.

In the cross-pollination experiment, 120 well-developed Hongshuijing flowers were randomly selected at 8:00 am on a clear day, removed the anthers, and bagged. The pollen of Dahong as the pollination variety was collected after blossoming at 8:30 pm, brushed on the emasculated stigma of Hongshuijing for cross-pollination. Then, the cross-pollinated flowers were bagged again.

Fluorescence microscopic observation of pollen growth in pistil

Self and cross-pollinated Hongshuijing flowers at 0, 12, 24, 36, 48, and 72 haps (hours after pollination) were collected. For each treatment, five flowers were sampled for each time point. FAA (50% alcohol: formalin: acetic acid = 85:10:5) was used to glue the pistils. Fluorescence microscopic observation was conducted as follows:

1. The pistils were kept in FAA (50% alcohol: formalin: acetic acid = 85:10:5) at room temperature for 48 h before being stored at 4 °C.

2. The pistils were rinsed with distilled water and immersed in 2 M NaOH for 1 h at 45 °C, then neutralized with glacial acetic acid for 30 s before being rinsed three times with distilled water.

3. The pistils were dyed with 0.1% aniline blue (1 g aniline blue was dissolved with 1000 mL of pre-made K₃PO₄ solution at a concentration of 0.1 M) in the dark for 12 h.

4. The sampled pistils were transversely cut into 1-cm-long segments and stored in order. Then, the segments of pistils were carefully slit longitudinally as pollen tubes clustered in the transmitting tissue. Finally, the transmitting tissue was viewed under a fluorescent microscope (ZEISS Axioplan 2; ZEISS, USA) at a wavelength of 360 nm after being separated and squished.

Effects of glycoproteins on in vitro pollen tube elongation

Nine well-developed flowers of Hongshuijing were randomly selected at 8:00 am on a clear day, and anthers were removed and bagged. The style of each flower was collected after blooming at 8:30 pm and used in protein extraction. Three styles were combined in one sample, and the experiment was conducted in three replicates. Concurrently, the pollen of two varieties, Hongshuijing and Dahong, were collected separately. For each dragon fruit variety, the pollens of three flowers were mixed as a sample, and the experiment was conducted in three replicates. Total soluble proteins were extracted from the styles of Hongshuijing as described in the protocol of the plant total protein extraction kit (Sangon, Shanghai, China). Glycoproteins were separated from the total soluble proteins using the FOCUS™ Glycoprotein kit (Sangon, Shanghai, China), following the manufacturer’s instructions. The soluble proteins without glycoproteins collected from the style were named glycoproteins-free protein fractions.

The pollen of Hongshuijing (self-incompatible) and Dahong (cross-compatible with Hongshuijing) were cultured in a liquid medium (20% sucrose + 300 mg/L H₃BO₃ + 100 mg/L Ca(NO₃)₂.4H₂O + 300 mg/L KNO₃ + 300 mg/L MgSO₄.7H₂O under 28 °C, pH 7.0) supplemented with total soluble proteins sample, isolated glycoproteins, or the total proteins fraction from which glycoproteins were removed (glycoproteins-free fraction), respectively. The liquid medium without any extra protein added was used as control. The protein concentration in either of the protein samples described above was10 ng L⁻¹. For in vitro germination experiment, 0.1 g of pollen grain was placed in a five mL liquid medium and incubated at a temperature incubator (28 °C and RH80%) for 12 h, and the pollen tube growth was investigated under a microscope. For each sample, three fields were observed, including no less than 30 pollen grains. The experiment was conducted in three replicates.

Sample preparation and sequencing

A total of 21 pistil samples from seven treatments were used in this experiment. The 21 samples consisted of the 12, 24, and 36 h after self-pollination samples, the 12, 24, and 36 h after cross-pollination samples, and the control sample (without pollination). Five pistils from distinct plants were combined in one sample for each treatment. The experiment was conducted in three replicates. All these samples were subjected to RNA-seq and qPCR analysis. Total RNA was extracted with CTAB plus using an OMEGA Plant RNA Isolation kit according to the protocol described by Ouyang et al. (2014). The library construction and sequencing were conducted using a standard protocol by Novogene (Beijing, China).

DEG identification and functional annotation

The raw data were filtered to remove adapter or low-quality sequences and assembled by the Trinity pipeline (Haas et al., 2013). The transcripts per million (TPM) approach was used to calculate unigene expression (Wagner, Kin & Lynch, 2012). The differential expression gene (DEG) analysis was performed using the R package “edgeR” 3.36.0 with the following thresholds: q-value 0.05 and |log2(foldchange)| > 1.5 (Robinson, McCarthy & Smyth, 2010). q-value was adjusted p-value using Benjamini and Hochberg’s approach (Benjamini & Hochberg, 1995). Power analysis was conducted with the R package RNASeqPower (Therneau & Stephen, 2014). The following databases were used to determine functional annotations: NCBI non-redundant (Nr) protein sequences, NCBI non-redundant nucleotide (Nt) sequences, protein family (Pfam), Clusters of Orthologous Groups of proteins (KOG/COG), Swiss-Prot, a manually annotated and reviewed protein sequence database, KEGG Ortholog (KO) database, and Gene Ontology (GO).

Co-expression analysis

A weighted gene co-expression network analysis (WGCNA) was used in this study (Langfelder & Horvath, 2008). The input expression matrix included all expression data of DEGs from the pistil RNA-Seq. This R program divides the DEGs into multiple modules (gene sets) based on direct and indirect gene correlations. The systematic strategy, which included GO and KEGG enrichment analyses, was used to investigate the functions of these gene sets. The significance cutoff of GO and KEGG enrichment analysis is q-value (false discovery rate adjusted p-value) <0.05. The correlation matrix generated by the WGCNA was used to build the co-expression network. The network consisted of the gene connections and their 50 most correlated genes. Cytoscape v3.8.0 was used to visualize the network, and the nodes were ordered using the yFiles Layout Algorithms (Tuikkala et al., 2012).

Real-time Quantitative PCR (qPCR) analysis

In total, 20 differentially expressed genes identified in RNA-seq datasets were randomly selected, and their expression was confirmed by qPCR. The qPCR was conducted using the SYBR Premix ExTaqII Kit (Takara) with the ABI StepOne Real-Time PCR system (Applied Biosystems, San Mateo, CA, USA). The thermal cycling conditions were as follows: an initial denaturation of 95 °C for 30 s, followed by 40 cycles of denaturation, annealing, and elongation (95 °C for 5 s and 60 °C for 30 s), then a melt-curve step (0.3 °C/s, from 60 to 95 °C). The relative expression level of each gene was calculated via the 2^−ΔΔCt method (Livak & Schmittgen, 2001) with three biological replicates. In addition, actin and EF1- α were used as reference genes (Chen et al., 2019).

Pollen growth in pistil in vivo and effects of style glycoproteins on pollen tube elongation in vitro

To observe in vivo pollen tube elongation in the pistil of Hongshuijing, we conducted self-pollination and cross-pollination of Hongshuijing with Dahong pollen. The pollen grain can germinate on the stigma of self- and cross-pollination (Fig. 1A, Table S1). After self-pollination, the pollen tube elongated to 17.5% of style after 12 h, reached the middle parts (48.5%) of style after 24 h, and then 50.2% and 50.4% of style after 36 and 48 h, respectively. It demonstrated that the pollen tube elongated to the middle parts of the style and almost ceased extension at 24 h after self-pollination. However, in cross-pollination, the pollen tube elongated to 19.9% of style after 12 h, extended to nearly 2/3 (69%) of the style length after 24 h, and then 99.3% and 100% of style after 36 and 48 h, respectively.

It is indicated that the elongation of the self-pollinated pollen tube in style was significantly slower than that of the cross-pollinated pollen tube at different time intervals (Fig. 1B). The pollen tube almost ceased elongation at 24 h after self-pollination, and the pollen tube in cross-pollination entered the ovary after 36 h. This result displays that dragon fruit exhibits gametophytic self-incompatibility (GSI).

Based on the hypothesis that S glycoprotein is the leading participant in GSI, we performed an in vitro pollen germination test using a germination medium supplemented with total soluble proteins, glycoproteins’ fraction isolated from total soluble proteins, and a glycoproteins-free fraction of total soluble proteins derived from the style of Hongshuijing. We discovered that the pollen tube elongation of Hongshuijing (self-incompatible) was sensitive to the total soluble proteins and glycoproteins but tolerated the non-glycoproteins compared to the control (Fig. 1C). However, the pollen tube elongation of Dahong (cross-compatible with Hongshuijing) was not affected by the three protein treatments (Fig. 1D).

The results of in vivo pollen tube elongation and the effect of the style glycoproteins on in vitro pollen tube growth demonstrated that dragon fruit exhibits the GSI, and glycoproteins are the players in the mechanism of self-incompatibility. This is consistent with the findings that the female S determinants are glycoproteins in the two known GSI mechanisms, either the S-RNase-based type or the poppy type of GSI.

Transcriptomes analysis of self- or cross-pollinated pistils

We sequenced RNA from pistil tissue from 21 different samples to investigate and compare the global transcriptome dynamics in Hongshuijing pistils in response to self and cross-pollination. These samples included control (without pollination), SP12, SP24, and SP36, CP12, CP24, and CP36 (corresponding to 12, 24, and 36 h after self or cross-pollination, respectively). The transcriptome of Hongshuijing pollen was investigated individually. A total of 155.83 GB of clean data was retrieved from 21 samples, with each sample containing at least 6.51 GB of clean data with a Q30 quality score ≥ 95.21% (Table S2).

The vast majority of the treatments were found to be highly reproducible. Only one CP36 sample replicate did not cluster with the other two samples (Fig. S1). In total, 43,954 unigenes were evaluated with an N50 of 1586 bp, and the guanine-cytosine (GC) content was 49.20%. The unigene length ranged from 201 to 5869 bp with an average length of 1,430 bp, of which 67.12% (29,503) were longer than 1,000 bp. The 42,194 unigenes out of the total 43,954 had functional annotations in the known database (Table 1), 91.8% of the annotated unigenes with E-value smaller than 10⁻³⁰ (Fig. S2), and the final reference transcriptome is 61.2 Mb (unigene of dragon fruit pistil, average sequencing depth =108). Our previously uploaded pollen transcriptome data on NCBI (PRJNA822018) was also used for male S factors identification in this study. To validate the accuracy of our RNA-seq data, qPCR was used to verify the randomly selected 20 DEGs (the sequence and corresponding primers are listed in Table S3). Our results confirmed that the expression patterns of these genes were very similar to those found in RNA-seq (Fig. 2), confirming that our RNA-seq expression data were highly reliable and can be used for subsequent experimental analysis.

Identification of DEGs and cluster analysis

The power analysis revealed that fold change (FC) >=1.5 could be consistently detected with a power of 0.9 and a false discovery rate (FDR) of 0.05 with a sample size of 7 when accounting for average sequencing depth and biological variance (Fig. S3). A total of 6,353 significant differentially expressed genes (DEGs), including 58 transcription factors (TFs), were identified from 15 pairwise comparisons. These included six comparisons of the self-pollinated samples (Fig. 3A, control vs. SP12, control vs. SP24, control vs. SP36, SP12 vs. SP24, SP12 vs. SP36, and SP24 vs. SP36), six comparisons of the cross-pollinated samples (Fig. 3B, control vs. CP12, control vs. CP24, control vs. CP36, CP12 vs. CP24, CP12 vs. CP36, and CP24 vs. CP36), and three comparisons of the self and cross-pollinated samples (Fig. 3C, SP12 vs. CP12, SP24 vs. CP24, and SP36 vs. CP36). In the case of self-pollinated pairwise comparisons, control vs. SP24 had the most DEGs (Fig. 3A). However, control vs. CP36 yielded the highest number of DEGs in CP comparisons (Fig. 3B). Furthermore, SP24 vs. CP24 exhibited more DEGs than the other two comparisons (Fig. 3C). These findings correlated with the result of the SP experiment which depicted that the pollen tube ceased elongating after around 24 h (Fig. 1A and 1B), indicating that more genes have been up or down-regulated to suppress the pollen tube elongation.

SI is the consequence of a complex interplay of biological processes, most of which are unknown in dragon fruit. To better understand these processes, we artificially divided them into two subprocesses: pollen recognition and rejection. We discovered that the pollen recognition occurred before 12 haps because the pollen tube length in SP treatment was already shorter than in CP treatment (Fig. 1A and 1B). In pairwise comparisons, 982 DEGs were common between control vs. SP12 and control vs. CP12 (Fig. 3D). We hypothesized that these pollen recognition-related genes would be included in the 982 DEGs because recognition occurs in both SP and CP treatment. Furthermore, the expression profiles of these 982 DEGs were similar, with expression values increasing or decreasing in both SP and CP samples from 0 to 12 hap. Following recognition, the stigma can assess whether the pollen tube is from self-pollen and determine whether to reject or accept it. Pollen tube rejection is a continual interaction between pollen and pistil until the pollen tube stops elongation (Chen et al., 2018). It is unclear when exactly the rejection occurs, but it is certain that it will only happen in SP pistils. Therefore, most of the rejection-related genes will be differentially expressed in SP pistils and thus will be included in the 214 DEGs that existed between SP and CP (Fig. 3C). Therefore, the two important DEGs sets, the 982 DEGs as candidate pollen recognition-related genes and the 214 DEGs as potentially rejection-related genes were identified.

Using WGCNA analysis with the soft threshold of 12 (Fig. S4 A, the approximate scale-free fit index can be attained when the soft-thresholding power β was set as 12) all 6,353 DEGs between the pistils at different intervals of time after self- or cross-pollination were assigned to ten modules. These modules are depicted in different colors: black, blue, brown, green, magenta, pink, red, turquoise, yellow, and grey. The grey module, which consisted of 88 genes, represented the group of genes that could not be classified. (Table S4; Fig. S4B). The biggest module depicted in turquoise contained 2483 genes, and the smallest module depicted in magenta contained 74 genes. Of the two crucial DEGs sets, in the 214 DEGs gene set related to pollen recognition, 132 (61.7%) of the 214 genes were assigned to the green module, which contains a total of 461 genes, and in the 982 DEGs gene set related to rejection, 534 (54.4%) of the 982 genes were assigned to the turquoise module, which contains the total of 2,483 genes (Table S5). A module is a collection of genes that work together to form a single or a series of interconnected processes. We hypothesized that the green and turquoise modules are the other two DEG sets containing SI-related genes. Therefore, four DEGs sets, the pollen rejection related 214 DEGs (Fig. 3C, the DEGs between CP12 and SP12, CP24 and SP24, CP36 and SP36), the recognition related 982 DEGs (Fig. 3D, the shared DEGs between control and CP12, control and SP12), the green module and the turquoise module (two of the modules generated by WGCNA with a large proportion of overlap with the 214 and 982 DEGs sets), were identified through two different approaches. Furthermore, the identified gene set overlap, the 132 DEGs overlapped between the 214 DEGs and the green module, and the 534 DEGs overlapped between the 982 DEGs and the turquoise module were narrower gene sets potentially involved in SI.

GO and KEGG functional enrichment analyses of the key DEGs sets

To study what could be the functions of the identified gene sets, we conducted a general survey with GO and KEGG enrichment. The results were particularly pronounced for the two narrowed-down gene sets, the overlapped 132 DEGs and 534 DEGs. For the overlapped 132 DEGs, there were 51 GO terms, including 24 for biological process, 23 for molecular function, and four for cellular components significantly enriched. Most of the genes in the biological process are involved in the carbohydrate metabolic process. The molecular function terms are mainly related to hydrolase activity and catalytic activity. The cellular component genes were located in the cell wall, external encapsulating structure, cell periphery, and microtubule (Fig. 4). KEGG pathway and enrichment analysis demonstrated four significantly enriched pathways, mainly in the pathways of pentose and glucuronate interconversions (Fig. 4).

For the overlapped 534 DEGs, there were 216 GO terms, including 156 for biological process, 55 for molecular function, and 7 for cellular components significantly enriched. The enriched GO terms mainly included small molecule metabolic processes, carboxylic acid metabolic process, oxoacid metabolic process, transferase and oxidoreductase activity, extracellular region, apoplast, and cell wall (Fig. 4). The KEGG pathway and enrichment analysis revealed that nine pathways were significantly enriched, primarily in the pathways of secondary metabolite and amino acid biosynthesis, cysteine, and methionine metabolism (Fig. 4). GO and KEGG functional enrichment analyses of the other 4 DEGs sets, 214, 982, the green module, and the turquoise module, are depicted in Figs. S5 and S6.

Identification of S genes and non-S locus genes

Searching for the homologous genes in other species is the standard method to identify the candidate genes. Based on the most common S-RNase-based GSI mechanisms, six candidates of S-RNase were identified in the Hongshuijing pistil database, and 158 F-box genes were identified in the Hongshuijing pollen database (Table S6). However, we did not find any homologous genes of PrsS and PrpS, so we exclude the possibility of P. rhoeas type SI in dragon fruit. As a result, we hypothesized that the identified S-RNase and F-box genes are candidate female and male S genes, respectively.

In addition to well-known S genes, non-S locus genes have also been discovered to participate in GSI (Zhang et al., 2017; Su et al., 2020; Ma et al., 2021).To simplify the transcriptome gene expression data, we used the principle of data dimensionality reduction. A gene co-expression analysis was carried out based on the TOM matrix generated by WGCNA analysis. A gene was designated as a node in each network in our analyses, and its co-expression with other genes was represented as lines (edges). The degree of connectedness of a node was defined as the number of genes associated with it. Hub genes, which act as representations of a gene set, had the highest degree of interconnectedness. Here, the top 50 genes of the co-expression network were defined as hub genes (Tables S7 and S8). By querying the Pearson’s Correlation Table (Degree of Freedom = 20, Pearson’s correlation >0.536800), we found that all correlations in the network reached a highly significant level (p < 0.01).

The network hub is defined by a mathematical algorithm rather than a biological method. However, the biological significance can be speculated based on functional gene annotations. Among the first 10 hub genes of co-expression network for the 534 DEGs, seven of them may directly or indirectly be involved in SI response (Table S7), which were annotated as F-box (c37611_f1p1), U-box E3 ubiquitin ligase (c104651_f1p2) (Kim et al., 2021), pistil-specific extension-like proteins III (PELPIII) (c28784_f2p0), 120 kDa (c96768_f1p3) (Schultz et al., 2002; Hegedűs, Lénárt & Halász, 2012), serine/threonine-protein kinase isoform X2 (c62321_f1p0) (Nie et al., 2019), cysteine-rich receptor-like protein kinase (c106266_f2p1) and aquaporin (c27095_f28p9) (Ikeda et al., 1997). The other three hub genes, YABBY4 (c30021_f1p0), which is involved in seed abortion in Vitis vinifera L. and elongates the pistil abnormally in tomatoes when overexpressed (Massonnet et al., 2020), xyloglucan galactosyltransferase XLT2 (c131771_f1p0), which is involved in plant cell wall polysaccharide xyloglucan biosynthesis (Schultink, 2013), and ferredoxin-like protein (c1725_f3p0) which functions in stress resistance (Huang et al., 2020), were among the 10 hub genes that had no obvious connection to known function in SI response.

For 132 DEGs, six of the first 10 hub genes may directly or indirectly be involved in SI response (Table S8). Those were annotated as beta-galactosidase (c161600_f1p9) (Wanyi et al., 2019), pectinesterase (c127787_f1p9, c70692_f1p1, c99737_f1p2) (Wanyi et al., 2019), 120 kDa (c28368_f3p4) (Lin et al., 2014) and arabinogalactan protein (c9808_f1p3) (Pereira, Lopes & Coimbra, 2016). Interestingly, the other four hub genes, c1947_f8p6, c7984_f1p0, c103690_f1p1, and c205802_f1p2, were not annotated in any databases (Table S9). However, we cloned them from dragon fruit pistils and fully verified their sequences. We hypothesized that they are new genes from dragon fruit and may participate in the GSI regulation.

In the present study, a total of 13 SI-related genes were identified as being non-S locus genes based on the results of gene co-expression analysis, which indirectly supports the validity of this analysis. Based on this, there is a definite possibility that the other hub genes (Tables S7 and S8) from two co-expression networks may also function in some unknown way in the SI mechanism. We also discovered 58 differentially expressed TFs among the DEGs in this investigation. Through the building of a co-expression network for the 58 TFs, we observed that these TFs clustered into three groups, with the first group having only one member (TF1) and the other two groups having 22 (TF22) and 35 (TF35) members, respectively (Fig. S7). Furthermore, the top 50 hub genes of the 534 DEGs network were fully included in the TF22 sub-network, which had a strong correlation with 13 of the 22 TFs (Fig. 5). However, we must point out that 2 of the 13 TFs are also the top 50 hub genes of the 534 DEGs. The workflow diagram of our candidate gene set screening is illustrated in Fig. 6.

The different expressions of the 13 TFs at different time points are depicted in Fig. S8. These TFs were annotated as YABBY4 (c30021_f1p0), ANL2 (c105690_f1p0), OBF1 (c28526_f6p1, c99327_f1p2), ERF43 (c64391_f1p0), ARF2 (c105910_f1p1), BLH 7 (c129877_f1p4), KNAT6 (c42568_f1p1), PIF3 (c166445_f1p1), HY5 (c36498_f1p1, c71338_f1p1) and LHY/CCA1 (c125807_f4p5, c137898_f1p7).

They could also potentially be the non-S locus genes of GSI in dragon fruit. This particularly applies to ERF43 (Su et al., 2020), ARF2 (Tantikanjana & Nasrallah, 2012), BLH7, and KNAT1 (Arnaud & Pautot, 2014), which have already been found to mediate SI responses in Brassica rapa, Ipomoea cairica, and Theobroma cacao.