Ping’ou hybrid hazelnuts are hybrids with excellent traits bred from several advanced C. heterophylla trees from northeast China and the mixed pollen of several C. avellana seedlings introduced from Italy in the 1980s (Liang et al., 2012). Being the main cultivars of Ping’ou hybrid hazelnut, both “Dawei” (breeding code: 84-254) and “Liaozhen No. 7” (82-11) were grown in the germplasm repository located in Yanqing District, Beijing, China. The specific S-alleles are currently unknown for both the cultivars. Nonetheless, “Dawei” can display SI and exhibits adequate compatibility in reciprocal crosses with “Liaozhen No.7” according to field pollination experiments. Therefore, here we designed experiments for the pollination of self-incompatible (“Dawei” × “Dawei”) and cross-compatible (“Dawei” × “Liaozhen No. 7”) plants to analyze the dynamic change of pollination and to identify the likely possible candidates for the SSI system of Corylus.

One-year-old branches, each ca. 60–80 cm in length, with well-developed male inflorescence (catkins) were collected from “Dawei” and “Liaozhen No.7” trees in February 2021. To avoid contamination, these branches were cultured in water and stored in two air-proofed incubators at ca. 25°C. Pollen grains were carefully collected from the catkins within 48 h of their shedding and stored at −20°C. In addition, another group of branches of “Dawei” was emasculated and water-cultured to boost female flowers at 25°C under 60% relative humidity without any contamination from other pollen. Artificial pollination events were then performed when the styles had protruded 3–4 mm by using a cotton-swab to brush the pollen of “Liaozhen No.7” and “Dawei,” respectively. The styles of the two combinations were collected at 10, 30, and 60 min after each pollination event. In this way, seven kinds of style pools were generated in this experiment: (1) unpollinated styles, serving as the control samples (CK); (2) styles at 10 min after a self-incompatible pollination (IC10); (3) styles at 30 min after self-incompatible pollination (IC30); (4) styles at 60 min after self-incompatible pollination (IC60); (5) styles at 10 min after cross-compatible pollination (C10); (6) styles at 30 min after cross-compatible pollination (C30); (7) styles at 60 min after cross-compatible pollination (C60). Each pool was prepared with three biological replicates. Their styles were carefully dissected from the flower bud, using two fine-tipped tweezers, and these samples were immediately frozen in liquid nitrogen. For each replicate, styles of at least 30 female flowers were mixed together for a total weight of 0.3 g. All the samples were stored at −80°C before their sequencing was carried out.

Total RNA was extracted from the seven style pools’ samples by using a cetyl trimethylammonium bromide (CTAB) method (Gambino et al., 2008). The quality of total RNA was assessed using a NanoDrop 2000 (Thermo Fisher Scientific, Wilmington, DE, United States) and further examined by agarose gel electrophoresis. Construction of cDNA libraries and their respective RNA sequencing were both performed by Biomarker Technologies Co., Ltd. (Beijing, China). The cDNA libraries were constructed using NEBNext Ultra™ RNA Library Prep Kit for Illumina (NEB, San Diego, CA, United States) and following the manufacturer’s recommendations. Sequencing was done on an Illumina NovaSeq sequence analyzer to yield 150-bp paired-end reads. To generate a high-quality clean data set, any reads of low-quality or having adapters or poly-N were first removed using in-house Perl scripts. Next, values of Q20, Q30, GC-content, and sequence duplication level were calculated. The clean reads were then aligned to the reference genome of C. heterophylla (PRJNA655406) (Zhao et al., 2021), using the Hisat2 tool (Kim et al., 2015). Only perfect matches or reads with one mismatched base were included in the subsequent formal analyses. All raw sequencing data have been deposited (accession no. PRJNA763748) into the NCBI Sequence Read Archive¹.

Relative expression levels of genes in each sample were estimated by the value of fragments per kilobase of transcript per million fragments mapped (FPKM). Comparisons between 13 groups (Table 1) were made to identify those differentially expressed genes (DEGs) possibly involved in the pollen–pistil interactions of Corylus. DEGs were designated using the edgeR method with dual criteria of a false discovery rate (FDR) ≤ 0.05 and a fold-change ≥ 1.5 (Ma et al., 2018).

Functional enrichment analyses were then performed for these DEGs. The Gene Ontology (GO) analysis was implemented using “Goseq” package for R, which is based on the Wallenius non-central hypergeometric distribution (Young et al., 2010). The Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of DEGs was conducted using the KOBAS software (2.0) (Mao et al., 2005). Heatmap analysis was performed using BMKCloud².

The S-locus region has been identified in the genome of C. avellana L. It spans 193.5 kb and contains 18 genes, including those encoding three transmembrane receptor-like serine/threonine protein kinases and five leucine-rich repeat receptor-like protein kinases (LRR-RLK) (Hill et al., 2021). Markers (Figure 1) that were previously identified from C. avellana (Gürcan and Mehlenbacher, 2010; Hill et al., 2021) were selected here as targets to perform a local BLAST analysis in the C. heterophylla genome. The regions covered by the start-end location of homologous sequences of markers were considered here as the predicted corresponding S-locus region in the C. heterophylla genome. We focused on screening genes predicted to occur in the homologous S-locus region and analyzing their respective expression levels. Meanwhile, genes expressed inconsistent pattern in comparisons of CK (CK1, CK2, and CK3) vs. IC (IC10, IC30, and IC60) groups and C (C10, C30, and C60) vs. IC (IC10, IC30, and IC60) groups were identified from a Venn diagram analysis.

A total of 18 DEGs were selected for qRT-PCR to validate the relative expression data generated from RNA-Seq. Their first-strand cDNA was synthesized using a BioRrad iScript cDNA Synthesis Kit (Bio-Rad, Hercules, CA, United States) and these were used as a template for the qPCR analysis. Each reaction was conducted in a 20-μl total volume mixture according to the manufacturer’s instructions of Bio-Rad iTaq™ Universal SYBR Green Supermix (Bio-Rad, Hercules, CA, United States). Thermal cycler parameters were set at 95°C for 3 min, 95°C for 30 s, then 40 cycles of 95°C for 5 s and 50–60°C for 5 s. The qRT-PCR analysis with three biological replicates was conducted with the Bio-Rad CFX Manager version 3.0 software (Bio-Rad, Hercules, CA, United States). The relative expression level was calculated using the 2^–ΔΔCt method (Livak and Schmittgen, 2001) with ChaActin and ChaEF1-α serving as the reference genes (Hou et al., 2021).

Transcriptome sequencing generated 144.95 Gb of high-quality clean reads. Among them, 94.38% had a Q-score larger than the Q30 value. The mapping ratio of each sample against the C. heterophylla reference genome ranged from 87.77 to 88.99%. These results indicated robust sequencing and mapping in our study. Finally, we obtained 27,591 genes annotated in the reference genome and 2,335 novel transcripts without direct annotation. We further evaluated the correlations for different biological replicates, which indicated a high within-group similarity for the samples (Supplementary Figure 1).

Through the comparison of 13 groups, we identified DEGs that met the threshold of FDR ≤ 0.05 and fold-change ≥ 1.5 (Figure 2A). The maximum number of DEGs was obtained in the CK vs. C60 comparison group, while the minimum number of DEGs was found in IC10 vs. IC30. In total, we tallied that 19,163 DEGs were upregulated and 13,314 were downregulated. Then we conducted a global comparison among the unpollinated (CK1, CK2, and CK3), crosscompatible pollination (C10, C30, and C60), and self-incompatible pollination (IC10, IC30, and IC60). The comparison group CK vs. IC generated 4,279 (2,765 upregulated and 1,514 downregulated) DEGs (Figure 2B), these mainly being pollen-specific ones and those induced by self-incompatible pollination. There were 3,927 (2,589 upregulated and 1,338 downregulated) DEGs from the comparison group CK vs. C (Figure 2C), these consisted of pollen-specific DEGs and those induced by cross-compatible pollination. Comparison group C vs. IC yielded 840 DEGs (101 upregulated and 739 downregulated) induced by the difference between incompatible and compatible pollination events (Figure 2D). We performed the analysis of changed expression of DEGs in these three comparison groups (Figure 3). Evidently, there were significant differences among the pollination treatments. Taken together, these results can provide new sources to better understand the pollination process.

In the GO annotation of all DEGs, the greatest number of enriched genes was for those related to cellular process, metabolic process, membrane, membrane part, binding, and catalytic activity (Figure 4A). We drew a topGO acyclic graph (Figure 4B) to show the 10 most enriched terms for each major category. Among them, we focused on the term “recognition of pollen” (0048544, Figure 4C) containing 71 DEGs, of which 59 are homologous with the G-type lectin S-receptor-like serine/threonine-protein kinase, another 7 DEGs are receptor-like serine/threonine-protein kinase, 3 DEGs are homologous to a putative receptor protein kinase (Zmpk1), 1 DEG is a hypothetical protein FH972_006764 from Carpinus fangiana, and 1 DEG is a DETOXIFICATION 34 (AtDTX34). In the differential comparison groups of CK vs. C (CK vs. C10, CK vs. C30, CK vs. C60), CK vs. IC (CK vs. IC10, CK vs. IC30, CK vs. IC60), and C vs. IC (C10 vs. IC10, C30 vs. IC30, C60 vs. IC60), “recognition of pollen” remained one of the most significantly enriched terms. Likewise, a vast majority of DEGs were also homologous with G-type lectin S-receptor-like serine/threonine-protein kinase. However, another GO term (0009860), “pollen tube growth,” was significantly enriched in the comparison of C and IC groups. It contained five protein kinesin light chain-related 2 (KLCR2) homologous genes, one CRIB domain-containing protein RIC1 (RIC1) homologous gene, two pollen-specific leucine-rich repeat extensin-like protein 3 (PEX3) genes, and one myosin-11 homologous gene.

We also carried out the pathway enrichment for all DEGs and respectively those from the comparison groups (CK vs. IC, CK vs. C, and C vs. IC). The plant–pathogen interaction, plant hormone signal transduction, and MAPK signaling pathway-plant were significantly enriched for all DEGs and the DEGs from the comparison of CK vs. C and CK vs. IC (Figure 5A); this suggested the pollination process was mainly associated with those pathways. Furthermore, the pentose and glucuronate interconversions, circadian rhythm-plant, phenylalanine metabolism, and phosphatidylinositol signaling system enriched the largest number of DEGs in the comparisons of C vs. IC (Figure 5B). DEGs involved in these major pathways included but were not limited to pathogenesis- and disease resistance-related genes, various receptor kinases, calcium-related proteins, and WRKY transcription factors. They were upregulated after self-incompatible and compatible pollination (Figure 6).

To further identify those DEGs associated with self-incompatible pollination in Corylus, we identified DEGs common to every comparison group involved in self-incompatible pollination (i.e., CK vs. IC10, CK vs. IC30, CK vs. IC60, C10 vs. IC10, C30 vs. IC30, and C60 vs. IC60). We designed six combinations of all different comparison groups (Table 2) and from them identified a total of four common upregulated and five common downregulated DEGs. The upregulated DEGs shared by all combinations were EVM0007329, EVM0011085, EVM0012084, and EVM0016921 (Figure 7A); the shared downregulated DEGs for all combinations were EVM0006306, EVM0008280, EVM0010363, EVM0025402, and EVM0025530 (Figure 7B). The expression changes of common genes are shown in Supplementary Figure 2, and their annotation information can be found in Supplementary Table 1. These genes may play a pivotal role in the process of incompatible pollination. However, none of these DEGs were located within the predicted S-locus. Additionally, in Figure 7 can be seen the unique DEGs for each combination.

A series of molecular markers narrowed the S-locus to 193.5 kb in C. avellana. The S-locus is now flanked by marker RH_SLOC07 on the left and KG847 on the right (Hill et al., 2021). Based on that study, we obtained all molecular markers used for fine-mapping the S-locus in Figure 1. Specifically, we used these molecular markers as targets to perform the local BLAST against the C. heterophylla genome. We then confirmed the ordered location of these molecular markers in the C. heterophylla genome (Supplementary Table 2). The homologous gene for the right-side marker KG847 of the S-locus region was located at 3,722,377–3,723,320 bp on LG05 of C. heterophylla. Yet, we failed to identify the left-side marker RH_SLOC07, in that no significant matches were found, perhaps because of interspecific differences. Finally, the predicted S-locus was located at 3,500,000–3,800,000 bp on LG05 by using the right-side marker KG847 and the left-side marker 846-HRM1. We further identified those genes located in the predicted S-locus in terms of their location, functional annotation, and gene expression intensities (Figure 8). The expression of EVM0026696 and EVM0001920 was zero in all treatments, and they were not included in the heat map analysis. Among S-locus genes, eight genes showed homology with those identified in the S-locus of C. avellana, namely two PIX7s (EVM0002129 and EVM0025536), three MIK2s (EVM0002422, EVM0005666, and EVM0009820), one aldose 1-epimerase (EVM0002095), one 3-dehydroquinate synthase II (EVM0021283), and one At3g28850 homologous gene (EVM0016149). They are the most promising candidate genes for explaining SSI functioning in Corylus.

Furthermore, we compared the expression of those eight genes in the different comparison groups, finding no significant differences for EVM0002129 and EVM0021283 in any groups. EVM0025536 showed downregulated expression pattern in the comparisons of CK vs. C60 and CK vs. IC60; however, it had no significant differences in other comparison groups. By contrast, EVM0002422 and EVM0016149 were each upregulated in CK vs. C and CK vs. IC comparison groups, although no difference was found for C vs. IC. For EVM0009820, although it was upregulated in CK vs. IC30/IC60 and CK vs. C30/C60, its expression was negligibly changed in other comparison groups. EVM0005666 showed consistent expression in CK vs. C10 and C vs. IC, but it was downregulated in CK vs. IC10, yet upregulated in CK vs. C30/C60 and CK vs. IC 30/IC60 comparisons. Further, upregulated expression of EVM0002095 was detected in CK vs. C and CK vs. IC, whereas it was downregulated in C vs. IC.

A previous study suggested that MIK2 and PIX7 homologous genes are the most likely candidates in Corylus for studies of its SSI (Hill et al., 2021). Our conserved domain analysis revealed that EVM0025536 had three characteristics. The first domain was homeodomain involved in the transcriptional regulation of key developmental processes (Alonso, 2002), the second domain was homeobox-associated leucine zipper, and the third was catalytic domains of serine/threonine and tyrosine protein kinases (STYKc) (Figure 9A). Yet, EVM0002129 has only a STYKc domain (Figure 9A). The three MIK2 homologous genes (EVM0009820, EVM0005666, and EVM0002422) all harbor leucine-rich repeats (LRR) and have the same STYKc domain that PIX7 has (Figure 9B).

Eighteen DEGs were used for a qRT-PCR analysis to independently check the relative expression levels detected by RNA-Seq. These genes included 8 S-locus genes common to C. heterophylla and C. avellana, EVM0000218 (serine decarboxylase), EVM0000333 (probable xyloglucan endotransglucosylase/hydrolase protein 23), EVM0002634 (Wound-induced protein 1), EVM0002247 (adenylate isopentenyltransferase 3), EVM0000979 (mitochondrial uncoupling protein 5), EVM0000443 (Expansin-B3), EVM0002638 (Cytochrome b561 and DOMON domain-containing protein At3g25290), EVM0008735 (stress-induced protein KIN2-like), EVM0005669 (sphinganine C4-monooxygenase 1 homologous gene), and EVM0000081 (pectate lyase homologous gene). These qPCR results indicated a strong consistency with those from the RNA-Seq analysis (Supplementary Figure 3), thus confirming the reliability of our RNA-sequenced data and its results.

The acceptance of compatible pollen and the rejection of self-incompatible pollen by pistil is an important step for pollen–pistil interactions at a very early stage. Pollen germination occurs within 15 min after pollen grains land on the stigma in Arabidopsis thaliana (Mayfield and Preuss, 2000); pollen germination begins within 2–3 min, and pollen tubes reach the bottom part of the ovary within 4 h in rice (Lan et al., 2004); and pollen grains germinate within 30 min after landing on the stigma in Chinese cabbage (Jiang J. J. et al., 2013). However, researchers did not observe any pollen attachment and pollen-tube growth on the stigma surface at 30 min and even 6 h after SI–pollination in W1 canola; however, pollen tubes could be observed in transmitting tissue at 6 h after compatible pollination, although no pollen attachment was discernible at 30 min (Sankaranarayanan et al., 2013). In contrast in non-heading Chinese cabbage (Brassica campestris ssp. chinensis Makino), a limited number of pollen grains could adhere to the stigmatic surface within 0.25 h after incompatible pollination. Abundant callose was observed in the incompatible stigma epidermal cell, which could prevent pollen tubes from penetrating the stigmatic surface (Wang et al., 2014). These findings are consistent with our recent study on pollen–pistil interaction in Corylus. We observed secretion and pollen adhesion on the stigma surface and blue fluorescence on the papillae within 60 min after compatible and incompatible pollination. While in incompatible pollination, bunched, short, or bulbous pollen tubes which do not penetrate the stigmatic surface were observed after 4 h (Li, 2020; Hou et al., 2022). These observations suggested that pollen adhesion, hydration, and penetration of pollen tube are the earliest steps blocked in self-incompatible pollination. Following the initial adhesive interaction of pollen on the stigmatic surface, there is a latent period of 30 min when signals are exchanged between pollen coat proteins and the stigmatic components (Chapman and Goring, 2010; Sankaranarayanan et al., 2013). Deciphering the transcriptional changes during this period would reveal genes involved in compatible and self-incompatible responses. Therefore, sampling pistil tissue at 10, 30, and 60 min after pollination is deemed crucial to rigorously explore this question through transcriptome analysis, as done in this study for hazelnut (Corylus).

We found several types of genes, including pathogenesis- and disease resistance-related genes, various receptor kinases, calcium-related proteins and WRKY transcription factors, that were up-regulated expression after self-incompatible and compatible pollination. The role played by these genes in Corylus is not yet known, but they are important candidates for regulating pollination of further studies.

In plants, Ca²⁺ is an important second messenger functioning during the entire pollination process (Denninger et al., 2014). During the SSI response of Brassica, the direct interaction between the male- and female-determinants, SP11/SCR and SRK results in a drastic increase of Ca²⁺ in the papilla cell (Iwano et al., 2015). This finding showed that Ca²⁺ influx into stigma papilla cells mediates SI signaling. In GSI species, for example, Pyrus pyrifolia (Jiang X. T. et al., 2013), calcium is required for the normal growth of pollen tubes. In our study, calcium-related genes were significantly enriched in the major pathway, while many calcium-related proteins were also upregulated after self-incompatible pollination. These results indicated that calcium-related genes could function as potential signal factors in regulating the SI pollination process.

The WRKY transcription factors are essentially integral parts of signaling webs that modulate many plant processes. The WRKY34 transcription factor has been proven to negatively mediate cold sensitivity in mature pollen of Arabidopsis, and the overexpression of WRKY34 and WRKY2 much reduced the pollen viability, pollen germination, and pollen tube growth (Zou et al., 2010; Guan et al., 2014). Cotton WRKY transcription factor, GhWRKY22, acts as a transcriptional repressor to regulate anther/pollen development possibly by modulating the expression of the JAZ genes (Wang et al., 2019).

Receptor-like protein kinases genes are important components of proper plant growth and development, playing key roles in defense against pathogens, morphogenesis of various tissues and organs, signal transduction, and cell-cell recognition (especially in the case of SI). In flowering plants, RLKs are essential not only for the recognition of pollen, but also for the guidance and reception of the pollen tube. In response to pollen adhesion, plants are able to rapidly generate new transcripts, and this generates a highly sensitive compensatory mechanism that can be activated to reject/accept pollen as needed. That there were a considerable number of upregulated RLK genes in our transcriptional analysis highlights the relevance of RLKs in pollen recognition and rejection. In plant species, there is a great diversity and abundance of RLKs, and they contain the well-known leucine-rich repeat (LRR), lysin-motif (LysM), and S-domain (SD) RLKs as potential pattern-recognition receptors (Antolín-Llovera et al., 2012). Some members of SD–RLKs are involved in SI, like the S-locus receptor kinase (SRK) and S-locus glycoprotein (SLG) for the SLG domain, this is a well-proven determinant for SSI in Brassica (Tantikanjana, 1993). Nevertheless, SD-RLKs also occur in species that do not possess SI, wherein they are upregulated in response to wounding, pathogen recognition, or enemy attacks. This suggests that evolution has driven the expansion of specific RLK families to serve multiple roles in different plant physiological processes (Sanabria et al., 2008). Here we found many SD-RLKs enriched in the KEGG pathway of plant–pathogen interaction in our transcriptome results. But these genes are not situated in the S-locus; hence, they may have other potential functions. Three LRR–RLKs and two serine/threonine protein kinases (RLK) were identified in the predicted S-locus. According to the SSI system of model plants, the involvement of RLKs and LRR-RLKs in pollen–pistil interactions might be through the recognition and binding of ligands to initiate both auto- and trans-phosphorylation of itself and other participating proteins.

Pollen–pistil interactions may be analogous to host–pathogen interactions (Wilkins et al., 2014). Plants have evolved unique defense mechanisms enabling them to perceive pathogens and initiate effective defense strategies. Both SI and the immunity system evolved under different selection pressures, namely, the former case prevents inbreeding and the latter parasitism. The recognition and rejection of self-pollen are remarkably similar to the enemy recognition and defense activation in plants. In the SI system, the pistil recognizes and responds to the self-pollen. By contrast, in the immunity system, non-self-ligands from the pathogen are recognized by RLKs. The SI response and plant–pathogen interactions have close similarities, such as the “host” penetration by a tubular cell emanating from a spore-like structure (Hodgkin et al., 1988) and the striking structural similarity between SRK genes and the wheat leaf rust kinase (WLRK) defense gene (Feuillet et al., 1998). Likewise, in our transcriptome analysis, many differential genes mediating the pollen–pistil interaction were also enriched in plant–pathogen interaction, such as numerous pathogenesis- and disease resistance-related genes and the WRKY transcription factor. This is consistent with pollination studies in other species, for example, Arabidopsis (Swanson et al., 2005) and Asteraceae (transcriptomic comparison of self-pollinated and cross-pollinated flowers of Erigeron breviscapus to analyze candidate SI-associated genes). These results indicated that both the immunity system and the SI systems rely on receptors to recognize ligands.

We identified two PIX7s, three MIK2s, one aldose 1-epimerase, one 3-dehydroquinate synthase II, and one At3g28850 homologous gene in the predicted S-locus of C. heterophylla. These results are consistent with the recent findings of Hill et al. (2021). Initially, we had expected that some of these genes could be upregulated in all differential comparison groups that involved incompatible pollination (i.e., CK vs. IC and C vs. IC). However, transcriptome results for the expression of these genes did not support our original hypothesis. Lack of upregulation of SI-specific genes and downregulation of 19 proteins following self-incompatible pollination were also found in Brassica (Samuel et al., 2011; Sankaranarayanan et al., 2013). Nevertheless, further research is still needed to verify the functioning of those genes.

An aldose 1-epimerase, EVM0002095, is a key enzyme of carbohydrate metabolism; it catalyzes the interconversion of the alpha- and beta-anomers of hexose sugars, such as glucose and galactose. The 3-dehydroquinate synthase II (EVM0021283) was isolated from the archaeon Methanocaldococcus jannaschii and demonstrated to play a key role in an alternative pathway for the biosynthesis of 3-dehydroquinate (DHQ) (White, 2004). EVM0016149 is an At3g28850 homologous gene whose protein product contains a glutaredoxin (GRX) domain. In Arabidopsis, GRXs are involved in petal development and salicylic acid signaling (Rouhier et al., 2008). PIX7 encodes receptor-like cytoplasmic kinases (RLCK) belonging to subfamily VII of receptor-like kinases (RLK), and it is capable of interacting with both wild-type and mutant XopAC forms in a yeast two-hybrid screen (Guy et al., 2013). The leucine-rich-repeat RLK MIK2 is a male perception to the LURE1 peptides for the specific binding between the extracellular domain of MIK2 and LURE1 (Wang et al., 2016). By contrast, tip-focused pollen-specific receptor-like kinase 6 (PRK6) is an important receptor for perceiving the LURE1 (Zhang et al., 2017). Moreover, the C-terminal loop of the LRR domain of AtPRK6 could bind to AtLURE1 (Takeuchi and Higashiyama, 2016). MIK2, being the receptor for the SERINE RICH ENDOGENOUS PEPTIDE (SCOOP) (Rhodes et al., 2021), directly binds with SCOOP12 and triggers a complex formation between BRASSINOSTEROID INSENSITIVE 1-ASSOCIATED KINASE 1 (BAK1) and MIK2. Notably, MIK2 fosters immunity against Fusarium pathogens via recognition of Fusarium-derived SCOOP-like sequences. Meanwhile, the serine-rich SCOOP as a MIK2 ligand implies an interaction between MIK2 and PIX7, as the latter also has a serine/threonine domain. These results provide insight into the participation of MIK2 and PIX7 in the SSI system of Corylus. Therefore, it is not surprising to find that those MIK2 and PIX7 homologous genes could be involved in the SSI of Corylus. Yet the respective functions of these genes in Corylus remain unclear. Nevertheless, we propose that they are promising candidates for controlling SSI.

In Brassica, the S-locus-related glycoprotein 1 (SLR1), a stigma-specific protein, was proved to be responsible for the adhesion of pollen grains to the stigmatic surface. Two SLR1-binding proteins, SLR1-BP1 and SLR1-BP2, were determined to operate as the interacting counterpart of SLR1 in pollen (Takayama et al., 2000). However, none of the sequences showed homology to SLR1 in C. heterophylla’s genome and there was negligible to zero expression of 3 SLR1-BP in every transcriptome sample. EXO70A1 is part of the exocyst complex, carrying out a compatibility function along with annexin and actin in delivering vesicles to the site of pollen attachment. The vesicle contains the majority of resources required for pollen germination. Accordingly, suppression of EXO70A1 and reduced expression of both annexin and actin would be implicated in the SI response, whereby vesicle delivery to the pollen attachment site is inhibited (Samuel et al., 2011). However, no expression of an annexin gene (EVM0010564) was detected in all IC samples, and the expression of an actin gene (EVM0010391) was neither upregulated nor downregulated in all comparison groups. Moreover, our previous study had selected ChaActin as a suitable reference gene to evaluate gene expression patterning for the pollen–pistil interaction in Corylus. Our results provide further compelling evidence that Corylus displays a unique SSI mechanism, and they should prove helpful to investigate this potential mechanism of the SSI system in Corylus spp.

In general, more remains to be explored concerning the plant genes involved in pollen–pistil interactions. Deciphering the underlying molecular mechanisms behind the events pertaining to SSI responses in Corylus and other plants is still a challenging field of research. In this study, a range of candidate genes was obtained from the analyzed transcriptome data, which enhances the essential data available for understanding the Corylus pollen–pistil interaction. Further research will involve whole-sequence analysis, characterization of temporal-spatial expression, and gene–gene interaction analyses.