Functional dynamic genetic effects on gene regulation are specific to particular cell types and environmental conditions

Genetic effects on gene expression and splicing can be modulated by cellular and environmental factors; yet interactions between genotypes, cell type and treatment have not been comprehensively studied together. We used an induced pluripotent stem cell system to study multiple cell types derived from the same individuals and exposed them to a large panel of treatments. Cellular responses involved different genes and pathways for gene expression and splicing, and were highly variable across contexts. For thousands of genes, we identified variable allelic expression across contexts and characterized different types of gene-environment interactions, many of which are associated with complex traits. Promoter functional and evolutionary features distinguished genes with elevated allelic imbalance mean and variance. On average half of the genes with dynamic regulatory interactions were missed by large eQTL mapping studies, indicating the importance of exploring multiple treatments to reveal previously unrecognized regulatory loci that may be important for disease

Functional dynamic genetic effects on gene regulation are specific to particular cell types and environmental conditions

Genetic effects on gene expression and splicing can be modulated by cellular and environmental factors; yet interactions between genotypes, cell type and treatment have not been comprehensively studied together. We used an induced pluripotent stem cell system to study multiple cell types derived from the same individuals and exposed them to a large panel of treatments. Cellular responses involved different genes and pathways for gene expression and splicing, and were highly variable across contexts. For thousands of genes, we identified variable allelic expression across contexts and characterized different types of gene-environment interactions, many of which are associated with complex traits. Promoter functional and evolutionary features distinguished genes with elevated allelic imbalance mean and variance. On average half of the genes with dynamic regulatory interactions were missed by large eQTL mapping studies, indicating the importance of exploring multiple treatments to reveal previously unrecognized regulatory loci that may be important for disease

CAG repeat expansion in the Huntington's disease gene shapes linear and circular RNAs biogenesis.

Alternative splicing (AS) appears to be altered in Huntington's disease (HD), but its significance for early, pre-symptomatic disease stages has not been inspected. Here, taking advantage of Htt CAG knock-in mouse in vitro and in vivo models, we demonstrate a correlation between Htt CAG repeat length and increased aberrant linear AS, specifically affecting neural progenitors and, in vivo, the striatum prior to overt behavioral phenotypes stages. Remarkably, a significant proportion (36%) of the aberrantly spliced isoforms are not-functional and meant to non-sense mediated decay (NMD). The expanded Htt CAG repeats further reflect on a previously neglected, global impairment of back-splicing, leading to decreased circular RNAs production in neural progenitors. Integrative transcriptomic analyses unveil a network of transcriptionally altered micro-RNAs and RNA-binding proteins (Celf, hnRNPs, Ptbp, Srsf, Upf1, Ythd2) which might influence the AS machinery, primarily in neural cells. We suggest that this unbalanced expression of linear and circular RNAs might alter neural fitness, contributing to HD pathogenesis

Differentially expressed genes presenting AS events in the striatum tissue of KI mouse models.

The bar graph presents the number of differentially expressed genes (DEG) which also are characterized by at least one significant AS event in the striatum of KI animal models of HD. Logarithmic fold change (LogFC) of their expression compared to Q20 control is reported in the y-axis. Transcripts are filtered based on significant p-value (p-value (TIFF)</p

Skipped exon (SE) events in the striatum (STR) are age and genotype specific.

The upset plot displays the number of SE events shared among different time points (2, 6 and 10 months) within the mouse striatum. For each time points, all the KI mice genotypes with expanded CAG tract were pooled together. The number of events within each intersection is presented in the vertical bars. Intersection groups (lines) or single time points (dots) are shown in the lower panel. Sample set size is indicated at the bottom left of the panel. (TIFF)</p

<i>Htt</i> CAG size dictates a specific linear AS signature in neural cells partly leading to NMD.

A) The Venn diagram reports the comparison between striatal (STR) and neuronal progenitors (NPC) SE events. All genotypes (Q20, Q80, Q92, Q111, Q140 and Q175) and time points (2, 6 and 10 months) for striatal districts and all genotypes (Q20, Q50, Q92, Q111) for NPC were combined. Shared SE events (18,85% of the total striatal SE events, p-value: 0, Fisher’s exact test) are indicated in the intersection. B) The histogram reveals shared GO&KEGG pathway between SE and NPC events, mainly involved in ‘cell junction’ and ‘synapse organization’. Terms are ordered by −log10(p−value). The number of genes in each term/pathway is indicated (numbers within columns). C) The Venn diagrams show the proportion of SE events in the striatum (STR) and neuronal progenitors (NPC) also annotated as NMD. More than one third of SE events (36.1% in STR and 36.2% in NPC) is annotated as NMD. The p-value of enrichment for each intersection is indicated. The tables for striatum and NPC, highlight the transcriptional dysregulation of Upf1, and Smg1, important regulators of the NMD pathway. D) A smaller proportion of SE events (7.8% in STR and 8.4% in NPC) intersected with transcripts subjected to m6A RNA modification (see methods). E-F) Motif analysis identified the splicing factors and/or RBP binding sites in the ± 100bp upstream and downstream adjacent regions to the alternatively spliced exons for the striatum (E) and for NPC (F). The p-value of enrichment testing for individual motifs in each data set is indicated. The candidate binding splicing factor and/or RBP family is shown.</p

Related to Fig 4.

The spreadsheets (n = 12) report 1. the circRNA and 2. Linear RNAs count in mouse neuronal progenitors (NPC) and 3. 4. embryonic stem cells (ESC) from the in vitro Htt CAG expanded alleles. Data are presented for all genotypes (Q20, Q50, Q92, Q111) for biological duplicate experiments. 5. Total list of dysregulated circRNAs in ESC and NPC. The comparison between the two most extreme Htt CAG expanded alleles (Q20 versus Q111) is reported. The list of 6. GO terms and 7. pathways associated to the genes originating differentially expressed circRNAs from the NPC Q20 versus Q111 comparison reported in Fig 4A is presented. 8. List of the 12 circRNAs satisfying the 3 stringent criteria of decreasing expression and negative correlation with Htt CAG, significantly different expression by circTest (see Methods for further details) is presented. CircRNA identification (ID), Chromosomal location (Chr), genomic coordinates, strand, overall circularization region and circRATIO [ratio circRNA / (circRNA+linear)] is reported for each molecule. 9. CISBP and 10. POSTAR motif analysis (see methods) of the sequences adjacent (100bp +/-) to circularizing points for the 12 stringently defined circRNAs (Fig 4D). 11. Mbnl1 CLIP biding sites on the 12 stringently defined circRNAs (Fig 4D). 12. The total list of transcriptionally dysregulated (DE) small RNA is presented for the ESC and NPC comparisons for the Q20 versus Q111. (XLSX)</p

Neuronal differentiation impact on small RNAs.

A) The coloured bar graphs report the number of small RNAs (>1 count/sample, see also Methods) comparing ESC and NPC conditions. Series of 4 Htt CAG expansion alleles (Q20, Q50, Q92 and Q111) are presented. Different classes of small RNAs are examined. Abbreviations as follows: microRNAs (miRNAs), Mitochondrial transfer RNAs (Mt-tRNAs), processed pseudogenes, ribosomal RNAs (rRNAs), small nucleolar RNAs (snoRNAs), small nuclear RNAs (snRNAs) small Cajal body-specific RNAs (scaRNAs), To be Experimentally Confirmed (TEC), transcribed processed pseudogenes. B) The heatmap presents the number of expressed (> 1 count per million, CPM) small RNAs among heterozygous Htt CAG knock-in HttQ20, HttQ50, HttQ92 and HttQ111 (Q20, Q50, Q92, Q111), comparing ESC and NPC lines. Different classes of small RNAs are examined. Abbreviations as in GENCODE transcript biotypes [86]. Color code bar (upper right) reports the number of expressed small RNAs in each condition. C) The bar chart shows the number of small RNAs differentially expressed between Htt Q111 versus Q20 genotypes. The comparison is presented for pluripotent (ESC) and neural committed progenitors (NPC). The number of small RNAs increasing (Increasing in Q111—upper part of the plot), and decreasing their expression in Q111 versus Q20 (Decreasing in Q111—lower part of the plot) is depicted. (TIFF)</p

Aberrant linear alternative splicing in the striatum of KI animal models of HD.

A) Bar graphs show the number of differential AS events in the striatum and cortex from 6 mouse KI models (Q20, Q80, Q92, Q111, Q140 and Q175) of HD, presenting different Htt CAG repeat lengths and 3 ages (2, 6, and 10 months). The number of events is shown for each genotype, time point and brain region. The inclusion level is calculated in comparison to Q20 controls and the positive or negative values are plotted. Source data by Langfelder P. et al (2016) [28]. Further details can be found in the Methods section. Each color of the bar chart represents a different AS type. B) Circos plot represents the number of transcripts within the striatum—showing differential AS events—shared between different genotypes (Q80, Q92, Q111, Q140 and Q175) and time points (2, 6, and 10 months). Conditions (genotypes and/or time points) sharing more than 50 transcripts are depicted in red. C) Weighted nodes graphical representation shows the functional enrichment analysis for transcripts displaying significant skipped exon (SE) events in the striatum. Highly expanded Htt CAG sizes (Q140 and Q175), the major contributors to aberrant SE in the striatum, are shown separated. Nodes’ size legend is depicted at the bottom. D) Representative agarose gel images (left) and quantification plots (right) report the RT-PCR results of AS validation for Crem, Adamts6, and Trpm4, selected transcript targets. RT-PCR assay and quantification were performed on striata RNA from an independent set of wild-type (WT), Q20 and Q111/Q175 mice. Transcripts isoforms with inclusion or exclusion of the variable exon (in yellow) are visualized and quantified. The plots report the PSI, percent-spliced-in. *p-value n = 3), error bars indicate standard deviation.</p

Validation of selected post-transcriptional regulators altered by <i>Htt</i> CAG expansion in mouse neuronal progenitors (NPC).

A) Dot plots and bar graphs report the RT-qPCR results of the validation for Upf1 and Mettl3 selected post-transcriptional splicing/back-splicing regulators (see also Fig 5). RT-qPCR assay and quantification were performed on RNA from Q20 and Q111 NPCs. Bar graphs plot the relative normalized mRNA expression levels (2-^ΔΔCt, see methods). Actβ Pgk1 was using as stable housekeeping gene. Error bars represent standard deviations from the mean of 3 biological replicates. *p-value n = 3). B) Representative western blot (WB) analysis reports the expression of Ptbp1, Ptbp2 and Ptbp3 proteins, normalized on Hsp90 housekeeping gene in Q20 and Q111 mNPCs. Molecular weights of proteins are indicated. C) The bar graph reports the quantification of the Ptbps proteins as detected by WB. While Ptbp1 and Ptbp3 are significantly less expressed in Q111 compared to Q20 NPCs, Ptbp2 is significantly more expressed in the HD condition. **p-value n = 4) D) The graph describes the relative normalized mRNA expression levels (2-^ΔΔCt, see methods). Pgk1 was used as stable housekeeping gene. Error bars represent standard deviations from the mean of biological replicates. *p-value n = 3). The results show that only Ptbp3 is significantly overexpressed in Q111 NPCs. E) The bar plot reports the RPM related to Ptbp1, 2 and 3 as detected in RNAseq data. These results confirmed the expression of Ptbps detected through RT-qPCR. (TIFF)</p