Analysis of pollen-specific alternative splicing in Arabidopsis thaliana via semi-quantitative PCR

Bioinformatic analysis of splicing patterns

Our previously published report described using RNA-Seq alignments to identify differential alternative splicing in pollen and leaf RNA-Seq data sets (Loraine et al., 2013). Briefly, alternative splicing (AS) events annotated as part of the TAIR10 gene models from Arabidopsis thaliana were identified by comparing regions of overlap between exons and introns in different gene models belonging to the same gene, as described in English, Patel & Loraine (2010). For each alternative splicing event, the number of spliced reads supporting mutually exclusive choices was identified for each of two leaf libraries and one pollen library. For simplicity, only annotated alternative donor and acceptor sites were considered, and alternative splice forms were compared as pairs. For each AS event, the percentage of gapped reads in each library that supported each form was calculated and compared between pollen and leaf samples. Alternative splicing events were identified in which percentage of reads supporting a given splicing choice varied more than 20% between the two sample types. The AS events were then ordered by this percent difference; the listing of these events is shown in Data S1. The ten most differentially spliced genes were identified for testing by PCR; the list included At-U2AF65A, which was tested previously and is not discussed again here (Loraine et al., 2013). Note also that previously we referred to the leaf samples as “seedlings,” but here refer to them as “leaf” samples. This is because the older terminology was confusing to some readers. For many people the term “seedling” means: whole seedlings that are two weeks old or younger. The samples used in here and in the previous study were from three-week old, soil-grown plants and included leaves and not roots.

PCR primer design for analysis of splicing patterns

Primers flanking the differentially spliced region were designed using the Primer3 software deployed on the NCBI Web site. Primer locations within the target gene were visualized using Integrated Genome Browser (Nicol et al., 2009) and amplicon sizes were calculated for each splice variant. Where possible, primers were designed to have predicted melting temperatures of 60 °C.

Pollen collection

Pollen RNA samples were collected and prepared as described in Johnson-Brousseau & McCormick (2004) and Loraine et al. (2013). Pollen samples corresponding to samples P1, P2, and P3 were collected by vacuuming flowers from different sets of pots growing under 16 h/8 h light/dark illumination. P1 and P2 were collected in 2012 and P3 was collected in 2008. P3 was also the same sample used to synthesize pollen cDNA that was used for sequencing (Loraine et al., 2013). Pollen samples were frozen following collection and stored at −80 °C prior to use.

Arabidopsis growth conditions

Leaf samples used in this study corresponded to non-stressed controls from pot “L” described in Gulledge et al. (2012). The above-ground parts of plants (leaves but not roots or stems) of plants growing on soil were collected after 21 days of growth under 16/8 light/dark illumination in an incubator (model AR36L3 from Percival Scientific, Perry, Iowa, USA) set to 22 °C, 45% relative humidity, and 150 photosynthetic photon flux density (PPFD). The plants had not yet bolted at the time of collection.

RNA extraction

Frozen pollen or leaf samples were ground in liquid nitrogen with a mortar and pestle and RNA was extracted using a RNeasy kit (Bio-Rad catalog number 74904: Bio-Rad; Hercules, California, USA) following manufacturer’s instructions. RNAs were DNaseI treated via an on-column digestion with the Plant RNeasy Kit. RNA quality was assessed with an Agilent BioAnalyzer and only high quality samples were used.

cDNA synthesis

cDNA samples were synthesized following the SMART MMLV reverse transcriptase kit protocol (Company is Clonetech, catalog number 639537; Clontech, Mountain View, California, USA), using the supplied 5X buffer and DTT supplied with the kit, 50 mM oligo dT primers (Bioline, catalog number BIO-38029; Bioline, Taunton, Massachusetts, USA) and 10 mM dNTPs (Bioline, cat. no. 39053; Bioline, Taunton, Massachusetts, USA). 1 µg of total RNA was used in 20 µL cDNA synthesis reactions.

PCR amplification of cDNAs from pollen and leaf cDNAs

PCR amplifications included an initial “hot start” 4 min denaturation incubation at 94 °C followed by 30 cycles of denaturation (94 °C), annealing, and elongation (72 °C) incubations. Annealing temperatures were selected based on Tm predicted by the Primer 3 software, and when primers failed to amplify a band, gradient PCR was performed to identify the optimal annealing temperature. Amplifications were typically performed in volumes of 20 µL per template/primer combination with 2 µL each of 5 µM forward and reverse primer stock solutions, .08 µL of taq polymerase (Bioline, catalog number 21086; Bioline, Taunton, Massachusetts, USA), and 0.5 µL aliquots of pollen or leaf cDNA reactions that had been diluted 1:20. Cycle numbers were estimated to be in logarithmic phase based on prior analysis of qPCR analysis of genes of varying starting RPKMs (Loraine et al., 2013).

Analysis of gel images for semi-quantitative PCR

The relative abundance of splice variants amplified by PCR was estimated using methods presented in Palusa, Ali & Reddy (2007) and Venables et al. (2012). Equal amounts of PCR reaction mixtures were fractionated by agarose gel electrophoresis containing ethidium bromide. Gels were visualized under short wave UV light in a UVP Biospectrum AC Imaging System GelDoc station and photographed using a digital camera. Each gel was typically photographed at different run times and using different exposure times to obtain high quality images for analysis.

Gel images were analyzed using the VisionWorksLS version 7.0.1 software program. Rectangular regions surrounding ethidium-stained bands in each lane were identified, and for each band, the software was used to estimate molecular weight along with a background-corrected band intensity value, called “I-vol” in the software. Visualization of lane profiles was also done to ensure that band intensity measurements were within the range of detection of the instrument. Gels also included lanes containing molecular mass standards (Hyperladder 1 kb, Bioline catalog number BIO-33053; Bioline, Taunton, Massachusetts, USA). Band intensities of the molecular size/mass markers were plotted against their masses and only gels where this relationship was close to linear (R² > = 0.90) were used. Band intensities and size estimates were then imported into the R statistical programming environment for analysis.

For each band, its abundance relative to other bands in the same lane was calculated and used to test for differences between pollen and leaf samples. The relative abundance of per band was calculated as its percent contribution the total amount of DNA in each lane. The percent contribution per band was calculated as 100 × b/T, where b was the band’s size-normalized intensity and T was the sum of all size-normalized band intensities in the lane, calculated as $\sum _{i=1}^N{b}_i$. Size-based normalization of band intensity was done to account for how ethidium staining (and band brightness) depended on both amplicon abundance and amplicon size, since larger amplicons contained more sites (base pairs) for ethidium to bind (Green & Sambrook, 2012). For each gel containing replicate pollen and leaf samples, a t test was used to determine whether there was a statistically significant difference in a band’s percentage between pollen and leaf samples. As the prior RNA-Seq data analysis indicated whether pollen or leaf percentages were greater and since there was no reason to expect unequal variance, one-sided t-tests assuming equal variance were used. R code used to calculate percentages and perform statistical testing are available as Markdown (.Rmd) in the folder named “GelAnalysis” in the project repository at http://www.bitbucket.org/lorainelab/pollenas. For each gene tested, a Markdown file and its corresponding HTML report file is available; readers can re-run the code or open the HTML file in a Web browser to view the code, the output of the code, and accompanying figures.

Sanger sequencing of PCR products

For higher-abundance amplicons that did not match annotated gene models or that suggested use of alternative transcription start sites, PCR products were extracted from the gels and sequenced by conventional Sanger sequencing to confirm band identity. Sequence reads obtained from Sanger sequencing were aligned on the Arabidopsis genome using blat (Kent, 2002) and visualized in Integrated Genome Browser for comparison with gene models. Blat was run with default settings except for the maxIntron parameter, which was set to 2,000, since 99% of Arabidopsis introns are 2,000 bases or smaller. Sequencing results are described in detail in the folder named “AmpliconSequencing” in the project repository at http://www.bitbucket.org/lorainelab/pollenas.

Data access

The RNA-Seq data are available from the Sequence Read Archive under SRA accessions SRR847501, SRR847502 (pollen library); SRR847503, SRR847504 (leaf library 1); and SRR847505, SRR847506 (leaf library 2). Read and junction feature alignments depicted in figures are available from the IGBQuickLoad data source named “pollen” under the Arabidopsis June 2009 (TAIR9/10) genome assembly release. A tutorial introducing Integrated Genome Browser and describing how to perform a visual splicing analysis using the pollen RNA-Seq data set is available at http://tinyurl.com/pk25nbm.

Splice variants from Arabidopsis SR genes At-SR30 (AT1G09140) and At-RS41 (AT5G52040) are differentially spliced in pollen and leaves

At-SR30 and At-RS41 encode two of eighteen serine/arginine-rich (SR) RNA-binding proteins in Arabidopsis (Barta, Kalyna & Reddy, 2010) thought to play a role in regulation of splice site choice. According to the RNA-Seq data, At-SR30 (Fig. 1A) and At-RS41 (Fig. 2A) were differentially spliced in pollen and leaf. Figures 1A and 2A show junction features aligned onto each gene, where features summarize the number of spliced reads in each sample that supported the corresponding intron in the gene models. In both genes, the relative abundance of spliced reads supporting alternative splicing events was clearly different in leaf versus pollen. We also observed a small number of unspliced reads that aligned to the differentially spliced intron, indicating a low level of intron retention (not shown), a common form of alternative splicing in plants (English, Patel & Loraine, 2010).

To validate the RNA-Seq results, we designed PCR primers flanking the differentially spliced region and performed semi-quantitative PCR using newly prepared pollen and leaf cDNA samples. For both genes, semi-quantitative PCR confirmed the differential splicing patterns found by RNA-Seq. Figures 1B–1C shows that two splice variants of At-SR30 were expressed in both pollen and leaf samples, but the relative abundance of the two variants was different. We also observed a low abundance, higher molecular weight amplicon in pollen and leaf samples consistent with intron retention. For At-RS41, we observed similar results; PCR amplified four splice variants, the relative abundance of the splice variants was different between pollen and leaves, and low abundance amplicons consistent with intron retention were apparent (Figs. 2B–2C).

In both genes, differential splicing affected the protein-coding region of the transcripts, but in different ways. In At-SR30, the leaf-preferred isoform introduced a premature stop codon that deleted 23 amino acids from the carboxyl end of the protein, including five SR repeats, which are possible phosphorylation targets and mediate interactions with other proteins (reviewed in Shepard & Hertel, 2009). For At-RS41, the pollen-preferred isoform included an alternative exon that introduced a new translation start site leading to loss of 35 N-terminal amino acids encoding part of an RNA recognition motif. Thus, the pollen-preferred form of At-SR30 encoded the full-length protein, while the pollen-preferred form of At-RS41 encoded a truncated form lacking part of the amino terminal region of the protein.

Alternative splicing in pollen remodels starch synthase gene AT5G65685

AT5G65685, encoding a starch synthase-like protein, is highly alternatively spliced in the 5′ region (Fig. S1A). Visualization of RNA-Seq alignments suggested there were many possible splice variants affecting the amino terminus of the protein and that the splicing patterns in pollen and leaves were different. According to the RNA-Seq data, the first intron (relative to isoform AT5G65685.1) contained two alternative donor sites that were differentially spliced in pollen compared to leaves, with leaf samples using the upstream site and pollen samples using the downstream site. The RNA-Seq data also suggested there were several unannotated variants containing retained introns. Results from semi-quantitative PCR were consistent with these observations (Figs. S1B–S1C). PCR primers designed around the alternatively spliced region amplified multiple bands corresponding to annotated variants and to many novel, lower abundance variants corresponding to retained intron isoforms.

Alternative splicing of AT4G21720, gene of unknown function

AT4G21720 encodes a conserved protein of unknown function with two annotated splice variants, with alternative acceptor sites affecting the 3′ end of the second intron (Fig. S2A). Visualization of RNA-Seq read alignments in Integrated Genome Browser indicated that the upstream acceptor site, which introduced a frameshift that eliminated 44 amino acids, was preferred in pollen, and the downstream acceptor site, which allowed translation of the full-length protein, was preferred in leaves.

As before, primers flanking the alternatively spliced region were designed and used to amplify cDNAs prepared from pollen and leaves (Figs. S2B–S2C). The PCR assay results supported differential splicing of AT4G21720 in pollen versus leaves. PCR products consistent with both acceptor sites were observed and the relative abundance of the products was similar to what was observed in the RNA-Seq data (Figs. S2B–S2C and Data S1). An additional higher molecular weight band was also observed in both samples, and its size suggested it was derived from occasional retention of the alternatively splice intron. However, its abundance was extremely low.

5′ alternative splicing coincided with alternative TSSs in four of nine genes tested

AT3G17120, another gene of unknown function, is annotated as producing two splice variants due to exon skipping in the 5′ UTR of the gene, and the RNA-Seq data analysis indicated that the relative abundance of splice variants was different in pollen and leaves (Data S1). However, visual inspection of RNA-Seq alignments alongside the gene models suggested the gene contained an alternative transcription start site (TSS) located within the first annotated intron (Fig. 3A). We tested this in PCR experiments using variant specific primers F1 and F2 and a common primer R. In leaves, F1 with R amplified a product, but F2 and R did not (Fig. 3B). In pollen, F2 and R amplified a product, but F1 with R did not. Combining all three primers in the same reaction amplified leaf and pollen-specific products, and sequencing confirmed the identity of products shown in Fig. 3B. This confirmed that pollen and leaves expressed distinct isoforms and strongly suggested that the alternative transcripts arose from alternative transcriptional start sites. Because the suggested pollen alternative TSS was downstream of the 5′ end of the first intron, transcripts arising from this alternative TSS would therefore require a different 5′ donor site to be used in splicing. In this case, alternative transcription and alternative splicing were linked.

Alternative splicing of AT3G06330, encoding a putative RING-type ubiquitin ligase, also appeared to coincide with alternative TSSs. AT3G06330 is annotated as having an alternative donor site in the 5′ end of the gene that changes the amino terminal region of the protein. According to the RNA-Seq data, the downstream donor site was preferred in pollen, while the upstream donor site was preferred in leaves (Fig. S3A). As before, we designed variant-specific forward primers and used them together with a common reverse primer in PCR reactions with pollen and leaf cDNAs. Reactions that included all three primers amplified an approximately 455 bp product in both leaves and pollen and an approximately 395 bp product in pollen only. Comparison with gene models indicated that the larger product was from RNAs transcribed from the annotated TSS and spliced using the upstream donor site, and Sanger sequencing confirmed this. However, comparison with gene models suggested that the lower-molecular weight product appearing only in pollen came from the downstream donor site (consistent with gene model AT3G06330.2) and an alternative TSS. If the annotated start site were used in conjunction with the downstream donor site, primers F1 and R would have amplified a larger product than the two we observed. These results suggested the models shown in Fig. S3B, and Sanger sequencing of the PCR products confirmed them. These results suggest two alternative TSSs were used in pollen but only the upstream site was used in leaves (Figs. S3B–S3C).

A similar result was obtained for RH11, an alternatively spliced gene encoding a DEAD/H-box RNA helicase. Inspection of the RNA-Seq alignments suggested there was at least one alternative TSS, and this alternative TSS enforced the use of an alternative donor site in the first intron (Fig. S4A). Variant-specific forward primers were used with a common reverse primer to amplify the differentially spliced regions. The variant-specific primers amplified product in both leaves and pollen, indicating that both samples used the alternative TSSs, but pollen and leaf samples each favored a different variant (Figs. S4B–S4C).

Lastly, we also found evidence for alternative TSSs in AT1G67480, which encodes a galactose oxidase/kelch repeat superfamily protein. The RNA-Seq data indicated that the annotated gene models were correct, but also suggested that other isoforms were expressed (Fig. S5A). The pollen data contained reads indicating a novel splice form, a retained-intron form, and two possible alternative TSSs. Semi-quantitative PCR amplification supported these alternative forms (data not shown). Primers designed to amplify the annotated, differentially spliced region produced amplicons with sizes that were consistent with predictions from the RNA-Seq data (Fig. S5B). Multiple variants were co-expressed in pollen and leaves, but their relative abundance was different between the sample types (Fig. S5C).

Semi-quantitative PCR analysis of AT2G19480, a NAP1 histone chaperone, failed to show pollen specific splicing

AT2G19480 (NAP1;4) encodes a NAP1 histone chaperone (Liu et al., 2009; Gao et al., 2012) with three splice variants that differ in the three prime region. According the RNA-Seq data, an exon-skipped form was preferred in pollen and the non-skipped, retained form was preferred in leaves (Fig. S6A). PCR primers were designed to amplify all forms (F and R2 in Fig. S6A) or just the exon-included form (F and R1 in Fig. S6A). However, the PCR amplifications contradicted the RNA-Seq results in that we found evidence for the exon-skipped form only in each sample (data not shown). Thus, the PCR assessment of alternative splicing confirmed existence and expression of the exon-skipped form in both pollen and leaves but was unable to identify the exon-retained form.

On-line data repository and RNA-Seq visualization resource

To enable readers to reproduce analyses, we assembled a git-based data and source code repository at http://www.bitbucket.org/lorainelab/pollenas. Readers familiar with the git version control program can clone the repository onto their local desktop, thus enabling them to view and retrieve older versions of the repository or use it as a basis for their own work. Tutorials on using git are available on the bitbucket.org Web site and many others. However, it is not necessary to use the git program to obtain a copy. Readers can also download the repository as a single “zip” by clicking the “download” icon on the project home page. Files can also be downloaded individually by clicking the “source” icon and then browsing the repository files and folders.

Analyses described above are organized into folders named for their contents. Typically, each folder contains a “results” sub-folder containing computer-generated results files and a “data” sub-folder containing primary data files obtained from the VisionWorks software or sequencing facility. Most analyses reported here were done using the R programming language and are documented in R “Markdown” (.Rmd) report-generation files that enable results and analysis code to be shown together in the same report file. Each “Markdown” file has a corresponding HTML report output file that was generated using the knitr R library, a tool for literate programming in R. As such, readers can either re-run the code themselves or can open the report file in a Web browser, after downloading the repository. Note that re-running analyses will also require installing some libraries that are not part of the base R distribution; this can be done using the “install packages” command within R.

We also provided RNA-Seq alignment files for visualization in Integrated Genome Browser as part of a pollen-specific IGBQuickLoad site at http://igbquickload.org/pollen. To view the data in IGB, readers first should obtain a copy of the software from the IGB web site at http://www.bioviz.org/igb. The software can be downloaded for free from the site; installers are available for Apple, Windows, and Linux computers. Once IGB is running, users can click the Arabidopsis image on the IGB home screen; this opens the Arabidopsis TAIR10 genome and loads the Arabidopsis reference gene annotations. To load RNA-Seq data shown in the figures, users can open the Pollen data source under the IGB Data Access tab. The data sets described here are named Pollen, Leaf 1, and Leaf 2. Alignments of sequences obtained from Sanger sequencing of amplicons can also be loaded. To view these data, readers can download the alignments file (named Amplicons.pslx) from the data and source code repository and open them in IGB as with any other file.

To illustrate the potential usefulness of this resource, we examined RNA-Seq alignments overlapping U2AF65B, one of two genes encoding U2AF65, the large subunit of the U2 snRNP accessory factor dimer, which helps to define 3′ splice sites by recruiting the U2 snRNP to the pre-mRNA. Genes that encode U2AF65 in Arabidopsis include U2AF65A (AT4G36690) and U2AF65B (AT1G60900). Previously, we showed that U2AF65A is differentially spliced in pollen as compared to leaf, but we did not report on alternative splicing of U2AF65B because there is only one model associated with this gene, and our splicing analysis pipeline detected differential splicing only in genes with multiple models. We therefore visually inspected RNA-Seq read alignments overlapping U2AF65B. Visualization of the RNA-Seq read alignments suggested that U2AF65B encodes at least one other splice variant thanks to alternative splicing of the fourth intron (Fig. 4). However, because the gene was not abundantly expressed in pollen, it was impossible to determine from RNA-Seq data alone whether the gene is differently spliced in pollen versus leaf.