The Molecular Dialog between Flowering Plant Reproductive Partners Defined by SNP-Informed RNA-Sequencing

A SNP Database to Determine Whether RNA-Seq Reads from Pollinated Pistils Are Derived from the Pollen or Pistil Parent

There are ∼100 to 200 pollen tube cells and thousands of pistil cells in a pollinated pistil, so we estimated ∼1% of the pollinated pistil transcriptome would be pollen tube derived. To maximize the likelihood that we would detect RNA-seq reads unambiguously defining pollen tube transcripts, we chose two A. thaliana accessions, Col and Cvi, with a high single nucleotide polymorphism (SNP) frequency (∼1 SNP every 200 bp; Schmid et al., 2003; Nordborg et al., 2005). We created a custom database of Cvi SNPs that had been validated by genome sequence (Gehring et al., 2011; Nodine and Bartel, 2012; Pignatta et al., 2014). To capture the greatest number of SNPs in expressed genes from Col and Cvi, we incorporated additional SNPs identified from RNA-seq of unpollinated and self-pollinated Cvi pistils (8 h after hand self-pollination; Figure 1A). Variant analysis of unpollinated and self-pollinated Cvi pistils added 15,995 new SNPs to the previously identified 576,722 Cvi SNPs (Pignatta et al., 2014), resulting in a comprehensive variant list used in all subsequent analyses (1 SNP every 201 bp; Supplemental Data Set 1). As 36.2% of these SNPs occur in exons (coding sequences and untranslated regions), they are useful for determining the parent of origin of individual RNA-seq reads from pollinated pistils (Figure 1B; Supplemental Figure 1). Importantly, 75% of annotated A. thaliana genes had at least one SNP in a transcribed gene region, providing the potential to quantify pollen tube and pistil transcript abundance unambiguously (24,441/32,678 = 74.7%).

• Download figure
• Open in new tab
• Download powerpoint

Figure 1.

Pollen- and Pistil-Specific Gene Expression Can Be Distinguished in Pollinated Pistils Using SNPs in RNA-Seq Reads.

(A) Experimental samples prepared for RNA-seq.

(B) The distribution of SNPs across indicated gene regions in our Cvi versus Col SNP database (percentage is given above bars; black represents transcribed regions; gray represents untranscribed regions of transposable elements [TE]).

(C) Reads from three biological replicates of a pollinated pistil RNA-seq experiment were analyzed to determine the number with Col or Cvi SNPs; the numbers of total and sex-specific reads are given. Venn diagram is not proportional to number of included elements.

(D) and (E) Scatterplots of mean reads/gene (three biological replicates).

(D) Pistil reads (y axis) are plotted against pollen reads (x axis) for all genes with SNPs detected in pollinated pistils. The y axis is 10-fold greater than the x axis.

(E) Total reads/gene are plotted for unpollinated pistils (y axis) and pollinated pistils (x axis). All genes are plotted as colored dots indicating sex-specific expression: blue, pollen specific; red, pistil specific, gray, genes without SNPs.

(F) Quantification of genes expressed in pollen, pistil, or both tissues. Supporting bioinformatic analysis is shown in Supplemental Figures 1, 3, and 4.

Eight Percent of Total mRNA within the Pollinated Pistil Is Derived from the Pollen Tube

We used bioinformatic approaches previously developed to define allele-specific gene expression in the developing embryo and endosperm (Gehring et al., 2011; Pignatta et al., 2014) to determine whether RNA-seq reads contained SNPs and whether reads with SNPs were derived from the pollen tube or the pistil (Supplemental Figure 2, analysis workflow; see Methods). We used methods that count reads mapping to a single genome location (HTseq) or that count reads mapping to multiple locations (EasyRNASeq; Supplemental Figure 2, Supplemental Data Sets 2 and 3, comparison of methods). In pollinated pistils, 12% of 131 million singly mapped reads analyzed contained SNPs (Figure 1C), and these reads were distributed across the transcribed regions of genes similarly to reads without SNPs (Supplemental Figure 3 and Supplemental Data Set 2). Inclusion of reads that were mapped to multiple locations did not significantly alter the proportion of reads with useful polymorphisms (Supplemental Data Set 2); however, multiply mapped reads facilitated enhanced detection of members of gene families expressed in pollen. So, multiply mapped reads were used for subsequent analyses.

These data suggest that in a typical experiment, ∼17 million reads with SNPs are available to quantify patterns of transcript abundance in the pollen tube and pistil. Importantly, we determined that ∼60% of genes expressed in the pollinated pistil have at least one SNP, allowing unambiguous determination of parent of origin (Figure 1) in crosses between Col and Cvi. This analysis shows that we will not be able to unambiguously determine the parent of origin for all genes using Cvi pistils pollinated with Col pollen. Further bioinformatic analysis of sequenced A. thaliana accessions showed that crosses between Col and Pu2-8 would provide the greatest number of additional useful SNPs (∼20% more genomic features; Supplemental Figure 4).

To determine the relative contributions of pollen tube and pistil to the pollinated pistil transcriptome (Figure 1A), we analyzed the total number of reads from Cvi pistils pollinated with Col pollen that contained SNPs (Figure 1C) and found that 92% of these reads had the Cvi allele (pistil reads), while 8% had the Col allele (pollen tube reads; Figure 1C). These data show that pollen tube RNA constitutes ∼8% of the total RNA extracted from a pollinated pistil. Moreover, a scatterplot comparing the number of pistil reads (Cvi, y axis) with pollen tube reads (Col, x axis) for all genes with SNPs shows that pistil transcripts are much more abundant for the vast majority of genes (y axis is 10-fold greater than x axis). This is not surprising given the abundance of pistil cells compared with pollen tubes. However, genes with expression values along the x axis point to a significant fraction only expressed from the pollen tube and readily detected using this method (Figure 1D). Analysis of these data defines groups of pollen-expressed and pistil-expressed genes.

Sex-Specific Transcript Analysis Reveals the Parent of Origin of Genes Expressed during Pollination

A scatterplot comparing the number of total reads/gene in unpollinated pistils with pollinated pistils identifies genes with expression patterns that changed upon pollination and highlights a group of genes that are only detected upon pollination (Figure 1E, points along x axis). The simplest hypothesis is that these genes are expressed in the pollen tube; however, it is possible that some of these genes are pistil-expressed genes induced in response to pollen tube growth. To distinguish between these possibilities, we color-coded the scatterplot as follows: pollen tube-expressed genes are shaded blue (>10 average pollen reads), pistil-expressed reads are shaded red (>10 average pistil reads), and genes expressed in both tissues (>10 average pollen and pistil reads) are assigned a color between red and blue based upon the ratio of pistil to pollen reads (Figure 1E; Supplemental Data Set 3). Pollen gene expression is primarily responsible for the group of genes that are highly induced during pollination (Figure 1E, blue dots).

To define pistil genes that responded to pollination, we focused on reads containing Cvi SNPs in unpollinated (Cvi) versus pollinated pistils (Cvi x Col; Supplemental Figure 2). We used the DESeq2 software package (Love et al., 2014), which uses read counts and estimates variance-mean dependence and tests differential gene expression using the negative binomial distribution and does not normalize expression to the length of genes within libraries. This feature of DESeq2 was advantageous because we wanted to compare the expression of individual genes across conditions using the subset of reads that included SNPs. While a majority of pistil-expressed genes remain unchanged by pollination (Figure 1E, red dots along diagonal), it is clear that a number of pistil-expressed genes (pistil reads only) are differentially expressed in response to pollination; 476 genes are significantly below the diagonal (induced by pollination) or above the diagonal (higher in unpollinated pistils, DESeq2; see Methods; Supplemental Data Set 4) and could play roles in modulating the pistil environment to promote directed pollen tube growth.

To identify pollen-specific and pistil-specific genes, we normalized the number of sex-specific reads over each gene to the number of SNPs/gene. Then, we applied a 2-fold enrichment criterion for pollen-specific genes: >20 Col reads/SNP/gene, and <10 Cvi reads/SNP/gene (average of three replicates). The reciprocal criteria were used to identify pistil-specific genes. Normalization by SNPs/gene accounts for the large variation in the SNP frequency across genes (1 to >300; Supplemental Data Set 3). These stringent thresholds identified 277 genes expressed specifically from the pollen tube and 6838 genes expressed specifically from the pistil (Figure 1F). The ∼25-fold higher number of pistil-specific genes is due, at least in part, to the vast excess of pistil tissue.

To determine whether SNP-informed determination of sex-specific gene expression identifies genes known to be expressed in the pollen tube or pistil, we analyzed genes defined in published microarray data. We graphed read coverage over genes and color-coded reads to indicate whether they contained a pollen- or pistil-derived SNP (red, pistil SNP; blue, pollen SNP) or whether they had no SNP (gray, Figures 2A to 2F). AT1G47530, AT2G13690, and AT5G46800 had been shown to be abundant in pistils (Boavida et al., 2011), but absent or minimally detectable in pollen grain and pollen tube microarray experiments (Becker et al., 2003; Honys and Twell, 2003; Pina et al., 2005; Qin et al., 2009). All of the reads with SNPs carried the Cvi allele for these genes (Figures 2A to 2C, red), indicating that in pollinated pistils all of the RNA transcribed from these genes is transcribed in pistil cells.

• Download figure
• Open in new tab
• Download powerpoint

Figure 2.

SNP-Informed RNA-Seq Defines the Male and Female Contributions to Gene Expression during Reproduction.

(A) to (F) Plots of reads over gene models at representative loci. Coverage (y axis) = log₁₀ number of reads counted over each nucleotide. Gray color represents reads without SNPs. In regions with SNPs (green line), reads are divided into sex-specific components: red = reads with the Cvi SNP (pistil parent) and blue = reads with the Col SNP (pollen parent).

(A) to (C) Pistil-expressed genes.

(D) to (F) Pollen-expressed genes.

(G) to (J) Genes that are abundant in pollinated pistils and are expressed in both pollen and pistil. Exons are represented by black bars, introns by black lines, and polarity of transcription by arrowheads.

(K) Pearson’s R correlation values are graphed with the indicated scale comparing pollen tube- and pistil-expressed genes defined by SNP-informed RNA sequencing of pollinated pistils with previously published microarray experiments (references given in text). Only genes with one unambiguous microarray ID were analyzed. Supporting bioinformatic analysis is shown in Supplemental Figure 5.

AT1G29140, AT1G55560, and AT2G26450 were shown to be expressed in SIV-grown pollen tubes (Qin et al., 2009) but were absent in unpollinated pistils (Boavida et al., 2011). We found that all of the reads with SNPs carried the Col allele for these genes (Figures 2D to 2F, blue), indicating that in pollinated pistils all of the RNA transcribed from these genes is transcribed in pollen tubes.

Finally, we took the 10 most highly expressed genes detected in microarray analysis of pollinated pistils (Boavida et al., 2011) and graphed the distribution of reads for SNP-containing genes (8/10 genes contained useful SNPs). Four representative genes are shown (AT4G10340, AT4G35100, AT5G02380, and AT1G72260; Figures 2G to 2J). These genes highlight one of the chief advantages of SNP-informed RNA-seq analysis of reproductive gene expression; they clearly show that the male and female components of gene expression can be determined unambiguously for genes expressed in both the pollen tube and the pistil.

SNP-Informed RNA-Seq Correlates with Previously Defined Reproductive Tissue Transcriptomes

Comparison of sex-specific gene expression defined by SNP-informed RNA-seq with previous microarray analysis of component cell and tissue types indicates overall similarity, but also identifies genes missed in prior analyses. One of the principal advantages of SNP-informed RNA-seq analysis of pollination is that the pollen tube and pistil components of gene expression can be determined without invasive procedures that potentially distort native expression patterns. Prior to this analysis, it was necessary to piece together male and female gene expression patterns from experiments that interrogated pollen (Honys and Twell, 2004; Wang et al., 2008; Qin et al., 2009) and pistil separately (Swanson et al., 2005; Boavida et al., 2011). We performed a Spearman’s correlation analysis to test whether sex-specific transcriptomes determined using SNP-informed RNA-seq correlate with those of component cell and tissue types previously determined by microarray analysis (Figure 2K; Supplemental Figure 5). Correlations with microarray experiments were limited because the lists of gene identifiers are not completely overlapping. Nonetheless, we found that in Cvi x Col pollinations, pollen-expressed genes (reads containing Col SNP) demonstrated a high correlation with pollen transcriptomes, especially SIV-grown pollen tubes (Figure 2K, left column; dark red is the highest Pearson coefficient observed). Conversely, genes from Cvi x Col pollinations identified as pistil-expressed (reads containing Cvi SNP) demonstrated a high correlation with pistil transcriptomes, especially pistils pollinated for 3.5 or 8 h (Figure 2K, right column). As expected, pistil-expressed genes defined by SNP-informed RNA-seq did not correlate with previously identified pollen development, tube growth (in vitro and SIV pollen tubes), or vegetative tissue profiles (Figure 2K, right column). Therefore, SNP-informed RNA-seq analysis provides a noninvasive, in vivo analysis of gene expression that occurs while the pollen tube is growing through pistil tissue and thus provides a more biologically relevant view of gene expression during pollen tube-pistil interactions.

We compared pollen tube-expressed genes defined by SNP-informed RNA-seq analysis with lists of genes previously defined by microarray analysis of SIV-grown pollen tubes (Qin et al., 2009), unpollinated and pollinated pistils (Boavida et al., 2011), and by translating ribosome affinity purification of pollen tubes growing through pistil tissue (Lin et al., 2014). All comparisons are reported in Supplemental Data Set 5. The SIV transcriptome included 7275 expressed genes, whereas SNP-informed RNA seq identified 4375 genes with RNA-seq reads that included the Col SNP (pollen reads); the overlap was 28% (2540 genes). This discrepancy was not due to a lack of SNP coverage: Only 31 SIV genes were not identifiable because they lack a Col-Cvi SNP. Conversely, we identified a much larger group of genes (4704) with useful SNPs that were not detected as pollen tube-expressed genes in vivo but that were detected as pollen tube-expressed genes using the SIV method. Some of these genes may have been induced by the manipulations of the SIV method in which pollen tubes exit pistil tissue and are exposed to air while growing on agar media (Qin et al., 2009). We also identified 1735 genes with clear evidence for pollen tube expression in the pollinated pistil (>10 Col reads) that were not detected by the SIV method. These newly identified pollen tube expressed genes include genes not present on the ATH1 microarray and genes that require native interaction with pistil tissue for their expression.

We found substantial overlap between the 267 pollen tube-specific genes identified by SNP-informed RNA-seq (Supplemental Data Set 5) and 328 previously identified pollination-induced genes (Boavida et al., 2011; overlap = 38%, 123 in common). This is not surprising given that both analyses began with RNA isolated from intact pollinated pistils. By contrast, there was no overlap between the genes defined as pollen tube specific by SNP-informed RNA-seq and by translating ribosome affinity purification (n = 469; Lin et al., 2014; Supplemental Data Set 5). Comparison of all genes identified as pollen expressed in our study (n = 3939) with the pollen tube-specific genes identified by translating ribosome affinity purification identified a modest overlap of 69 genes. The higher degree of overlap between our method and microarray analysis of SIV-grown pollen tubes (Qin et al., 2009) and pollinated pistils (Boavida et al., 2011) is likely because these techniques begin with total RNA rather than RNA associated with a tagged ribosomal subunit that must compete with endogenous ribosomal subunits (Lin et al., 2014).

The Sex-Specific Components of the Pollination Transcriptome Comprise Distinct Functional Categories

We found that pollen tube-specific genes were enriched for Gene Ontology (GO) terms suggesting functions in pectin modification (14-fold; P value = 2.04 × 10⁻⁸, all GO enrichment P values were determined using the Benjamini method), cell-cell signaling (25-fold; P value = 1.4 × 10⁻⁷), calcium-mediated signaling (14-fold; P value = 1.23 × 10⁻⁶), and pollen tube growth (8-fold; P value = 5.5 × 10⁻⁴). Analysis of GO cellular component annotations revealed 27.9% of transcripts encoded genes predicted to enter the endomembrane system (n = 74, P value = 2.72 × 10⁻¹⁵), which includes nine members of the rapid alkalinization factor gene family (27-fold; P value = 1.18 × 10⁻⁷; Supplemental Data Set 6), which have been found to participate in extracellular signaling (Haruta et al., 2014; Du et al., 2016; Stegmann et al., 2017). On the other hand, there were 476 pistil-expressed genes that were differentially regulated upon pollination: 216 were upregulated upon pollination, 259 decreased abundance upon pollination (Supplemental Data Set 7). Secreted glycoproteins (n = 38, P value = 7.25 × 10⁻¹²) and secreted cell wall localized proteins (n = 26, P value = 6.31 × 10⁻⁹) both increase upon pollination as well as a small but significant number of genes involved in response to chitin (n = 4, P value = 0.001) and response to auxin (n = 8, P value = 0.005). These enriched GO categories may indicate parallels between the pistil’s response to pollen tube growth and previously characterized biotic stress responses. By contrast, there is widespread decrease in expression of cytochrome p450 (n = 12, P value = 4 × 10⁻⁶) and oxidative stress genes (n = 10, P value = 4 × 10⁻⁴) upon pollination. Of the 476 pistil-expressed genes, 168 were found to be pistil specific by our criteria, of which 46 are predicted to be secreted (Supplemental Data Set 7). These observed changes in expression likely indicate a major shift in the transcriptional program of female tissue during pollination.

SNP-Informed RNA-Seq Analysis Defines the Role of MYB Transcription Factors in Regulating Pollen Tube Gene Expression and the Pistil’s Response to Pollination

Pollen tubes lacking three closely related MYB transcription factors (Stracke et al., 2001), MYB97, MYB101, and MYB120, fail to release their cargo of two sperm upon entering an ovule. Consequently, myb triple-mutant pollen tubes continue to grow within the female gametophyte producing a coil of pollen tube that fails to fertilize female gametes (Leydon et al., 2013; Liang et al., 2013). It was proposed that pollen tube MYBs respond to the pistil environment and in turn activate the expression of genes critical for pollen tube reception signaling (Leydon et al., 2013, 2014), the interaction between the male and female gametophyte that results in sperm release (Kessler and Grossniklaus, 2011). This model was explored using microarray analysis of triple-mutant pollen grains (Liang et al., 2013) or pollinated pistils (Leydon et al., 2014). It is impossible to know the parent of origin of transcripts using microarray analysis, so neither of these approaches was able to determine the pollen- and pistil-expressed genes with altered transcript abundance in mutant compared with wild-type pollinations. Importantly, microarray analysis failed to identify differential expression of the pollen component of a gene expressed in both the pollen and the pistil because changes in the pollen tube component were masked by the much larger pool of pistil transcripts. SNP-informed RNA-seq provides an opportunity to determine the changes in gene expression due to loss of critical transcriptional regulators of the pollen tube cell as it interacts with the pistil tissue it has evolved to navigate.

First, we compared Cvi x Col pollinations to Cvi x myb triple mutant (accession = Col) pollinations and analyzed differential expression using all mapped reads (Figures 3A and 3B). This analysis is similar to the previous microarray analysis because it ignores the parent-of-origin of transcripts. As expected, RNA-seq identified the majority of the differentially expressed genes previously identified by microarray analysis of pollinated pistils (∼63%, 30/48, DESeq2; Leydon et al., 2013); however, the overlap with microarray analysis of MYB-dependent gene expression in pollen grains was much lower (16%, 4/25, DESeq2; Liang et al., 2013). This difference between the MYB regulatory networks of the pollen grain and the pollinated pistil is consistent with our previous observation that pollen tube MYBs are induced by growth through the pistil and likely regulate unique targets as the pollen tube grows through pistil tissue (Leydon et al., 2013). Importantly, RNA-seq analysis uncovered an additional 288 MYB-dependent genes (318 total, DEseq2; Supplemental Figure 6 and Supplemental Data Set 8; Figure 3B). Genes with significantly lower expression in myb triple mutants were enriched in several functional categories including transmembrane transporters, carbohydrate active proteins, and secreted proteins, as was previously defined (Supplemental Data Set 8; Leydon et al., 2013).

• Download figure
• Open in new tab
• Download powerpoint

Figure 3.

SNP-Informed RNA-Seq Defines MYB-Responsive Genes in the Pollen Tube and Pistil.

(A) Experimental samples prepared for RNA-seq.

(B) and (C) Scatterplots of all average reads (ignoring SNPs) in wild type versus myb triple-mutant pollinations.

(B) Differentially expressed genes identified before division of reads based upon SNPs are indicated (purple, 318 genes, adjusted P value < 0.05).

(C) The color code of the scatterplot in (B) was altered to indicate genes with significant changes in the pollen (reads with Col SNP, light blue, n = 278, P value < 0.05, DESeq2) or pistil (reads with Cvi SNP, pink, n = 34, adjusted P value < 0.05) components of their expression profile.

(D) A comparison of the differentially expressed genes identified with or without SNP analysis (P value < 0.05, DESeq2).

(E) Differential expression analysis was performed again comparing genes identified after analysis of reads without SNPs with those identified based on differential abundance (P value < 0.05, DESeq2) of pollen (Col) or pistil (Cvi) reads. Differentially expressed genes identified by both analyses are highlighted in dark colors, and genes only identified by SNP-informed RNA-seq are highlighted in light colors (blue, pollen reads; red, pistil reads). The total number of differentially expressed genes is shown next to each category, along with the number that were downregulated (down) and upregulated (up) in mutant compared with wild-type pollinations.

(F) Representative read coverage graphs (as in Figure 2) of MYB-dependent genes (AT1G04310 and AT1G33470). Top row, Cvi x Col-0 (wild-type) pollinations; bottom row, Cvi x myb triple-mutant pollinations. Supporting bioinformatic analysis is shown in Supplemental Figure 2. Supporting bioinformatic analysis is shown in Supplemental Figure 6.

Next, we implemented sex-specific read identification on RNA-seq data from wild-type and mutant pollinations to define the patterns of differential expression resulting from loss of MYB97, MYB101, and MYB120 in pollen tubes and the pistil. This analysis identified 278 pollen-expressed genes and 34 pistil-expressed genes that are differentially expressed in myb triple-mutant pollinations (Supplemental Data Set 9). When the mutant versus the wild type total read scatterplot (Figure 3B) is color-coded to highlight genes with significant changes in the pollen tube (blue) or pistil (pink; Figure 3C) component of their expression profile, it becomes clear that the total abundance of many potentially MYB-regulated transcripts does not change even though the pollen component does (blue dots along diagonal). Interestingly, this analysis also defines a small number of pistil-expressed genes that change significantly when the pistil is pollinated by mutant compared with wild-type pollen tubes (pink dots; Figure 3C).

By comparing the differentially expressed genes identified using total reads (ignoring parent-of-origin) with those determined following analysis of sex-specific reads, we find that 218 differentially expressed genes were unmasked when differential gene expression was determined for pollen and pistil reads independently (Figure 3C). This finding shows that SNP-informed RNA-seq is successful at defining changes in the pollen- and pistil-specific components of gene expression for those genes expressed in both tissues and indicates that many genes are under regulatory control in the pollen tube that is distinct from that of the pistil. However, since we used a different bioinformatic procedure to initially align total reads than that used to process reads with SNPs (see Methods), we repeated this analysis and compared sex-specific reads with all reads lacking a SNP; the result was similar to our previous analysis (Figure 3D). Many pollen tube-expressed genes were determined to be differentially expressed in myb triple mutants whether or not reads with SNPs were analyzed (Figure 3E, dark-blue intersection, n = 68); however, an additional 210 genes are identified by determining whether pollen reads (Col SNP) are differentially abundant (P < 0.05, Wald test, DESeq2) in the mutant compared with the wild type. Pollen tube-specific MYB-dependent genes comprise several functional annotation categories, the most enriched being transmembrane transporters (n = 26, P value = 4.46 × 10⁻⁵, Benjamini test; Supplemental Data Set 9).

By contrast, only a few pistil-expressed genes (34 genes; Figure 3E) were found to be differentially expressed in mutant compared with wild-type pollinations. Thirteen of these were also identified upon analysis of pistil (Cvi) reads. This small but interesting gene set is enriched in chitin-responsive genes (n = 6, P value = 9.08 × 10⁻⁵, Benjamini test) and may represent pistil genes that respond to myb mutant pollen tubes, which are perceived as different from the wild type (Supplemental Data Set 9).

To illustrate the consequence of MYB loss of function on the pollen tube component of gene expression, we plotted read coverage for two genes identified as differentially regulated in myb triple-mutant pollen tubes (Figure 3F; Supplemental Data Set 9). In wild-type pollinations (Figure 3F, top graph) the RNA-seq reads covering SNPs comprise both pistil (Cvi, red) and pollen (Col, blue) reads, with a 10-fold greater contribution from the pistil. In myb triple-mutant pollinations, the pollen component (Col, blue) disappears, yet the much larger pistil component is unchanged. This analysis shows the utility of SNP-informed RNA-seq analysis of reproduction, which facilitates independent analysis of the pollen tube and pistil components of gene expression as the pollen tube interacts with the pistil.

Genetic Reporter Constructs Corroborate SNP-Informed RNA-Seq Analysis

We produced transgenic plants expressing GFP fusion constructs to determine whether bioinformatic separation of pollen tube reads from pistil reads accurately identifies pollen tube-specific genes. We chose three MYB-dependent genes for which only the pollen (Col) allele was detected in reads over SNPs and were thus deemed pollen tube specific (Figures 4A, 4H, and 4O). We proposed that MYB expression, and thus target gene activation, is induced by growth through pistil tissue (Leydon et al., 2013, 2014). To test this hypothesis and to determine whether RNA-seq analysis was validated by reporter construct analysis, we analyzed GFP fusion proteins in the SIV pollen tube growth assay and in wild-type pistils.

• Download figure
• Open in new tab
• Download powerpoint

Figure 4.

Pollen-Specific Gene Expression of Three MYB-Dependent Genes Detected by SNP-Informed RNA-Seq.

(A), (H), and (O) Sequencing reads over SNPs identifying potential MYB target genes as pollen specific (as in Figure 2).

(B) to (F) SUC9:SUC9:GFP (AT5G06170) and LAT52:DsRed.

(I) to (M) PL or PL:PL:GFP (AT5G09280) and LAT52:DsRed.

(P) to (T) ACA4:ACA4:GFP (AT4G20990) and LAT52:DsRed.

(B), (I), and (P) Representative SIV pollen tube growth experiments on ms1 pistil explants (epifluorescent micrographs). In each, the left image is the red channel (LAT52:DsRed) signal, which is present in the pollen grain and the tube. The right image is the green channel (GFP) signal. The autofluorescence of the ms1 stigma explant is visible in each panel. Bars = 500 µm.

(C), (J), and (Q) Confocal micrographs of SIV-grown pollen tubes: left panel, GFP; middle panel, DsRED; right panel, merge. Bars = 20 µm.

(D), (K), and (R) Representative confocal optical sections of unpollinated pistils; the merge of differential interference contrast and GFP channels is shown. Arrows indicate expression of the GFP fusion protein in the female gametophyte. Bars = 100 µm.

(E), (L), and (S) Close up of unpollinated ovules from (D), (K), and (R). Bars = 100 µm, projection of red and green channels.

(F), (M), and (T) Projections of confocal micrographs of pollinated pistils. Red and green channels are merged, and arrows indicate regions of strong coexpression of DsRed and the indicated GFP fusion protein (yellow) in the style. Twelve hours after pollination. Bars = 100 µm.

(G), (N), and (U) Close-up of ovules targeted by GFP reporter constructs from (F), (M), and (T). Bars = 100 µm.

Supporting confocal micrographs are shown in Supplemental Figure 7. Supporting RT-qPCR analysis is presented in Supplemental Figure 8.

Full genomic loci containing the predicted promoter regions and coding sequences of SUCROSE-PROTON SYMPORTER9 (SUC9), PECTIN LYASE (PL), and ALPHA CARBONIC ANHYDRASE4 (ACA4) were cloned in frame with GFP to produce translational fusions (see Methods). We introduced these reporter constructs into A. thaliana (Col) carrying the pollen-specific LAT52:DsRed marker (Twell et al., 1989; Francis et al., 2007), which drives expression of soluble DsRED in the pollen tube cytoplasm. Mature pollen from each reporter was analyzed by confocal microscopy, but GFP was not detectable, consistent with the hypothesis that these MYB-dependent genes are induced by growth through pistil tissue. By contrast, after analyzing at least 21 independent primary transformants (T1) for each construct in SIV growth experiments, we found that pollen-tube expression was observed in the vast majority of transformants (SUC9:SUC9:GFP, 36/36; PL:PL:GFP, 18/21; ACA4:ACA4:GFP, 32/36; Figures 4B, 4I, and 4P).

We determined the subcellular localization of each GFP fusion protein in growing pollen tubes and each displayed a localization pattern consistent with the predicted localization. SUC9:GFP appeared to localize predominantly to the plasma membrane and to vesicles in the pollen tube tip (Figure 4C; Supplemental Movie 1). PL:GFP was observed in the endomembrane system (puncta) and what is likely the pollen tube cell wall (Figure 4J; Supplemental Movie 2). ACA4:GFP was observed in the endomembrane system (puncta, Figure 4Q; Supplemental Movie 3). In all SIV experiments where GFP fluorescence was observed (86/93), expression was apparent in extending pollen tubes and not pollen grains.

To determine whether these genes are pollen tube specific, as they had been designated by SNP-informed RNA-seq, we examined three independent transformants for the expression of GFP unpollinated pistils using confocal microscopy. For each transgenic line, we observed only autofluorescence in the stigma, style, transmitting tract, and all other somatic tissues of the pistil (Figures 4D, 4K, and 4R). Surprisingly, we observed expression of SUC9:GFP in the synergid cells (Figures 4D and 4E; Supplemental Figures 7A, 7B, and 7I) and ACA4:GFP in all cells of the female gametophyte (Figures 4R and 4S; Supplemental Figures 7C, 7D, and 7I); as expected, PL:GFP was undetectable in any pistil tissue. We then analyzed self-pollinated pistils for each reporter construct (Figures 4F, 4M, and 4T). Each reporter was most strongly detected in pollen tubes, especially in the style region where pollen tubes are concentrated (arrows). Additionally, each reporter fusion was detected in ovules clearly targeted by a pollen tube (Figures 4G, 4N, and 4U), supporting a possible role for these MYB-dependent genes in pollen tube reception (Leydon et al., 2013).

Expression of SUC9:SUC9:GFP and ACA4:ACA4:GFP reporter constructs in the female gametophyte was surprising in light of our SNP-based RNA-seq results, where we detected no SUC9 or ACA4 Cvi SNP-containing transcripts. Transgenic lines were generated in the Col accession, whereas our RNA-seq used Cvi as the pistil donor. So, we tested whether expression of SUC9, PL, and ACA4 in the unpollinated pistil was accession-specific using RT-qPCR. Of the three genes tested, only SUC9 transcripts were detected at low levels in Col unpollinated pistils; none were detected in either Cvi or Ler (ms1) unpollinated pistils (Supplemental Figure 8). This observation is consistent with SUC9:SUC9:GFP expression being detected in synergid cells of transgenic Col plants, but not in Cvi, the female used in our RNA-seq experiments. Interestingly, ACA4 transcripts were undetectable in unpollinated pistils of all three accessions, suggesting that gametophytic expression of ACA4:ACA4:GFP may be due to inclusion or exclusion of cis-regulatory regions in the promoter region of this transgene. Pectate lyase transcripts were also undetectable in unpollinated pistils, consistent with our RNA-seq results.

The observation that SUC9 was expressed in the synergid cells of Col but not other accessions (Ler and Cvi), prompted us to test whether accession affected SUC9, PL, and ACA4 expression in self- and cross-pollinations. Col self-pollinations result in higher levels of each transcript compared with Cvi (Supplemental Figure 8). All three genes behaved similarly in Col and Ws self-pollinations, except for SUC9, which was induced at lower levels in Ws. Furthermore, Cvi x Col inter-accession pollinations induce SUC9, PL, and ACA4 transcripts at lower levels than Col self-pollinations, suggesting that pollen-pistil interactions, not pollen genotype alone, affect the strength of MYB-dependent pollen tube transcriptional activation in vivo (Supplemental Figure 8). This accession-specific activation of MYB-dependent genes led us to examine the role of accession in reproductive gene expression in our RNA-seq data sets.

Gene Expression Changes Are Induced by Inter-Accession Pollination

Col and Cvi are interfertile; however, inter-accession crosses show subtle defects in pollen tube reception that resemble pollen tube overgrowth defects in fer, lre, or myb triple-mutant pollinations (Escobar-Restrepo et al., 2007; Capron et al., 2008; Leydon et al., 2014; Müller et al., 2016). These observations along with the finding that pollen tube-expressed genes are differentially expressed between accessions prompted us to determine whether there are genome-wide, accession-specific changes in pollen and pistil gene expression. First, we used hierarchical clustering to assess whether the transcriptomes of inter-accession pollinations were distinct from intra-accession pollinations or unpollinated pistils. As expected, hierarchical clustering demonstrates that pollinated pistils are distinct from unpollinated Cvi pistils (Figure 5A). We compared the transcriptomes of self-pollinations (Cvi x Cvi and Col x Col) to inter-accession pollinations (Cvi x Col) and found a core set of 16,838 genes expressed at detectable levels in each pollination type (>100 reads; Figure 5B). However, differential expression analysis identified hundreds of genes with significantly different levels of expression between Col and Cvi self-pollinations (Figure 5C; Supplemental Data Set 10).

• Download figure
• Open in new tab
• Download powerpoint

Figure 5.

Reproductive Gene Regulation Is Highly Divergent between Col and Cvi.

(A) Hierarchical clustering of all RNA-seq samples identifies accession-specific clusters (DESeq2).

(B) Intersection of all expressed genes (mean expression >100 reads) from self and inter-accession crosses. Venn diagram is not proportional to number of included elements.

(C) Differential gene expression observed between Col and Cvi 8 h after self-pollination. Differentially expressed genes are colored purple (adjusted P value < 0.005, DESeq2).

(D) to (F) Heat maps of Col-specific gene expression; colors represent read volume, with scales below indicating the number of reads for color value.

(D) Defensin-like genes.

(E) Low-molecular-weight cysteine-rich genes.

(F) Toll-interleukin receptor-nucleotide binding site-leucine rich repeat (TIR-NBS-LRR) type disease resistance proteins.

(G) Differential gene expression between self-pollinated Col and Cvi pistils. A total of 2877 differentially expressed were identified following pollination; 1640 transcripts were higher in inter-accession (Cvi x Col) pollinations, and 1238 transcripts were higher in self-pollinations (Cvi x Cvi). Adjusted P value < 0.005 (DESeq2).

(H) Sex-specific transcript identity of (A). Red dots, Cvi-specific reads; blue dots, Col-specific reads; gray dots, no SNPs for assignment.

Accession-specific transcriptional responses have been previously observed (Autran et al., 2011; Nodine and Bartel, 2012; Del Toro-De León et al., 2014), so it is not surprising that Col and Cvi accessions display variation in transcriptional profiles of reproductive cells and tissues. Interestingly, we identified members of gene families that were enriched in the Col pollination transcriptome that were absent from Cvi pollinations such as small secreted defensin-like proteins (Figure 5D), low molecular weight cysteine-rich proteins (Figure 5E), and NBS-LRR genes (TIR-NBS-LRR; Figure 5F). By contrast, analysis of genes enriched in Cvi pollinations, but absent in Col, did not show enrichment for any particular GO enrichment or gene families (see Methods). However, despite the lack of enriched GO terms suggesting roles in extracellular signaling, further analysis revealed that 44% of Cvi-specific pollination genes were predicted to contain signal peptides, similar to ∼31% found in Col-specific pollinations (SignalP). These data suggest that Cvi may express its own repertoire of secreted proteins during pollination, but that unlike Col, the gene families being expressed are not as well annotated.

The Pistil Transcriptome Responds to the Pollen Genotype during Pollination

We set out to determine whether the pistil transcriptome changes in response to pollen accession in the pollinated pistil. Our analysis of myb triple-mutant pollinations suggested the pistil responds differently to mutant versus wild-type pollen tubes (Figures 3C and 3E), and we reasoned the pistil may also mount unique responses to different pollen accessions. We defined pollen tube and pistil changes in gene expression (differentially expressed genes) that occur during intra-accession (Cvi x Cvi) versus inter-accession (Cvi x Col) pollinations (Figure 5F, purple dots). First, we ascertained genes that were differentially expressed between self-crosses (Cvi x Cvi) and inter-accession crosses (Cvi x Col; Figure 5G; Supplemental Data Set 10). This analysis identified 2877 differentially expressed genes; 1238 are significantly more abundant in intra-ecotypic pollinations, 1640 are significantly higher in inter-ecotypic pollinations (632 up in inter-accession pollinations and 418 up in intra-accession pollinations). We then determined the parent of origin for these genes and identified 1050 pistil-expressed genes that are differentially expressed depending on the pollen genotype (Figure 5H, red dots; Supplemental Data Set 10). Most of the accession-sensitive transcriptional changes were subtle (<4-fold). The 632 pistil genes that were differentially upregulated in inter-accession pollinations were enriched for responses to various molecules (22.38 enrichment score for this cluster of responses), including chitin (n = 37, P = 2.6 × 10⁻²⁵), carbohydrates (n = 39, P = 2.26 × 10⁻²⁰), and organic substances (n = 83, P = 1.2 × 10⁻¹⁵). Transcription factors were also enriched (n = 61, P = 1.6 × 10⁻⁶). Interestingly, genes induced in the inter-accession (Cvi x Col) pistil, included 14 ethylene-responsive transcription factors that are implicated in pathogenesis-related transcriptional response (n = 14, P = 4.7 × 10⁻³, Supplemental Data Set 11). The 418 pistil-expressed genes that were enriched in intra-accession crosses (Cvi x Cvi) were not strongly enriched for these or other functional categories. This subtle difference in transcriptional profiles indicates that pollen accession is perceived by the pistil, activating genes associated with defense responses during inter-accession pollinations.

Expression of a Large Family of Thionin-Like CRPs Is MYB Dependent and Accession Specific

Previous analysis of myb triple-mutant pollinations using microarray analysis had identified a single thionin-like gene (AT5G36720) as being among the most differentially regulated between mutant and wild-type pollinations (Leydon et al., 2013). This prompted us to propose that these proteins may be secreted by pollen tubes to female gametophytes to mediate pollen tube reception signaling. However, analysis of this family (cysteine-rich protein 2460 [CRP2460]; Silverstein et al., 2007) indicated that many members of this large group of genes were not annotated and therefore were not represented on the ATH1 array. We analyzed this family using our RNA-seq data to determine whether multiple members were potentially MYB regulated and to determine whether their expression varied among accessions. We plotted the RNA expression of each family member alongside the CRP2460 phylogenetic tree (Figure 6; Supplemental Data Set 12). Differential expression analysis identifies 22/28 CRP2460 transcript levels as statistically lower in pollinations with myb triple-mutant pollen tubes (Figure 6B, asterisk) compared with the wild type. CRP2460 thionin proteins were expressed accession preferentially: High expression was observed in Col pollen tubes, while Cvi pollen tubes expressed fewer members and at lower levels (Figure 6B). Of the 28 thionin-like proteins analyzed, 15 contained polymorphisms that allowed us to determine whether expression of these genes was from the pollen tube or pistil (Figure 6B, +). In each case, all reads over SNPs carried the Col allele, indicating that thionins are pollen tube expressed (Figure 6C). Thionin proteins therefore represent a pollen tube expression signature that differentiates A. thaliana accessions from each other and may be involved in signaling between the pollen tube and the female gametophyte important for sperm release.

• Download figure
• Open in new tab
• Download powerpoint

Figure 6.

MYB Transcription Factors Control Expression of Thionin CRPs (CRP2460).

(A) A phylogenetic tree of the CRP2460 protein coding sequences (neighbor-joining, jukes-cantor). Two MYB-dependent members of other CRP groups (CRP2465 and CRP2530) are included for comparison.

(B) Expression heat map of the CRP2460 transcripts in pistils. Cvi Unpoll, unpollinated Cvi pistils (collected 8 h after time when pollination would have been done); Cvi x Col, Cvi pistils pollinated with Col pollen for 8 h; Cvi x myb, Cvi pistils pollinated with myb triple-mutant pollen for 8 h; Col x Col, Col pistils pollinated with Col pollen for 8; Cvi x Cvi, Cvi pistils pollinated with Cvi pollen for 8 h. Sequencing data are shown in three columns per condition, with one column per biological replicate, and color intensity is based on read number (see scale). Asterisks indicate statistically significant gene expression differences between Cvi x Col and Cvi x myb triple-mutant pollinations (adjusted P value < 0.005, DESeq2). Green plus sign indicates the presence of SNPs in these CRP genes. Ka/Ks, numbers indicate Ka/Ks ratio between CRP and homolog in A. lyrata. Gray box indicates no homolog present for this CRP.

(C) Representative sequencing data for two MYB-dependent CRP2460 genes, AT1G58242 and AT3G24465.

MYB-Dependent, Pollen Tube-Expressed Genes Show Elevated Polymorphism Rates, Gene Family Expansion, and Positive Selection

myb triple-mutant pollen tubes fail to release sperm and continue to grow within the female gametophyte (Leydon et al., 2013; Liang et al., 2013). This defect, first identified in the A. thaliana feronia mutant (Huck et al., 2003; Rotman et al., 2003), had been characterized as a prezygotic barrier to hybridization in interspecific cross-pollinations (Williams et al., 1986; Escobar-Restrepo et al., 2007; Müller et al., 2016) and suggests that MYB-dependent genes may be involved in shaping the species identity of the pollen tube that is recognized by the female gametophyte (Leydon et al., 2014). Molecules involved in gamete:gamete recognition have been shown to evolve rapidly and to be responsible for species-level discrimination in mating interactions (Metz et al., 1998; Clark et al., 2006). Interestingly, it was recently shown that the pollen tube expresses a larger fraction of species-limited genes than other tissues (Arunkumar et al., 2013; Gossmann et al., 2014), which is consistent with the idea that pollen tubes may express species-limited proteins that determine their species-level identity. To address whether groups of genes expressed in reproductive tissues maintain higher rates of sequence diversity, we calculated Col versus Cvi exonic SNP rates and compared them to control rates simulated using 100 sets of randomly selected genes of the same number. We used the ratio of the measured SNP rate to the simulated rate to identify groups of genes with elevated SNP rates (Table 1; Supplemental Data Set 13).

Analysis of groups of genes defined by their expression pattern during pollination without determining parent of origin (analysis of all reads) revealed that the ∼10,000 genes expressed in the pollinated pistil have SNP rates that are significantly below the simulated rate. However, genes induced by pollination contain significantly higher than control rates (Table 1). The 292 genes that are differentially expressed in myb triple-mutant pollinations (identified from all mapped reads) do not show elevated SNP rates. When analyzed in bulk, we only observe increased polymorphism rates in genes induced by pollination.

However, when sex-specific reads are analyzed to focus on changes occurring in pollen and/or pistil gene expression, three groups were found to have polymorphism rates elevated compared with the simulated rate: pollen-expressed genes, pistil-specific pollination-induced genes, and MYB-dependent genes (Table 1). To gain insight into MYB-dependent genes with elevated polymorphism rates, we focused on a MYB-dependent CRP gene family (Figure 6). We started by analyzing the NBS-LRR gene family, which was previously found to have high rates of major effect SNPs in A. thaliana accessions indicative of diversifying selection (Clark et al., 2007; Schmitz et al., 2013) and found a ratio of measured SNP rate to control rate (1.24; Table 1) that was significantly above one. A significantly high SNP rate (P < 0.01, Poisson distribution) was also observed for a group of 825 CRP encoding genes, but not for the thionin subgroup (n = 71). We also determined the frequency that SNPs in these groups of genes resulted in missense mutations. Interestingly, the thionin subclade of CRPs had the highest rate of missense mutations (Table 1; Supplemental Data Set 14), suggesting this group of genes may be under diversifying selective pressure.

To test for positive selection within thionins (CRP2200-2610), we first determined whether there were homologs for each gene in Arabidopsis lyrata, Populus trichocarpa, grape (Vitis vinifera), rice (Oryza sativa), and Physcomitrella patens using a recently published genome-wide analysis (Supplemental Data Set 15; Lloyd et al., 2015). Appreciable fractions of Actin and NBS-LRR gene families have homologs in distantly related species such as the monocot, rice, and the moss, P. patens (Figure 7A). By contrast, homologs were only identified in closely related A. lyrata for the vast majority of CRP2200-2610 genes (Figure 7A). For example, 30% of NBS-LRR genes have homologs in P. patens, whereas no homologs for CRP2200-2610 genes were identified (Lloyd et al., 2015). Moreover, only a small subset of the CRP2200-2610 genes have homologs in angiosperm species outside the Brassicaceae (P. trichocarpa, grape, and rice), suggesting a recent expansion of this group.

• Download figure
• Open in new tab
• Download powerpoint

Figure 7.

Elevated Ka/Ks Rates in CRP Genes.

(A) Percentage of genes in indicated gene families with homologs identified in given species identified by Lloyd et al. (2015).

(B) Average Ka/Ks of all homologs found between A. thaliana and the indicated species. Bars are sd (sample sizes vary and are defined in [A]; Lloyd et al., 2015).

(C) Percentage of homologous gene pairs out of the total number of family members with a Ka/Ks >1 (indicator of positive selection) for each respective species. The key pertains to (A) to (C).

To profile diversification in these families, we plotted the mean Ka/Ks for genes with homologs; we found that Actin genes had lower rates than NBS-LRRs, consistent with the hypothesis that NBS-LRR genes are under diversifying selection, while the Actin genes are under purifying selection (Figure 7B). CRPs exhibit a Ka/Ks profile that is similar to NBS-LRR genes (Figure 7B). However, this analysis of mean Ka/Ks failed to highlight the elevated proportion of CRP gene family members (CRPs in general and CRP2200-2610 in particular) with a Ka/Ks above one in the comparison with closely related A. lyrata, which had the largest number of homologs (Figure 7C).

This analysis revealed a dramatic expansion of the CRP2200-2610 thionin clade of CRPs in the Brassica lineage. Of particular interest is the CRP2460 (Figure 6) group, which has undergone an expansion in A. thaliana compared with the closely related A. lyrata. Only six of 28 genes within this group have homologs in A. lyrata, and none have homologs in any of the other species tested (Figure 6; Ka/Ks indicated). Furthermore, the Ka/Ks rate for At1g21835 compared with its homolog from A. lyrata was found to be above one (Ka/Ks = 1.06; Figure 6), suggesting selective pressure to maintain diversity of this gene. At1g21835 appears to be the founder of an expanded group of A. thaliana thionins because it is sister to 14 genes without homologs in A. lyrata (Figure 6; Supplemental Figure 9). Analysis of this expanded group revealed evidence for diversifying selection at specific sites within these ∼100-amino acid proteins (Supplemental Figure 9 and Supplemental Data Set 15). These results are in keeping with our proposal that this group of secreted proteins convey species-unique information to the female gametophyte required for pollen tube reception.