RNA-seq of Arabidopsis pollen uncovers novel transcription and alternative splicing

Abstract

Pollen grains of Arabidopsis (Arabidopsis thaliana) contain two haploid sperm cells enclosed in a haploid vegetative cell. Upon germination, the vegetative cell extrudes a pollen tube that carries the sperm to an ovule for fertilization. Knowing the identity, relative abundance, and splicing patterns of pollen transcripts will improve our understanding of pollen and allow investigation of tissue-specific splicing in plants. Most Arabidopsis pollen transcriptome studies have used the ATH1 microarray, which does not assay splice variants and lacks specific probe sets for many genes. To investigate the pollen transcriptome, we performed high-throughput sequencing (RNA-Seq) of Arabidopsis pollen and seedlings for comparison. Gene expression was more diverse in seedling, and genes involved in cell wall biogenesis were highly expressed in pollen. RNA-Seq detected at least 4,172 protein-coding genes expressed in pollen, including 289 assayed only by nonspecific probe sets. Additional exons and previously unannotated 5' and 3' untranslated regions for pollen-expressed genes were revealed. We detected regions in the genome not previously annotated as expressed; 14 were tested and 12 were confirmed by polymerase chain reaction. Gapped read alignments revealed 1,908 high-confidence new splicing events supported by 10 or more spliced read alignments. Alternative splicing patterns in pollen and seedling were highly correlated. For most alternatively spliced genes, the ratio of variants in pollen and seedling was similar, except for some encoding proteins involved in RNA splicing. This study highlights the robustness of splicing patterns in plants and the importance of ongoing annotation and visualization of RNA-Seq data using interactive tools such as Integrated Genome Browser.

Figure 1.

Concordance between RNA-Seq and microarray present calls. Each diagram reports the size of different subsets of single-target (no cross hybridization) ATH1 probe set/target pairs. The array subset contains probe sets that were called as present in at least five of six mature, dry pollen ATH1 microarrays from Gene Expression Omnibus expression series GSE6696 and GSE17343, where the sample collection protocol was the same as that used for the RNA-Seq library. The RNA-Seq subset contains probe set target genes with pollen RPM expression values as high or higher than the indicated threshold. The intersections of RNA-Seq and array subsets are genes/probe sets called present by both methods. Values reported in the lower right corner of each diagram are probe set/targets called absent by both methods.

Figure 2.

Reads mapping outside any known gene or transposable element. A, Distribution of log-transformed (base 10) intergenic region (blue) and annotated gene sizes (red), in bp. Values on the y axis indicate the distribution density. B, Heat maps indicate the density of reads mapping to the indicated chromosome regions of the Arabidopsis genome. Darker (blue) hues indicate greater coverage, and yellow hues indicate less coverage. Individual gene models are shown as black marks above each chromosome; base pair positions are indicated. In most cases, stacked gene models that occupy the same horizontal location are annotated transcript variants arising from alternative splicing, alternative promoters, or alternative 3′ end processing.

Figure 3.

Relationships between number of reads and region size. The x axis shows the number of reads per region, and the y axis shows the size per region in bp. Both axes are on logarithmic (base 10) scale. There were 35 regions 500 bases or larger with 150 or more reads.

Figure 4.

Putative new genes expressed in pollen. Screen captures from IGB show regions of novel transcription outside previously annotated loci, including annotated transposable elements. Read alignments are in the top track. Thin lines indicate reads mapping across introns. Graphs indicate the number of reads that map at positions indicated in the coordinates track. The bottom track contains TAIR10 gene model annotations, including protein-coding genes, transposable element genes, noncoding genes, and pseudogenes. The chromosome number is indicated in the upper right corner of each image. Regions are as follows: chr2:3198768 to 3200220 (A), chr2:7632907 to 7634053 (B), chr1:5355736 to 5356794 (C), chr2:9256704 to 9258245 (D), chr5:23406941 to 23407752 (E), chr4:7262280 to 7263282 (F), chr1:4339027 to 4339918 (G), and chr5:9461451 to 9462429 (H). [See online article for color version of this figure.]

Figure 5.

Visualization from IGB depicting additional exons for RIC5. Single-mapping reads from pollen and seedling RNA-Seq data sets are shown aligned onto the reference Arabidopsis genome. The top row of the pollen reads track contains reads that are drawn on top of each other due to space limitations. [See online article for color version of this figure.]

Figure 6.

Visualizations from IGB depicting novel junction features. The number of spliced reads supporting junctions predicted from spliced read alignments appear above each junction feature. The annotated gene model is shown in black above the sequence axis. Arrows indicate the direction of transcription. [See online article for color version of this figure.]

Figure 7.

Splice pattern correlation between ESTs, pollen, and seedling RNA-Seq data sets. The percentages of ESTs or RNA-Seq reads supporting the shorter form (G_A) variant for all annotated alternative splicing events in TAIR10 are shown as scatterplots. Events with 20 or more overlapping reads are shown.

Figure 8.

Pollen-specific alternative splicing in U2AF65A. A, Three gene models for U2AF65A. Blocks represent exons, and taller blocks indicate translated regions within an exon. The gene is transcribed from the minus strand of chromosome 4, and transcription is shown as proceeding from right to left. Not shown is variant 4, which is identical to variant 1 in the variably spliced 3′ region. Arrows indicate PCR primer locations. B, Closeup of the 3′ differentially spliced region with scored junction features from pollen and seedling RNA-Seq data sets, where scores represent the number of spliced reads supporting the junction. C, PCR products from amplifying pollen and seedling cDNA samples. Pollen samples included the same pollen RNA used to create the pollen library (P09) and RNA extracted from two independent pollen collections (P11 and PC). S1 and S2 are seedling cDNA samples prepared independently from the RNA-Seq seedling libraries. D, Relative amounts of variants 1, 2, and 3 PCR product amplified from pollen and seedling cDNA. Relative amounts were calculated from C. [See online article for color version of this figure.]