Genome-wide analysis of alternative splicing in Zea mays: landscape and genetic regulation

Abstract

Alternative splicing enhances transcriptome diversity in all eukaryotes and plays a role in plant tissue identity and stress adaptation. To catalog new maize (Zea mays) transcripts and identify genomic loci that regulate alternative splicing, we analyzed over 90 RNA-seq libraries from maize inbred lines B73 and Mo17, as well as Syn10 doubled haploid lines (progenies from B73 × Mo17). Transcript discovery was augmented with publicly available data from 14 maize tissues, expanding the maize transcriptome by more than 30,000 and increasing the percentage of intron-containing genes that undergo alternative splicing to 40%. These newly identified transcripts greatly increase the diversity of the maize proteome, sometimes coding for entirely different proteins compared with their most similar annotated isoform. In addition to increasing proteome diversity, many genes encoding novel transcripts gained an additional layer of regulation by microRNAs, often in a tissue-specific manner. We also demonstrate that the majority of genotype-specific alternative splicing can be genetically mapped, with cis-acting quantitative trait loci (QTLs) predominating. A large number of trans-acting QTLs were also apparent, with nearly half located in regions not shown to contain genes associated with splicing. Taken together, these results highlight the currently underappreciated role that alternative splicing plays in tissue identity and genotypic variation in maize.

Figure 1.

Overlap of Computationally Predicted Transcripts Generated from IBM Syn10 Mapping Population Libraries Compared with Public B73 Tissue Libraries. (A) Overlap of new isoforms arising from known genes predicted from separate analyses of IBM Syn10 and public libraries. (B) Overlap of entirely new genes predicted from separate analyses of IBM Syn10 and public libraries.

Figure 2.

Examples of Proteins Encoded by Known Transcripts (Top) Compared with Those Encoded by Novel Transcripts (Bottom). (A) 3BETAHSD/D2’s loss of reticulon ER localization domain. (B) IRE1’s gain of unfolded protein-responsive luminal dimerization domain. (C) WAK3’s gain of EF hand calcium binding activation domain. (D) PRP8’s loss of RNA and U5/U6 interaction domains. (E) GRMZM2G119248’s switch from asparagine synthase to a putative bromodomain-containing transcription factor. [See online article for color version of this figure.]

Figure 3.

Relative Expression Levels of Isoforms for Genes with Substantial Protein or miRNA Binding Modifications in the Novel Transcript Set. Relative gene expression differences between tissues are denoted in blue. Relative isoform abundance within each tissue is denoted in red. (A) One new isoform of GRMZM2G119248 encoding a putative transcription factor. (B) One new isoform of GRMZM2G086430 is targeted by miR827. (C) One new isoform of GRMZM2G166976 loses the target binding site for miR827. (D) One new isoform of GRMZM2G034876 loses the target binding site for miR396a. (E) Novel pollen-specific gene targeted by miR164c/h.

Figure 4.

Gain and Loss of miRNA Target Sites in the Novel Transcript Set. Novel transcripts (gray) that gained or lost miRNA target sites relative to their most similar known isoforms (black) are shown. (A) Two new isoforms of GRMZM2G357595 are targeted by miR159c. (B) One new isoform of GRMZM2G086430 is targeted by miR827. (C) One new isoform of GRMZM2G166976 loses the target binding site for miR827. (D) One new isoform of GRMZM2G034876 loses the target binding site for miR396a. (E) Novel pollen-specific gene targeted by miR164c/h.

Figure 5.

Single Tissue-Specific Transcripts. Number of known and novel transcripts that are expressed at or above 1.3 FPKM in only one tissue and <0.1 FPKM in all others. Gray areas represent the portion of transcripts that are only expressed in one tissue whose genes are also only expressed in that tissue (gene expression dependent). Black areas represent the portion of transcripts that are only expressed in one tissue, with other isoforms present in other tissues (alternative splicing dependent). Transcript abundances were normalized by average transcript expression per library to account for variable tissue transcriptome complexity.

Figure 6.

Relative Expression Level of Isoforms for Genes with Differential Genotype-Dependent Alternative Splicing Variations in B73, Mo17, and IBM Syn10 Lines. Relative gene expression between genotypes is denoted in blue. Relative isoform abundance within each genotype is denoted in red. (A) GRMZM2G154278, which has a strong cis-acting QTL with clean segregation. (B) GRMZM2G040587, which has a strong cis-acting QTL with clean segregation. (C) GRMZM2G155232, which has a strong cis-acting QTL. (D) GRMZM2G136455, which has a strong trans-acting QTL.

Figure 7.

Mapping Results Reveal QTLs Regulating Genotype-Dependent Alternative Splicing Variations. Genes with statistically significant (q < 0.05) splicing differences between B73 and Mo17 were selected for mapping using Spotfire. The relative percentage of total gene expression that each isoform represents was used as the phenotype and compared with marker data from the Illumina MaizeSNP50 DNA analysis kit (Illumina). (A) GRMZM2G154278, exemplifying the case of a typical cis-acting QTL, with one clear peak surrounding the gene’s annotated position on chromosome 8. (B) GRMZM2G136455, exemplifying a trans-acting QTL. The QTL is located on chromosome 5, while the gene is located on chromosome 1.

Figure 8.

Splicing Factor EMB2769 Is Differentially Spliced in B73 Compared with Mo17, Which Has a 36-bp Deletion That Removes the AG 3′ Consensus Acceptor Sequence. Translation start site for the longest possible protein is denoted with blue arrows while the stop codon for that protein is denoted in red.