Levin M et al. (2020), Gene expression dynamics are a proxy for select...

XB-ART-57343

Nucleic Acids Res 2020 Jun 19;4811:5926-5938. doi: 10.1093/nar/gkaa359.

Show Gene links Show Anatomy links

Gene expression dynamics are a proxy for selective pressures on alternatively polyadenylated isoforms.

Levin M , Zalts H , Mostov N , Hashimshony T , Yanai I .

???displayArticle.abstract???
Alternative polyadenylation (APA) produces isoforms with distinct 3'-ends, yet their functional differences remain largely unknown. Here, we introduce the APA-seq method to detect the expression levels of APA isoforms from 3'-end RNA-Seq data by exploiting both paired-end reads for gene isoform identification and quantification. We detected the expression levels of APA isoforms in individual Caenorhabditis elegans embryos at different stages throughout embryogenesis. Examining the correlation between the temporal profiles of isoforms led us to distinguish two classes of genes: those with highly correlated isoforms (HCI) and those with lowly correlated isoforms (LCI) across time. We hypothesized that variants with similar expression profiles may be the product of biological noise, while the LCI variants may be under tighter selection and consequently their distinct 3' UTR isoforms are more likely to have functional consequences. Supporting this notion, we found that LCI genes have significantly more miRNA binding sites, more correlated expression profiles with those of their targeting miRNAs and a relative lack of correspondence between their transcription and protein abundances. Collectively, our results suggest that a lack of coherence among the regulation of 3' UTR isoforms is a proxy for selective pressures acting upon APA usage and consequently for their functional relevance.

???displayArticle.pubmedLink??? 32421815
???displayArticle.pmcLink??? PMC7293032
???displayArticle.link??? Nucleic Acids Res

???attribute.lit??? ???displayArticles.show???

	Figure 1. APA-seq measures expression levels of distinct 3′ UTR isoforms in individual C. elegans embryos. (A) APA-seq identifies gene expression levels for distinct alternatively polyadenylated isoforms. APA-seq is an adaptation of the CEL-Seq method which utilizes paired-end reads: Read 1 contains a sample-specific barcode while Read 2 identifies the transcript. If sequenced long enough (100 bp in our case), Read 1 also provides information on the exact location of the polyadenylation site. APA-seq thus uses Read 2 to identify the expressed gene and then maps Read 1 to the gene specific region, thus enabling unique mapping in spite of the low sequencing quality that results from sequencing through the low-complexity poly-T region. (B) Plotting the 3′ends of Read 2 sequences results in a wide distribution within the gene (red peak) due to the fragmentation step of library preparation. However, the 3′ends of Read 1 sequences all map to the site immediately upstream of the poly-A tail thus producing a clear peak (black peak) when mapped to the gene sequence and revealing exact polyadenylation sites. The white boxes at the bottom of the distribution plots mark the coding sequence, while the green lines indicate the determined 3′ UTR regions. (C) Using the APA-seq method enables detection of the exact location of polyadenylation sites in C. elegans. The green lines at the bottom of each plot mark previously annotated 3′ UTRs (41) for the indicated gene, showing good agreement with the APA-seq Read 1 peaks (black). (D) Global comparison of the detected polyadenylation sites using APA-seq to the C. elegans UTRome annotation (41) shows high consistency between the two datasets. (E) The expression of two unique 3′ UTR isoforms for the C. elegans rps-29 gene throughout embryogenesis. Although in total (leftmost panel) the 3′ UTR variants show equal expression, examining expression across developmental stages shows predominant expression of the shorter 3′ UTR variant during early embryogenesis, and the inverse later in development.
	Figure 2. 3′ UTR isoform expression throughout C. elegans embryogenesis. (A) Heatmap showing the relative expression of 3′ UTR isoforms for 305 genes which passed overall expression and 3′ UTR isoform dynamics threshold throughout C. elegans embryogenesis. Each row in the heatmap corresponds to a gene and indicates the ratio between the expression of its short and long isoforms. Red and blue indicate the maximum and minimum 3′ UTR ratio for each gene, respectively. White and grey shadowed boxes indicate sets of developmental stages with similar isoform usage (based on Ward clustering, see clustergram on top of the heatmap) and identify the periods of major isoform switches. mpfc = minutes past four-cell stage. (B) The plots indicate the fold-change and P-value of difference in the isoform ratios for pairs of successive periods as defined by Ward clustering in (A). Red and blue colouring indicates significant elongation and shortening events, respectively. Grey colouring indicates insignificant changes. (C) Most of significant 3′ UTR isoform switches occur between the identified developmental switch periods. 89% (271) of the genes show significant changes in 3′ UTR usage in at least one of the identified period switches. (D) Clusters of significant 3′ UTR isoform changes indicate four main patterns—almost constant elongation or shortening over time (Clusters 1 and 4, respectively) or peaking shortening during mid-developmental transition and proliferative periods (Clusters 3 and 4, respectively). The thick black line indicates the mean 3′ UTR isoform ratio profile of all genes in the cluster and individual genes profiles are shown as thin grey lines.
	Figure 3. Genes with constitutive overall expression have dynamically expressed 3′ UTR isoforms. (A) Expression heatmaps indicating 3′ UTR isoform expression for 16 genes with multiple isoforms displaying varying levels of correlation between the expression of their 3′ UTR isoforms. The Pearson correlation coefficient r for each of the genes displayed is indicated at the top of each heatmap. (B) Relationship between 3′ UTR isoform expression correlations and the overall expression dynamics. Genes whose 3′ UTR isoform expression levels are correlated are more dynamic in their overall expression (r = 0.94, P = 0.03, second degree polynomial regression test). Dynamics for each gene is defined as the fold differences between its maximum and minimum expression values throughout the time-course. Red and gray horizontal bars represent the median and the interquartile ranges of the data, respectively. The last bin with highest 3′ UTR isoform correlations exhibit significantly higher overall expression dynamics than the preceding bin (P< 10−20, Mann–Whitney test). (C) Relationship between 3′ UTR isoform expression correlations and the respective overall total mRNA expression levels. Genes whose 3′ UTR isoform expression levels are uncorrelated show significantly higher overall expression levels (r = 0.99, P = 0.002, second degree polynomial regression test). Red and grey horizontal bars represent the median and the interquartile ranges of the data, respectively. (D) Heatmaps in the two leftmost panels show standardized expression of the 3′ UTR isoforms of 545 genes belonging to the group of genes whose 3′ UTR isoform's expression show high correlation (HCI genes). The right panel shows a heatmap of the total mRNA expression (on a log10 scale) of the same genes confirming that genes with highly correlated isoform expression are also dynamically expressed. The time points are the same as indicated in Figure 2A. (E) Same as D for 79 genes belonging to the group of genes whose 3′ UTR isoform's expression show low correlation (LCI genes). Genes with distinct expression profiles of 3′ UTR isoforms appear constitutively expressed when examined at the total gene expression level.
	Figure 4. Genes with multiple, lowly correlated 3′ UTR variants show evidence for increased miRNA regulation. (A) The left panel shows a schematic representation of the miRNA analysis. The length of the 3′ UTR region, the number of basic miRNA seed matches, and the 3′ UTR region unique to the longer 3′ UTR isoform were considered. The right panel shows two examples of genes and their miRNA binding sites. (B) Boxplots indicating the number of miRNA binding sites in the full 3′ UTR regions of genes with single or multiple 3′ UTR isoforms, highly correlating (HCI) or lowly correlating isoforms (LCI). Genes with multiple 3′ UTR variants have significantly more miRNA targets than genes with a single variant (P< 10−12, Mann–Whitney test). Between multiple 3′ UTR isoform genes, LCI genes have significantly more miRNA binding sites than HCI genes (P< 0.02, Mann–Whitney test). (C) Same as B for the full 3′ UTR lengths of genes. Genes with multiple 3′ UTR variants have significantly longer 3′ UTRs than single isoform genes (P< 10−28, Mann–Whitney test). Within multiple isoform genes, LCI genes have significantly longer 3′ UTRs than HCI genes (P< 10−3, Mann–Whitney test). (D) Boxplots indicating the number of miRNA binding sites in the unique 3′ UTR region of the longer 3′ UTR isoform across HCI and LCI genes. LCI genes have significantly more miRNA binding sites in their unique 3′ UTR region than HCI genes (P< 10−3, Mann–Whitney test). (E) Boxplots indicating the length of the unique 3′ UTR region across HCI and LCI genes. LCI genes have significantly longer unique 3′ UTR regions than HCI genes (P< 10−4, Mann–Whitney test). (F) Correlating expression of miRNA expression with expression ratio dynamics between 3′ UTR isoforms. The black line shows the expression profile of cel-miR-1–5p miRNA throughout the developmental time-course. Depicted in red is the ratio between the expression profiles of the two 3′ UTR isoforms (short/long) of the F07C6.4 gene. The Pearson correlation coefficients between 3′ UTR isoform ratio and miRNA expression of all target genes were used for the analysis shown in G and H. (G) For the indicated miRNAs, the boxplots show the distribution of correlations between miRNA expression and 3′UTR isform ratio showing a significant difference between the HCI and LCI gene groups (out of 64 miRNAs with dynamic expression, see Supplementary Table S2 for statistics for all miRNAs). All ten exhibit a positive median correlation between miRNA expression and the 3′ UTR isoform expression ratio of genes it is predicted to bind (as in F). Further, this correlation is significantly higher in LCI than in HCI genes. (H) Heatmap of the standardized expression dynamics of the ten indicated miRNAs which show differences between HCI and LCI genes across development.
	Figure 5. Genes with lowly correlated 3′ UTR isoforms (LCI genes) have a lower correspondence between total mRNA and protein expression. (A) Violin plots indicate the distribution of Pearson correlation coefficients between total mRNA and protein expression levels across developmental stages for C. elegans genes with one and multiple 3′ UTR isoforms, highly correlating (HCI), and lowly correlating isoforms (LCI). (B, C) Same as A, for LCI and HCI genes in Drosophila melanogaster (B) and Xenopus laevis (C) embryonic development.
	Figure 2. 3′ UTR isoform expression throughout C. elegans embryogenesis. (A) Heatmap showing the relative expression of 3′ UTR isoforms for 305 genes which passed overall expression and 3′ UTR isoform dynamics threshold throughout C. elegans embryogenesis. Each row in the heatmap corresponds to a gene and indicates the ratio between the expression of its short and long isoforms. Red and blue indicate the maximum and minimum 3′ UTR ratio for each gene, respectively. White and grey shadowed boxes indicate sets of developmental stages with similar isoform usage (based on Ward clustering, see clustergram on top of the heatmap) and identify the periods of major isoform switches. mpfc = minutes past four-cell stage. (B) The plots indicate the fold-change and P-value of difference in the isoform ratios for pairs of successive periods as defined by Ward clustering in (A). Red and blue colouring indicates significant elongation and shortening events, respectively. Grey colouring indicates insignificant changes. (C) Most of significant 3′ UTR isoform switches occur between the identified developmental switch periods. 89% (271) of the genes show significant changes in 3′ UTR usage in at least one of the identified period switches. (D) Clusters of significant 3′ UTR isoform changes indicate four main patterns—almost constant elongation or shortening over time (Clusters 1 and 4, respectively) or peaking shortening during mid-developmental transition and proliferative periods (Clusters 3 and 4, respectively). The thick black line indicates the mean 3′ UTR isoform ratio profile of all genes in the cluster and individual genes profiles are shown as thin grey lines.
	Figure 3. Genes with constitutive overall expression have dynamically expressed 3′ UTR isoforms. (A) Expression heatmaps indicating 3′ UTR isoform expression for 16 genes with multiple isoforms displaying varying levels of correlation between the expression of their 3′ UTR isoforms. The Pearson correlation coefficient r for each of the genes displayed is indicated at the top of each heatmap. (B) Relationship between 3′ UTR isoform expression correlations and the overall expression dynamics. Genes whose 3′ UTR isoform expression levels are correlated are more dynamic in their overall expression (r = 0.94, P = 0.03, second degree polynomial regression test). Dynamics for each gene is defined as the fold differences between its maximum and minimum expression values throughout the time-course. Red and gray horizontal bars represent the median and the interquartile ranges of the data, respectively. The last bin with highest 3′ UTR isoform correlations exhibit significantly higher overall expression dynamics than the preceding bin (P< 10−20, Mann–Whitney test). (C) Relationship between 3′ UTR isoform expression correlations and the respective overall total mRNA expression levels. Genes whose 3′ UTR isoform expression levels are uncorrelated show significantly higher overall expression levels (r = 0.99, P = 0.002, second degree polynomial regression test). Red and grey horizontal bars represent the median and the interquartile ranges of the data, respectively. (D) Heatmaps in the two leftmost panels show standardized expression of the 3′ UTR isoforms of 545 genes belonging to the group of genes whose 3′ UTR isoform's expression show high correlation (HCI genes). The right panel shows a heatmap of the total mRNA expression (on a log10 scale) of the same genes confirming that genes with highly correlated isoform expression are also dynamically expressed. The time points are the same as indicated in Figure 2A. (E) Same as D for 79 genes belonging to the group of genes whose 3′ UTR isoform's expression show low correlation (LCI genes). Genes with distinct expression profiles of 3′ UTR isoforms appear constitutively expressed when examined at the total gene expression level.
	Figure 4. Genes with multiple, lowly correlated 3′ UTR variants show evidence for increased miRNA regulation. (A) The left panel shows a schematic representation of the miRNA analysis. The length of the 3′ UTR region, the number of basic miRNA seed matches, and the 3′ UTR region unique to the longer 3′ UTR isoform were considered. The right panel shows two examples of genes and their miRNA binding sites. (B) Boxplots indicating the number of miRNA binding sites in the full 3′ UTR regions of genes with single or multiple 3′ UTR isoforms, highly correlating (HCI) or lowly correlating isoforms (LCI). Genes with multiple 3′ UTR variants have significantly more miRNA targets than genes with a single variant (P< 10−12, Mann–Whitney test). Between multiple 3′ UTR isoform genes, LCI genes have significantly more miRNA binding sites than HCI genes (P< 0.02, Mann–Whitney test). (C) Same as B for the full 3′ UTR lengths of genes. Genes with multiple 3′ UTR variants have significantly longer 3′ UTRs than single isoform genes (P< 10−28, Mann–Whitney test). Within multiple isoform genes, LCI genes have significantly longer 3′ UTRs than HCI genes (P< 10−3, Mann–Whitney test). (D) Boxplots indicating the number of miRNA binding sites in the unique 3′ UTR region of the longer 3′ UTR isoform across HCI and LCI genes. LCI genes have significantly more miRNA binding sites in their unique 3′ UTR region than HCI genes (P< 10−3, Mann–Whitney test). (E) Boxplots indicating the length of the unique 3′ UTR region across HCI and LCI genes. LCI genes have significantly longer unique 3′ UTR regions than HCI genes (P< 10−4, Mann–Whitney test). (F) Correlating expression of miRNA expression with expression ratio dynamics between 3′ UTR isoforms. The black line shows the expression profile of cel-miR-1–5p miRNA throughout the developmental time-course. Depicted in red is the ratio between the expression profiles of the two 3′ UTR isoforms (short/long) of the F07C6.4 gene. The Pearson correlation coefficients between 3′ UTR isoform ratio and miRNA expression of all target genes were used for the analysis shown in G and H. (G) For the indicated miRNAs, the boxplots show the distribution of correlations between miRNA expression and 3′UTR isform ratio showing a significant difference between the HCI and LCI gene groups (out of 64 miRNAs with dynamic expression, see Supplementary Table S2 for statistics for all miRNAs). All ten exhibit a positive median correlation between miRNA expression and the 3′ UTR isoform expression ratio of genes it is predicted to bind (as in F). Further, this correlation is significantly higher in LCI than in HCI genes. (H) Heatmap of the standardized expression dynamics of the ten indicated miRNAs which show differences between HCI and LCI genes across development.
	Figure 5. Genes with lowly correlated 3′ UTR isoforms (LCI genes) have a lower correspondence between total mRNA and protein expression. (A) Violin plots indicate the distribution of Pearson correlation coefficients between total mRNA and protein expression levels across developmental stages for C. elegans genes with one and multiple 3′ UTR isoforms, highly correlating (HCI), and lowly correlating isoforms (LCI). (B, C) Same as A, for LCI and HCI genes in Drosophila melanogaster (B) and Xenopus laevis (C) embryonic development.

References [+] :

Alt, Synthesis of secreted and membrane-bound immunoglobulin mu heavy chains is directed by mRNAs that differ at their 3' ends. 1980, Pubmed