Two days after flowering, the plants were self- and cross-pollinated. D. chrysanthum pistils were collected 12, 24, 48, 72, and 96 h after self- and cross-pollination (HASP and HACP, respectively), and D. longicornu pistils were collected 2, 3, 4, 5, 7, and 9 days after self- and cross-pollination (DASP and DACP, respectively). The pistils were immersed in a fixing solution (formalin-acetic acid-80% alcohol [1:1:8]) for at least 24 h, rinsed in 70% alcohol, softened in a strong (8 N) sodium hydroxide solution for 3 h, and cleared in distilled water. The pistils were stained with 0.1% water-soluble aniline blue dye dissolved in 0.1 N K₃PO₄ for 12 h and then observed under a fluorescence microscope (Axioskop 40, Zeiss, Germany). The ultraviolet light used (wavelength: approximately 356 μm) facilitated the examination of the growth of pollen tubes on the style because the pollen tubes were lined with callose, which fluoresced bright yellow-green and contrasted strongly with the bluish or grayish fluorescence of the style and ovary tissues. The lengths of the majority of the pollen tubes in compatible and incompatible styles at specific times after pollination were measured using fluorescence microscopy software. Statistical analyses were performed using the software package GraphPad Prism 6 for Mac OS X (version 6.0c; GraphPad Software, Inc., La Jolla, CA, United States). All values were reported as the mean and SEM. The data were presented as the average of three independent measurements with error bars (SEM) indicated. The data were analysed by two-way repeated measures analysis of variance (ANOVA) using Sidak’s post hoc test (^∗∗∗P < 0.001, ^∗∗∗∗P < 0.0001).

Mature pollen, styles (containing stigma) and leaves of D. chrysanthum and the 3 DACP styles (containing pollen), 3 DASP styles (containing pollen) and leaves of D. longicornu were simultaneously collected during the pollination of the same plants, placed in 2-ml tubes in liquid nitrogen, and stored at –80°C until RNA extraction. Three biological replicates of each D. chrysanthum tissue were collected.

Total RNA was extracted from Dendrobium spp. tissues using an RNA prep Pure Plant Kit, and genomic DNA contamination was removed using RNase-free DNase I (both from Tiangen, Beijing, China). The integrity of the RNA was evaluated on a 1.0% agarose gel stained with ethidium bromide (EB), and its quality and quantity were assessed using a NanoPhotometer^® spectrophotometer (IMPLEN, Westlake, CA, United States) and an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, United States). Because the RNA integrity number (RIN) was greater than 7.0 for all samples, the samples were used for cDNA library construction and Illumina^® sequencing, which was completed by Beijing Novogene Bioinformatics Technology Co., Ltd. (China). The cDNA library was constructed using the NEBNext^® Ultra^TM RNA Library Prep Kit for Illumina^® (NEB, United States) with 3 μg of RNA per sample, following the manufacturer’s recommendations. The polymerase chain reaction (PCR) products obtained were purified (AMPure XP system), and the library quality was assessed on the Agilent Bioanalyzer 2100 system. Library preparations were sequenced on an Illumina^® HiSeq 2000 platform, generating 100-bp paired-end reads. Raw sequence reads were deposited in the National Center for Biotechnology Information (NCBI) Sequence Read Archive¹ (SRA) under accession number SRP097204.

Before assembly, adaptor sequences were removed from the raw reads, and FASTQC² reports were generated. Based on this information, the reads were trimmed at both ends. Nucleotide positions with a QC score lower than 20 were masked (replaced by an N), and the resulting high-quality reads were de novo assembled and annotated with TRINITY (Haas et al., 2013). The commands and parameters used for running TRINITY were as follows: Trinity –seqType fq –JM 200G –left sample_1.fq –right sample_2.fq –normalize_by_read_set –CPU 32 –output –min_kmer_cov 2 –full_cleanup. The transcript abundance level was normalized using the fragments per kilobase per million reads (FPKM) method, and FPKM values were computed as proposed by Mortazavi et al. (2008). In genes with more than one transcript, the longest was used to calculate the transcript abundance and coverage.

We first collected S-RNase-related RNase-T2 protein sequences from previous studies (Aguiar et al., 2015; Kubo et al., 2015). Then, we combined the genomic data of P. equestris (SC) (Cai et al., 2015), D. catenatum (partial SI) (Zhang et al., 2016), V. shenzhenica (SC) (Liu et al., unpublished), and A. shenzhenica (SC) (Liu et al., unpublished) with six transcriptome sequences in the proteome datasets.

These S-RNase-related RNase-T2 protein sequences were used as seeds to search against proteome datasets using basic local alignment search tools (tBLASTn and BLASTp) with default parameters. To confirm the identities of orchid S-RNase-related RNase-T2 candidate genes, the hidden Markov model (HMM)-based HMMER program 3.1b2 (Eddy, 2011) was used to identify all proteins containing a Ribonuclease_T2 domain. This domain (PF00445.16 in the Pfam database) (Finn et al., 2016) was then used to perform local searches in the proteome datasets. Sequences obtained from both methods were then aligned and manually adjusted in Multiple Alignment using Fast Fourier Transform (MAFFT) (Katoh and Standley, 2013) using the E-INS-I alignment strategy for sequence integrity analysis. Sequences with obvious errors were excluded from subsequent analyses. Each predicted sequence was subsequently verified using BLAST searches against public databases, including NCBI, Pfam (Finn et al., 2016) and Simple Modular Architecture Research Tool (SMART) (Letunic et al., 2012), to confirm its reliability. The HMM profiles of the Skp1 gene family (PF01466.17 in the Pfam database) obtained in the analysis were used in local searches. To maximize confidence, putative Skp1 sequences were also aligned and manually adjusted in MAFFT using the E-INS-I strategy for sequence integrity analysis. The IPs of all peptides were calculated in ExPASy (Gasteiger, 2003).

As monocots, most orchids are SC species; therefore, the SSK1 genes may be non-functional or not involved in the SI pathway. Thus, we chose orchid species of several SI degrees and phylogenetic positions and allowed for some variability regarding the “WAFE” and “GVDED” motifs when retrieving sequences. Angiosperms have two types (I and II) of SKP1-like genes, with type II being much longer than type I and encoding chimeric proteins (Kong et al., 2007). Multiple Skp1 homologs from the same species were shown to have evolved at highly heterogeneous rates, indicating that they have different evolutionary histories (Kong et al., 2004, 2007). Due to the effect of long-branch attraction, when both gene types are included, the analysis typically gives unstable results; in addition, SSK1-like genes belong to type I SKP1-like genes (Kong et al., 2004, 2007; Huang et al., 2006; Zhao et al., 2010; Matsumoto et al., 2012; Xu et al., 2013). For these two reasons, type II genes were excluded from the analysis.

The amino acid sequences of SSK1 genes, SSK1-like genes from selected land plant species (Kong et al., 2004, 2007; Huang et al., 2006; Zhao et al., 2010; Matsumoto et al., 2012; Xu et al., 2013), and Orchidaceae SSK1-like genes were used for phylogenetic analysis. Multiple sequence alignment was carried out using MAFFT with the E-INS-I strategy and adjusted manually as necessary, and the phylogenetic tree was generated by the maximum likelihood method using PhyML 3.0 (Guindon et al., 2010). The approximate likelihood-ratio test (aLRT) branch support, which was based on a Shimodaira-Hasegawa-like procedure, was estimated with a Whelan and Goldman (WAG) model.

The amino acid sequences of Orchidaceae S-RNase like genes and S-RNase-related RNase-T2s from species with S-RNase-based SI, maize and barley (Igic and Kohn, 2001; Hillwig et al., 2010; Kubo et al., 2015), were used for phylogenetic analysis the same method as that described for SSK1-related genes.

The data sets (the alignments file and the tree file of SKP1 and S-RNase protein sequences) are available through TreeBASE³, and also in the Supplementary Files S3, S4.

The development of pollen tubes that were stained with water-soluble aniline blue in the pistils of D. longicornu and D. chrysanthum was observed 12 h after pollination (HAP) using fluorescence microscopy. In D. longicornu, pollen tubes started to develop in the stigma 2 DASP or DACP (Figures 1, 7). At 3 DAP, when the pollen tubes started to grow into the style, the cross-pollination pollen tubes were longer than the self-pollination pollen tubes; however, approximately 4 DAP, the self-pollination pollen tubes stopped growing at the top of the style. There was a highly significant difference between the self- and cross-pollination pollen tubes beginning at 5 DASP or DACP (Figure 7). Approximately 7 DACP, the pollen tubes developed rapidly and reached the ovary. In D. chrysanthum, the pollen tubes started to develop in the stigma 24 h after self-pollination or cross-pollination (Figures 2, 8). Although the self-pollination pollen tubes grew slowly into the style during the 72 h after pollination and stopped growing before flower (containing ovary) senescence (120 h after pollination), the cross-pollination pollen tubes developed rapidly and reached the ovary 72 h after pollination. There was a highly significant difference between the self- and cross-pollination pollen tubes beginning at 48 HASP or HACP (Figure 8).

Overall, the germination of pollen tubes inside the stigma of D. longicornu occurred later than that in D. chrysanthum, and their growth took longer in D. longicornu styles than in D. chrysanthum styles. The self-pollination pollen tubes of D. longicornu grew to the top of the style before flower senescence, while in D. chrysanthum, they stopped growing before flower senescence. Thus, SI was identified in two orchid species by observing the development of pollen tubes using fluorescence microscopy.

We retrieved 42 sequences from orchids, including five A. shenzhenica (SC), two V. shenzhenica (SC), six D. catenatum (partial SI), and five P. equestris (SC) genome sequences; four sequences from the leaves and self- and cross-pollinated styles (containing pollen) of D. longicornu (SI); and four sequences from the leaves, five sequences from the pollen, and three sequences from the styles (containing stigma) of D. chrysanthum (SI). Two genes had a “WAFGE” motif rather than the conventional “WAFE” motif in the 3’ region; that is, a glycine was inserted between phenylalanine and glutamic acid. One of these genes was found in the leaves, pollen, and styles of D. chrysanthum (Dch1_L, Dch1_P, and Dch1_S, respectively) (Supplementary File S1), and the other was expressed in the styles (containing stigma and pollen) of D. longicornu (Dlo1_O and Dlo1_S) (Supplementary File S1). Similarly, the conserved C-terminus with the “WAFE” motif was replaced by “WAFAE” in one D. catenatum sequence (Dca004172) (Supplementary File S1), and a 16-amino acid deletion was found near the C-terminus of this sequence. Another D. catenatum sequence (Dca007256) presented a “WAFDLICL” motif, and one A. shenzhenica sequence (Ash003413) presented a “WAFEPQQ” motif at the C-terminus (Supplementary File S1). The Dca007256 gene was expressed in stigma, and the Ash003413 gene was also expressed in leaves, stems and tubers (data not shown). The FPKM analysis of the combined transcriptome data of the three D. chrysanthum tissues revealed that all of the abovementioned genes were expressed in D. chrysanthum leaves and, therefore, were not pollen-specific, except Dch3_P (Figure 3, Supplementary File S2, and Table S1).

In the phylogenetic analysis, SSK1 genes, SSK1-like genes from selected land plant species, and Orchidaceae SSK1-like genes were combined. The SSK1-like genes from Orchidaceae were divided into group I (35 genes) and group II (6 genes), and one gene was clustered with one A. thaliana gene (ASK2) (Figure 5) shown in blue. Within group I, all genes expressed in the style and/or pollen were also expressed in the leaf and thus were not pollen-specific genes (Supplementary File S2 and Table S1); this was the case for Dlo2_L, Dlo2_O, and Dlo2_S or Dch2_L, Dch2_P, and Dch2_S. In group II, which was a sister to that formed by the SSK1 genes, no genes were expressed in leaf tissues, and with a few exceptions, all genes in this group were expressed in the pollen of self- and cross-pollinated orchids, suggesting that these might be pollen-specific genes; the exceptions were PEQU_01554 and Ash012058, which were expressed in the leaf, stem, root, and pollen, and Dca025410, which was also found in the root tissue (data not shown).

The phylogenetic tree positions and tissue-specific expression analyses revealed that the genes Dch3_P in D. chrysanthum and Dlo3_S and Dlo3_O in D. longicornu may be Orchidaceae SSK1-like genes, despite their lack of a “GVDED” motif. Indeed, in SI Orchidaceae, the conserved “WAFE/D” motif might have evolved into a sequence different from that generally observed in S-RNase-based GSI.

As RNase-based GSI might be present in Orchidaceae, we attempted to identify an S-RNase gene in this family. Based on four criteria, S-RNase-like genes were identified in A. shenzhenica, V. shenzhenica, D. catenatum, and P. equestris genomes and the transcriptomes from the leaves, pollen, and styles of D. chrysanthum and the leaves and self- and cross-pollination styles of D. longicornu. The criteria were as follows: (1) The genes were similar, at the amino acid level, to those involved in S-RNase-based GSI in Rosaceae, Solanaceae, and Plantaginaceae. (2) The genes encoded a protein in which amino acid pattern 4 was absent; this pattern is found in proteins encoded by non-S-RNase lineage genes only (Vieira et al., 2008a; Nowak et al., 2011). (3) The genes encoded a protein with an IP higher than 7.5 because S-RNases are always basic proteins (Igic and Kohn, 2001; Roalson, 2003). Finally, (4) the genes should be mainly expressed in the style. Given that little is known about the physiology and genetic control of SI in Orchidaceae, we chose SC species (A. shenzhenica, V. shenzhenica, and P. equestris), a partial SI species (D. catenatum), and SI species (D. chrysanthum and D. longicornu) to identify putative S-genes, although mutations that disrupt the coding region might exist in the putative S-locus region (Tao et al., 2007).

We identified 34 RNase-T2-like genes in Orchidaceae: five in A. shenzhenica (SC); four in V. shenzhenica (SC); three in D. catenatum (partial SI); five in P. equestris (SC); two in the leaves, three in the self-pollination styles, and three in the cross-pollination styles (containing pollen) of D. longicornu (SI); and three in the leaves, three in the pollen, and three in the styles (containing stigma) of D. chrysanthum (SI). The features of all 34 gene sequences, including the IP, sequence motifs, intron numbers, and gene location, are summarized in Table 1. The gene intron numbers in D. longicornu and D. chrysanthum were determined by an alignment with D. catenatum and P. equestris genes. Six genes had IPs higher than 7.5: Ash010024, Ash010025, Dlo9_L, Dlo9_S, Dch6_S, and VN_GLEAN_10017549. Twenty-five genes presented amino acid pattern 4 ([R]) and, thus, had S-like amino acid patterns rather than S-RNase amino acid patterns. The patterns “CGS” and “CSS” were found in eight and one gene sequence, respectively. Among the eight “CGS” genes, Dlo9_L, Dlo9_O, and Dlo9_S were the same gene expressed in different D. longicornu tissues, suggesting that this gene was not involved in RNase-based GSI because no style-specific expression was observed. The genes Dch8_L, Dch8_P, and Dch8_S in D. chrysanthum were also the same gene, demonstrating that this gene is also not involved in RNase-based GSI processing. The other two genes (PEQU_20427 and VN_GLEAN_10005747) were found in the SC species P. equestris and V. shenzhenica; we investigated the functions of the genes located at the two gene scaffolds and analyzed the expression levels of the genes with SLF function (data not shown). The results indicated that these genes were not at the putative S-locus region and not primarily expressed in the style, which also suggested that they are not involved in RNase-based GSI processing. The gene Ash013339, which contained the sequence “CSS” in the SC species A. shenzhenica, was also determined to not be involved in RNase-based GSI based on the results of adjacent gene function and expression analyses (data not shown).

Phylogenetic analyses of these 34 genes and S-RNase-related RNase-T2s from species with S-RNase-based SI, maize and barley, were performed (Figure 6). The classification of RNase-T2s into Class-I RNase-T2s, Class-II RNase-T2s, and the S-RNase clade followed that of previous studies (Igic and Kohn, 2001; Hillwig et al., 2010; Kubo et al., 2015). The 34 genes were clustered into Class-I RNase-T2s or Class-II RNase-T2s and were therefore not involved in RNase-based GSI. Three orchid clusters were found in Class-I RNase-T2s: Class-I orchid RNase-T2s I; Class-I orchid RNase-T2s II; and Class-I orchid RNase-T2s III. The genes in Class-I orchid RNase-T2s I clustered with Zea mays kin1 (AAB37265.1), which is an S-RNase-like gene, and with Hordeum vulgare subsp. vulgare rsh1 (AAF45043.1), which is exclusively expressed in young leaf tissues. The genes within Class-I orchid RNase-T2s II and Class-I orchid RNase-T2s III grouped with Nicotiana alata RNase NE (AAA21135.1), which is not linked to the SI locus and not specifically expressed in styles, and Prunus dulcis RNase PD2 (AAG31930.1), which is predominantly expressed in petals, the pistils of open flowers, and leaves. Class-II orchid RNase-T2s were a sister group to Solanum lycopersicum RNase LER, which is not specifically expressed in styles (Kothke and Kock, 2011); Antirrhinum majus × Antirrhinum hispanicum AhSL28 (an S-RNase-like gene), which is not only expressed in pistils; and Nicotiana glutinosa RNase NGR2, which is constitutively expressed in leaves. Overall, the orchid S-RNase-like genes clustered with S-RNase-like genes that were not specifically expressed in styles (Figure 4, Supplementary File S2, and Table S2), suggesting that these genes might not be involved in orchid SI.