Dynamic regulation of alternative splicing and chromatin structure in Drosophila gonads revealed by RNA-seq

Abstract

Both transcription and post-transcriptional processes, such as alternative splicing, play crucial roles in controlling developmental programs in metazoans. Recently emerged RNA-seq method has brought our understanding of eukaryotic transcriptomes to a new level, because it can resolve both gene expression level and alternative splicing events simultaneously. To gain a better understanding of cellular differentiation in gonads, we analyzed mRNA profiles from Drosophila testes and ovaries using RNA-seq. We identified a set of genes that have sex-specific isoforms in wild-type (WT) gonads, including several transcription factors. We found that differentiation of sperms from undifferentiated germ cells induced a dramatic downregulation of RNA splicing factors. Our data confirmed that RNA splicing events are significantly more frequent in the undifferentiated cell-enriched bag of marbles (bam) mutant testis, but downregulated upon differentiation in WT testis. Consistent with this, we showed that genes required for meiosis and terminal differentiation in WT testis were mainly regulated at the transcriptional level, but not by alternative splicing. Unexpectedly, we observed an increase in expression of all families of chromatin remodeling factors and histone modifying enzymes in the undifferentiated cell-enriched bam testis. More interestingly, chromatin regulators and histone modifying enzymes with opposite enzymatic activities are coenriched in undifferentiated cells in testis, suggesting that these cells may possess dynamic chromatin architecture. Finally, our data revealed many new features of the Drosophila gonadal transcriptomes, and will lead to a more comprehensive understanding of how differential gene expression and splicing regulate gametogenesis in Drosophila. Our data provided a foundation for the systematic study of gene expression and alternative splicing in many interesting areas of germ cell biology in Drosophila, such as the molecular basis for sexual dimorphism and the regulation of the proliferation vs terminal differentiation programs in germline stem cell lineages. The GEO accession number for the raw and analyzed RNA-seq data is GSE16960.

Figure 1. Generation of gonadal transcriptomes using RNA-seq

(A) The four samples used for RNA-seq and four pair-wise comparisons in two categories to investigate: 1. sex-biased; 2. stage-specific patterns of gene expression and alternative splicing; (B) The anti-RNA Pol II ChIP-seq results using the bam testis sample. The five gene groups were classified according to their RPKM value based on RNA-seq results. The thresholds for low, moderate and high expression level were set up to reach the same number of genes in each group.

Figure 2. Enrichment of the chromatin regulators in undifferentiated cell-enriched bam testis

(A) Heat map showing genes that are uniquely expressed in both bam testis vs. wt testis, and bam ovary vs. wt ovary comparisons. Unique expression means that a certain gene is expressed in one sample (RPKM ≥1) but silent in another sample (RPKM <0.5). The RPKM cutoffs for expressed and silent genes are based on. Heat scale: log₂RPKM. To calculate the log₂ RPKM values of individual gene, all their original RPKM values were added by a pseudo-count of 1. (B) Enriched chromatin remodeling/modifying factors and histone modifying enzymes in bam testis. Abbreviations: PcG - Polycomb group complex, PRC1 and PRC2- Polycomb repressive complex 1 and 2, TrxG- Trithorax group complex, HMT- histone methyltransferase, HAT- histone acetyltransferase, HDAC- histone deacetylase, tTAFs- testis TAF homologs, tMAC- a testis-specific version of the MIP/DREAM complex, meiotic genes and TD- terminal differentiation genes for spermiogenesis. (C) In situ data using antisense riboprobes that recognize the trx (HMT) gene and lid (HDMT) gene in wt testis, respectively. There is a co-enrichment of the trx and lid transcripts in undifferentiated cells located at the tip of the wt testis, indicated by the black arrows. (D) In situ data using antisense riboprobes that recognize the Pcaf (HAT) gene and Rpd3 (HDAC) gene in wt testis, respectively. There is a co-enrichment of the Pcaf and Rpd3 transcripts in undifferentiated cells located at the tip of the wt testis, indicated by the black arrows.

Figure 3. Enhanced splicing activities in undifferentiated cell-enriched bam testis

(A) Expression of splicing factors is up-regulated in bam testis compared to wt testis. The red and green lines are the 2-fold cutoff lines. 177 splicing factors are 2-fold more enriched in bam testis relative to wt testis; while 26 splicing factors are 2-fold more enriched in wt testis relative to bam testis. The enrichment of splicing factors in bam testis is significant (P< 10⁻¹⁵). P-value was calculated using the chi-square test. (B) Enhanced overall splicing activities in bam testis compared to wt testis, shown by the CDF plot. P-value was calculated using one-sided Kolmogorov-Smirnov test. (C) Percentage of differentiation genes as single-isoform vs. multi-isoform genes [Differentiation genes are defined as silent in bam testis (RPKM in bam testis < 0.5); but are expressed in wt testis (RPKM in wt testis ≥ 1).]. Only 13% of the differentiation genes are multi-isoform genes, which is significantly (P< 10⁻¹⁵) lower than the 23% genome-wide proportion of multi-isoform genes, based on the Ensembl annotation. The P value was calculated using one-sided Fisher exact test.

Figure 4. Visualization of sex-specific isoforms of individual genes in bam and wt testis vs. ovary comparisons

(A) Genes with sex-specific isoforms in bam testis vs. bam ovary. All genes with ≥2 uniquely detectable isoforms (Methods) were analyzed for the presence (red) and absence (blue) of sex-specific isoforms. Presence (red) in “both” category indicated that the corresponding gene has at least one isoform expressed in both bam testis and bam ovary samples. The gene numbers of each class were labeled on the Y-axis of the panel. The RRR, RBR, BRR and RRB classes had at least one sex-specific isoform, therefore were labeled as black numbers (see Table S4-1 for details.). (B) Genes with sex-specific isoforms in wt testis vs. wt ovary. The presence (red) and absence (blue) of sex-specific isoforms for each individual gene were analyzed in wt testis and wt ovary samples; the presence (red) in both somatic samples indicated a lack of sex specificity for a particular isoform(s) of an individual gene (Table S4-2). (C) The wt ovary-biased p53 isoform is CG33336-RA or FBtr0084359. The wt testis-biased p53 isoform is CG33336-RB or FBtr0084360. (D) Real-time RT-PCR using isoform-specific primer sets showed that p53-RA is about 2-fold more enriched in wt ovary compared to wt testis; while p53-RB is approximately 14.6-fold more enriched in wt testis compared to wt ovary. The level of both p53-RA and p53-RB were normalized to the total p53 level, using a primer set that amplifies a common region of both isoforms (Methods). The total p53 level in wt ovary and wt testis was also shown, which has been normalized to the sample with a smaller RPKM (wt testis in this case, Table S1). (E) The wt ovary-biased Rab14 isoform is CG4212-RA or FBtr0080626. The wt testis-biased Rab14 isoform is CG4212-RB or FBtr0080627. (F) Real-time RT-PCR using isoform-specific primer sets showed that Rab14-RA is about 13.9-fold more enriched in wt ovary compared to wt testis; while Rab14-RB is approximately 171.4-fold more enriched in wt testis compared to wt ovary. The level of both Rab14-RA and Rab14-RB were normalized to the total Rab14 level, using a primer set that amplifies a common region of both isoforms (Methods). The total Rab14 level in wt ovary and wt testis was also shown, which has been normalized to the sample with a smaller RPKM (wt ovary in this case, Table S1).

Figure 5. Examples of characterized or putative splicing factors that have sex- or stage-specific isoforms

(A) UCSC snapshots show that the exu gene has sex-specific and sex-biased isoforms: exu-RC is wt testis-specific, exu-RA is wt testis-biased and exu-RB is wt ovary-biased isoforms, respectively. (B) Real-time RT-PCR results showed that exu-RA is 4.9-fold more enriched in wt testis compared to wt ovary; while exu-RB is approximately 94.3-fold more enriched in wt ovary compared to wt testis. The level of both exu-RA and exu-RB were normalized to the total exu level, using a primer set that amplifies a common region of all three isoforms (Methods). The total exu level in wt ovary and wt testis was also shown, which has been normalized to the sample with a smaller RPKM (wt ovary in this case, Table S1). (C) UCSC snapshots show that the imp gene is expressed in both undifferentiated spermatogonia-enriched bam testis and differentiating spermatids-containing wt testis. The imp gene has stage-specific isoforms: imp-RA/RB/RC are bam testis-specific, and imp-RG/RH are wt testis-specific. (D) Real-time RT-PCR results showed that imp-RA/RB/RC are 40.2-fold more enriched in bam testis compared to wt testis; while imp-RG/RH are approximately 1,966-fold more enriched in wt testis compared to bam testis. The level of both subsets of imp isoforms were normalized to the total imp level, using a primer set that amplifies a common region of all isoforms (Methods). The total imp level in bam testis and wt testis was also shown, which has been normalized to the sample with a smaller RPKM (bam testis in this case, Table S1).

Figure 6. Preferential gene distribution of testis terminal differentiation genes on the 2L chromosomal arm

Preferential chromosomal distribution of stage-biased genes in testis on the X chromosome and the 2L chromosomal arm: Percentage was calculated as the differentially expressed genes (≥ 2-fold change)/ total genes on a particular chromosome arm (X, 2L, 2R, 3L and 3R). To calculate the ratio of RPKM of bam testis/ wt testis or wt testis/ bam testis, the RPKM value was set to 0.5 if it less than 0.5. P-value was calculated with the Pearson’s chi-squared test (Method).

Figure 7. New features of the Drosophila gonadal transcriptome

(A) Distribution of identified CFTRs at the intronic and intergenic regions, based on the Ensembl annotation (http://www.ensembl.org/index.html) and the most update Flybase version (r5.19). *: Retained intronic CFTR must have a RPKM value greater than 10% of the RPKM(s) of the neighboring gene(s), to avoid contamination from pre-mRNAs. (B) Box plot of the expression level (log₂RPKM) of annotated expressed genes and CFTRs in all four samples.