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