Genetics of cell-type-specific post-transcriptional gene regulation during human neurogenesis

Abstract

The function of some genetic variants associated with brain-relevant traits has been explained through colocalization with expression quantitative trait loci (eQTL) conducted in bulk postmortem adult brain tissue. However, many brain-trait associated loci have unknown cellular or molecular function. These genetic variants may exert context-specific function on different molecular phenotypes including post-transcriptional changes. Here, we identified genetic regulation of RNA editing and alternative polyadenylation (APA) within a cell-type-specific population of human neural progenitors and neurons. More RNA editing and isoforms utilizing longer polyadenylation sequences were observed in neurons, likely due to higher expression of genes encoding the proteins mediating these post-transcriptional events. We also detected hundreds of cell-type-specific editing quantitative trait loci (edQTLs) and alternative polyadenylation QTLs (apaQTLs). We found colocalizations of a neuron edQTL in CCDC88A with educational attainment and a progenitor apaQTL in EP300 with schizophrenia, suggesting that genetically mediated post-transcriptional regulation during brain development leads to differences in brain function.

Keywords: RNA editing; alternative polyadenylation; genome-wide association studies; missing regulation; neurogenesis; quantitative trait loci.

Figure 1

Study design and cortical identity of hNPCs (A) Study design to identify cell-type-specific edQTLs and apaQTLs. (B) UMAP plot of 44,816 nuclei integrated across 4 phNPC lines from different donors differentiated for 4 and 8 weeks. Each dot represents a single nucleus, colored by its corresponding cluster. RG, radial glia (1–3 indicates different clusters based on Seurat algortihm); oRG, outer radial glia; tRG, truncated radial glia; PgS, S-phase progenitors; PgG2M, G2M phase progenitors; IPC, intermediate progenitor cell; EN, excitatory neuron; IT, inter-telencephalic; N, neuron; IN, inhibitory neuron; CGE, caudal ganglionic eminence. (C) Expression feature plots of select cell-type marker genes. HES1, radial glia; NRP1, excitatory newborn neurons; STMN2, pan-neuronal marker; GRIN2B, maturing neurons; CUX2, excitatory upper layer IT neurons; IL1RAPL2, excitatory L6 cortico-thalamic; GAD1, pan-inhibitory marker; CALB2, inhibitory CGE. (D) Cell type composition across all nuclei and for individual donor lines at 4 and 8 weeks of differentiation. (E) Heatmap showing ORA between phNPC and in vivo mid-gestation cell type markers. Enrichments reaching FDR-corrected significance (adjusted p value < 0.05) are colored and labeled with the odds ratio.

Figure 2

Features of RNA-editing sites (A) Proportions of RNA-editing events discovered in each cell type. (B) Pie chart showing proportions of RNA-editing events within Alu repeats, other repeats, or non-repeat regions in the genome per cell type. (C) Proportions of RNA-editing events across genomic regions. CDS, coding sequence; 3′ UTR, three prime untranslated region; and 5′ UTR, five prime untranslated region. (D) Overlap of RNA-edit sites with GTEx Cortex or BrainVar datasets. The edit sites overlapped are defined as annotated and the sites which do not overlap are defined as novel. (E) Local motif enrichment for annotated and novel edit sites per cell type. Number of edit sites (n) in each category is reported.

Figure 3

Cell-type-specific RNA editing during human neurogenesis (A) Comparison of Alu editing index (AEI) between progenitors and neurons. t test p value was reported. (B) Differential expression of ADAR1, ADAR2, and ADAR3 between progenitors and neurons. Adjusted p value (adj.pval) and log fold change (logFC) from limma were reported. (C) Overlap of RNA-editing sites discovered in progenitors, neurons, and fetal bulk data. (D) Enrichment of cell-type-specific RNA-editing sites within edit sites dysregulated in brain-relevant diseases. Number of edit sites overlap is reported, and the proportion of disease-specific edit sites overlapped with cell-type-specific edit sites is shown on the y axis. Asterisks indicate significant enrichment (FDR < 0.05). Autism (autism spectrum disorder), fragile X₁ (fragile X syndrome from UC Davis database), fragile X₂ (fragile X syndrome from NIH biobank dataset), GBM (glioblastoma), SCZ ACC (schizophrenia from anterior cingulate cortex), and SCZ DLPFC (schizophrenia from dorsolateral prefrontal cortex). (E) Editing rate of the edit site (chr1:3813767:A>G) at the 3′ UTR of CEP104 was positively correlated with the proportion of TUJ1+ neurons. Gene model for CEP104 and the genomic position of the edit site are given at the left. Scatterplot illustrating the correlation between edit rate and the relative abundance of TUJ1+ neurons, and correlation coefficient (r) and p values are shown. Representative immunocytochemistry images for TUJ1 (in green) and DAPI (in blue) staining of donors 18 and 321 (D18 and D321) are shown, scale bar is 100 μm.

Figure 4

Cell-type-specific edQTLs (A) Overlap of edSites detected in progenitors, neurons, and fetal bulk data. (B) Primary edQTLs were more significantly associated with editing as they were closer to the edit sites. Association p values at −log10 scale are shown on the y axis and distance from edit sites are shown on the x axis for progenitors (left) and neurons (right), density of the data points are indicated by color density per cell type. (C) Gene ontology results for edGenes found in neurons. (D) Enrichment of significant edQTLs within different RNA secondary structures illustrated at left. Enrichment p values at −log10 scale are shown on the y axis across structures and data are colored by cell types. The number of variants (n) significantly associated with edit sites and located within each structure are shown.

Figure 5

Cell-type-specific alternative polyadenylation and apaQTLs (A) Principal component analysis for alternative polyadenylation (APA) site usage colored by cell type. (B) Differentially expressed genes between progenitors and neurons that encode proteins playing a role in alternative polyadenylation. Fold change (logFC) is given on the x axis and data points are colored by adjusted p value at −log10 scale. logFC >0 indicates the genes upregulated in neurons and logFC < 0 indicates the genes upregulated in progenitors, and logFC = 0 is shown with dashed vertical line. (C) Two APA sites of CALM1 which were differentially expressed between progenitors and neurons are shown. Gene model for CALM1 is provided above, and relative read count per APA sites for each cell type are shown where relative read coverage calculated as ratio of number of count supporting each APA site to total number of reads including both APA sites are shown on the y axis; genomic position of the reads are shown on the x axis. Differential expression of APA₂ (longer 3′ UTR isoform) is shown between cell types. (D) Overlap of APA sites regulated by primary aSNPs across progenitors, neurons, and fetal bulk brain data. (E) Primary apaQTLs were more significantly associated with APA as they were closer to APA sites. Association p values at −log10 scale are shown on the y axis and distance from APA sites are shown on the x axis for progenitors (left) and neurons (right), density of the data points are indicated by color density per cell type. (F) apaQTL overlapping with canonical polyadenylation signal motif AAUAAA for RPL221 is shown for each cell type. Read coverage per genotype is shown.

Figure 6

Comparison of cell-type-specific molecular QTLs (A) Distribution of primary eQTLs, sQTLs, edQTLs, and apaQTLs from transcription start site (TSS) per cell type; the distance from TSS is shown on the x axis. (B) Distribution of primary eQTLs, sQTLs, edQTLs, and apaQTLs from transcription termination site (TTS) per cell type; the distance from TTS is shown on the x axis. (C) Distribution of primary eQTLs, sQTLs, edQTLs, and apaQTLs from splice sites per cell type, the distance from splice is shown on the x axis. Two distances were calculated relative to intron start and end sites, and the shortest distance was used for comparison for each QTL data. (D) Overlap of primary e/s/ed/apaQTLs via π₁ statistics (progenitors in purple and neurons in green). Matrices are colored based on the proportion of progenitor and neuron primary edSNP-edGene, aSNP-aGene, sSNP-sGene, and eSNP-eGene pairs that were non-null associations (π₁) in each of QTL datasets. We selected primary QTLs for each category (rows) to calculate π₁ using significance values from the molecular dataset listed in the columns. (E) pLOUEF values for aGenes, edGenes, eGenes, and sGenes per cell type are shown. p values from t test are reported.

Figure 7

Colocalization of cell-type-specific edQTLs and apaQTLs with brain-relevant trait GWAS loci (A) Genomics tracks illustrating that an edQTL within the CCDC88A locus in neurons was colocalized with an index variant for an education attainment GWAS, whereas there was not any significant eQTL for CCDC88A. Data points were colored based on the pairwise LD r² with rs563320407; p values for each association are shown on the y axis and genomic positions of genetic variants are shown on the x axis. (B) Boxplot illustrating distribution of editing rate across rs56320407. (C) The predicted secondary structure of the IRAlu hairpin for T and C alleles of rs56320407. Edit sites and genetic variants are indicated by red and blue colors, respectively. (D) Genomics tracks illustrating that an apaQTL within EP300 locus in progenitor was colocalized with an index variant for a schizophrenia GWAS, whereas there was not any significant eQTL for EP300. Data points were colored based on the pairwise LD r² with rs35508493; p values for each association are shown on the y axis and genomic positions of genetic variants are shown on the x axis. (E) Boxplot illustrating distribution of editing rate across rs56320407. (F) microRNA binding sites are shown for the genomic region which differs between two potential 3′ UTR isoforms of EP300 gene (genomic coordinates for APA site 1: chr22:41,178,956–41,180,079 and APA site 2: chr22:41,178,956–41,179,495).