Genome-wide Analysis of Gene Expression in Cortex Using High-Throughput Total RNA Sequencing.

To examine transcriptional regulation of AD specific gene expression, we performed high-throughput total RNA sequencing (RNA-seq) analysis of cortices from 6-mo-old wild-type (WT) and 5xFAD mice. Since the expression levels of many genes are programed and regulated in tissues and blood of AD mouse models (2021–22), genes with increased or decreased expression levels in AD might be functionally important in the brain. We isolated the cortex region from mouse whole brain and extracted RNA for total RNA-seq (Fig. 1A). Data from the cortices of WT and 5xFAD mice were obtained to identify differentially expressed genes.

Fig. 1.

High-throughput total RNA-seq profile for PLCγ1 and PLCβ subfamily gene expression in the cortex. (A) High-throughput total RNA-seq analysis workflow. (B) Heat map and hierarchical clustering of total RNA-seq analysis from mouse cortex. The 17,109 genes are DEGs between WT and 5xFAD. (C) Normalization of total RNA-seq using scatterplot matrices from WT and 5xFAD model mouse cortex. The x axis indicates gene expression of WT; y axis indicates gene expression of AD. (D) Genome browser view of the PLCγ1, PLCβ1, PLCβ3, and PLCβ4 genes genomic locus and expression display with total RNA-seq data in cortex. ER, End Repair; ssCirc, single-strand Circle; DNB, DNA NanoBall; PE 100, Paired end 100; Cons, Consensus.

We confirmed that the expression levels of 17,109 genes were 5xFAD specific (Fig. 1B). Overall, 1,472 genes were significantly up-regulated, whereas 653 genes were down-regulated in the 5xFAD mouse cortex (Fig. 1C and SI Appendix, Fig. S1). Previous studies have reported that dysregulations of PLC genes are closely associated with many brain disorders, such as schizophrenia, bipolar disorder, Huntington’s disease, depression, and AD (10). Although the importance of the correlation between PLC genes and AD has been suggested, supporting evidence such as genes exhibiting differential expression or variation between healthy individuals and those with AD have not been identified. To demonstrate the expression levels of PLCβ1, PLCβ2, PLCβ3, PLCβ4, and PLCγ1 genes in the brains of WT and 5xFAD mice, we performed GWAS for PLCβ1, PLCβ2, PLCβ3, PLCβ4, and PLCγ1 genes using high-throughput total RNA-seq data to measure their transcription levels. The expression levels of PLCβ1, PLCβ2, PLCβ3, PLCβ4, and PLCγ1 were slightly up-regulated in the 5xFAD cortex, while PLCβ2 expression was up-regulated by ∼2.4 times in 5xFAD (Fig. 1D and SI Appendix, Fig. S2B). However, endogenous expression level of PLCβ2 was significantly lower than those of PLCβ1, PLCβ3, PLCβ4, and PLCγ1 in the 5xFAD cortex (SI Appendix, Fig. S2). In addition to the analysis of PLCβ1, PLCβ2, PLCβ3, PLCβ4, and PLCγ1 transcription levels, Western blot analysis was performed to confirm changes at the translational level, revealing that there were no significant differences among the levels of the proteins encoded by these genes in WT and 5xFAD cortices (SI Appendix, Fig. S3). Notably, PLCβ2 protein was not detected in the Western blot analysis because its endogenous level in neurons was very low (23).

Single-nucleotide polymorphisms (SNPs) in APOE or TREM2 have been proposed to play pathogenic roles in AD. Moreover, somatic mutations in brain have been suggested to be associated with genetic architecture in AD (2425–26). RNA-seq data are benefitted by the detection of unidentified alternative variations and specific gene expression (27). Though the expression levels of PLC subfamily proteins showed no significant differences between healthy individuals and those with AD, the evidence led us to analyze the genomic variations of PLC genes in AD. To obtain the profiles of PLC subfamily genes based on SNV data from total RNA-seq, we proposed a pipeline for RNA-seq−based SNV and insertion/deletion (InDel) analyses (Fig. 2). We performed gene expression and SNV profiling of cortex from brains of WT and AD mice using total RNA-seq to obtain unknown information on SNVs in PLCγ1, PLCβ1, PLCβ2, PLCβ3, and PLCβ4, and provide insight into genetic variation in the pathogenesis of AD.

Fig. 2.

Total RNA-seq−based SNV identification workflow for AD. Single and paired-end reads are calculated from NGS (Next-generation sequencing). Mapping was performed with reference genes and each of bam files was indexed with SAMtools. Preprocessing formed Ras sequencing data. SAMtools mpileup was performed with UCSC (University of California Santa Cruz) mm10 genome as reference, and then analyzing the SNV/InDel for AD. VCF files were generated with vcfutils.pl varFilter, and functional annotation of each variant was performed with ANNOVAR (ANNOtate VARiation) software.

Identification and Characterization of AD-Specific Mouse PLCγ1 and PLCβ SNVs in the Cortex.

High-throughput total RNA-seq is one of the most potent techniques for determining SNVs and genetic variants of target genes. A total of 163 variations (132 SNVs and 31 InDels) were identified in five genes. Of these, most were located in introns (152 variants), and only five variations were located in exons; the remaining six variants were located in splicing regions. In the 5xFAD cortex, PLCγ1 had 13 variants, while PLCβ1, PLCβ2, PLCβ3, and PLCβ4 had 66, 6, 3, and 75 variants, respectively (SI Appendix, Tables S1–S5). In this study, we mainly focused on variations in the exons of PLCγ1. There were five functional variations in the exons, including three synonymous SNV, one frameshift insertion in exon 27 of PLCγ1, and one stop gain in exon 11 of PLCβ3 (SI Appendix, Tables S1 and S4). A frameshift insertion, wherein a single “A” nucleotide was inserted at position 160,759,682 in chromosome 2, gave rise to an abnormal codon that substituted the isoleucine at the 970th amino acid position to asparagine (Fig. 3A) in PLCγ1 protein of the AD mouse model. Stop-gain mutation occurred due to the substitution of the reference “G” nucleotide with “A” at position 6,963,413 in chromosome 2, encoding a terminal codon due to a change of glutamine at the 352nd amino acid of PLCβ3 in the 5xFAD cortex (SI Appendix, Table S4).

Fig. 3.

Identification and characterization of AD-specific sequenced mouse PLCγ1. (A) Mapping for AD model mouse-dependent frameshift at exon 27 of PLCγ1 gene body. One single-nucleotide insertion was caused by abnormal codon, by breaking 970th amino acid, isoleucine, into asparagine. (B) UCSC genome browser sequence alignment of various species, including human for PLCγ1 gene. Both exon 27 of PLCγ1 for mouse and exon 26 of human PLCγ1 was well conserved in isoleucine site, and full length of amino acid sequence was completely matched. Red vertical bars indicate the locations of the isoleucine for both human and mouse PLCγ1 gene.

However, we found no differences in the translation and termination of PLCγ1 and PLCβ3 between WT and the 5xFAD cortex (SI Appendix, Fig. S3). Eukaryotic AS events are closely related to adult tissue functions, organ development, and tissue homeostasis (28, 29); alterations in any of these can affect cell proliferation, methylation, and migration of cancer cells. Aberrant splicing is implicated as an important factor contributing to diseases (30). To determine whether exon 27 of PLCγ1 affects AS after mRNA processing, we performed RT-PCR using various primer sets targeting the region near this exon. Exons 26 to 30 of PLCγ1 were deleted in the mature mRNA of 5xFAD cortex. Conversely, exons 27 to 32 of PLCγ1 were conserved in the pre-mRNA of both WT and 5xFAD cortex (SI Appendix, Fig. S4), suggesting that exon skipping or AS occurs at the PLCγ1 SNV region in AD. Single amino acid substitutions have been linked to many human diseases; these disease-causing mutations tend to cause changes in hydrogen-bonded linkages or bridges, resulting in harmful amino acid mutations (31).

The frequency of AS differs among species. For example, although both human and mouse genomes have similar numbers of conserved motifs, the average alternative pre-mRNA splicing rates of human (>95 to 100%) and mouse (∼63%) genes are different (32). This comparison between human and mouse genes has also revealed that alternative exons of the reading frame are conserved in human and mouse, indicating that genetic changes are caused by AS (32). We further examined whether the 970th amino acid, isoleucine, at exon 27 of mouse PLCγ1 was conserved in other species. We obtained the PLCγ1 sequences of various species using the UCSC genome browser. Sequence alignment of PLCγ1 from various species, including human, revealed that exon 27 of mouse PLCγ1 was completely conserved across species. Moreover, exon 26 of human PLCγ1 was well conserved with respect to the isoleucine site. Surprisingly, the full-length amino acid sequences encoded by exon 26 of human and exon 27 of mouse PLCγ1 matched completely (Fig. 3B). Collectively, our findings suggest that the SNV in exon 27 of mouse PLCγ1 plays an important role in AS during mRNA processing, and may have potential as a biomarker for prediction of AD.

High-Resolution Profiling of Histone Modification at PLCγ1 in the Forebrain during Brain Development.

Histone acetylation facilitates transcriptional activity and plays an important role in the development of human brain diseases through epigenetic modification events, such as altering chromatin structure or accessibility of transcription factors (33, 34). To determine whether epigenetic changes may occur in the PLCγ1 gene during mouse brain development, we performed H3K27ac enrichment profiling using chromatin immunoprecipitation sequencing (ChIP-seq) dataset to generate high-resolution histone modifications. Histone acetylation gradually accumulated in PLCγ1 gene in mouse forebrain regions ranging from E11.5 to P56. Histone acetylation increased in the introns of PLCγ1, and AD-associated SNVs were also concentrated in this region (Fig. 4). In addition, AD-specific nucleotide alternation sites were distributed in the noncoding region of PLCγ1, and we found several transcription factor binding motifs at AD-specific AS sites (SI Appendix, Fig. S5). These results suggest that SNVs in PLCγ1 are concentrated in the introns in AD, indicating that they may interfere with the epigenetic role associated with AD.

Fig. 4.

High-resolution binding profiles of histone acetylation in the forebrain during the brain development. H3K27ac binding profiles throughout the mouse forebrain development. Histone acetylation was highly accumulated at PLCγ1 gene body in P56 mouse forebrain. AD-dependent expressed SNVs were concentrated in the intron region of the PLCγ1 gene, which is enrichment of histone acetylation. Blue vertical bars indicate the locations of the AD-dependent SNVs region of the PLCγ1 gene. Black bold letters indicate the nucleotides for WT, red bold letters indicate the insertion or substitution nucleotide for AD, and gray letters indicate the deleted nucleotide for AD.

Prediction of AD-Specific Nucleotide Alteration Sites in the Human Genome Using Deep Learning.

Dysregulation of AS has been implicated in AD (35, 36). Recently, a new deep learning-based tool called SpliceAI, which can predict splice junctions of target genes from pre-mRNA nucleotide sequences with high accuracy, was introduced (4). Total RNA-seq−based SNV analysis confirmed the profile of AD-specific SNVs in PLCγ1 in 5xFAD mouse cortex. Then, we further examined whether a single-nucleotide alteration in human PLCγ1 may be correlated with the 5xFAD cortex. After SpliceAI analysis predicted a total of 14 splicing sites in the human PLCγ1 gene, accurate delta scores and positions were analyzed (Fig. 5), and a novel splicing site in exon 26 of human PLCγ1 was identified. The single nucleotide “G” at position 41,172,421 in chromosome 20 was substituted with “A,” “C,” or “T” at position 41,172,423, which was completely correlated with exon 27 of mouse PLCγ1 in the AD model cortex. SpliceAI analysis accurately revealed that the SNVs in exon 26 of human PLCγ1 and exon 27 of mouse PLCγ1 were correlated. These AD-associated SNVs occurred at the same position in human and the 5xFAD cortex. Taken together, these results clearly demonstrated the prediction accuracy of SpliceAI in humans with respect to genetic mutations and splicing in the AD model.

Fig. 5.

Prediction of AD specific nucleotide alteration sites from human genomic sequence with deep neural network. (Top) The schematic diagrams of the deep neural network procedure from human genomic sequence database through the SpliceAI analysis. (Bottom) The total 14 accurate prediction splicing sites in the human PLCγ1 gene, with the delta scores and position information.

Animals.

The 5xFAD (6-mo-old) transgenic mice were used as model of AD. The 5xFAD mice were obtained from the Jackson Laboratory. All animal experiments performed in this study were reviewed and approved by the Institutional Animal Care Use (IACUC) committee at KBRI (Korea Brain Research Institute) (IACUC-20-00018).

Total RNA-seq.

Total RNA-seq was performed as previously published (20). Total RNA extraction from mouse cortex was performed using commercial methods based on TRIzol (Invitrogen). Libraries construction was performed using TruSeq Stranded Total RNA LT Sample Prep Kit (Human Mouse Rat) by manufacturer’s instructions (Illumina). Firstly, to make a short fragment of mRNA, we added the fragmentation buffer. The oligo dT-primer was used to synthesize the first-strand complementary DNA (cDNA) to take short fragments as templates. Preparation of synthesis for second-strand cDNA was performed using buffer, 2′-deoxynucleoside 5′-triphosphate (dNTPs, containing 2′-deoxyuridine 5′-triphosphate instead of thymidine 5′-triphosphate), RNaseH, and DNA polymerase I, respectively. Double-stranded cDNA was purified with QIAQuick PCR extraction kit (Qiagen), and cDNA was eluted by elution buffer. Following the synthesis of second strand, end repair, addition of single A base, and adaptor ligation, cDNAs were connected with sequencing adaptors. The library concentration was measured by real-time PCR, and Agilent 2100 Bioanalyzer was used for profiling the distribution of insert size. Library constructions were sequenced by Illumina HiSeq 4000 based on the manufacturer’s instructions (Illumina), and were sequenced for 100 cycles. The HiSeq Control Software HCS (HiSeq Control Software) (v3.3) with RTA (Real-Time Analysis) (v2.7.3) was used to provide the management and execution of the HiSeq 4000 experiment runs.

Sequencing Data Analysis.

Total RNA-seq data analysis was performed as previously published (20). Images generated by HiSeq4000 were converted into nucleotide sequences by base calling and stored in FASTQ format utilizing Illumina package bcl2faseq (v2.16.0.10). To filter the dirty reads, which contain adaptors, unknown or low Phred quality-scored bases were obtained from raw reads, and clean reads were generated. Clean reads were mapped to reference UCSC hg19 genome and gene sequences using Tophat2 (v2.1.0). No more than five bases mismatched reads were allowed in the alignment. To annotate gene expression, fragments per kb per million reads values of each gene were calculated using Cufflinks package (v2.2.1). DEGs (differentially expressed genes), fold-change analysis was performed by Agilent GeneSpring (v7.3). Hierarchical clustering analysis for expression pattern was performed using Agilent GeneSpring (v7.3). Functional enrichment analysis was done using Gene Ontology functional classification system (geneontology.org/) and DAVID website (https://david.ncifcrf.gov/) as well. Raw reads and data are accessible in the Gene Expression Omnibus (GEO) (accession nos. GSE 151270 and GSE 147792).

SNV Analysis.

We performed SNV calls using the SAMtools (v.1.3) and bcftools (v.1.3) with bam files from two samples (WT and 5xFAD cortex). Each of bam files was indexed with SAMtools, and SAMtools mpileup was performed with UCSC mm10 genome as reference. To obtain variant calling format (VCF) files for specific target genes, option –r was used to set the position during mpileup. To perform variant calling for five PLC genes, five specific regions were set (in detail, chr2:160,731,310 to chr2:160,761,923 for PLCγ1, chr2:134,786,164 to chr2:134,988,192 for PLCβ1, chr2:118,728,319 to chr2:118,709,202 for PLCβ2, chr2:135,741,830 to chr:136,013,068 for PLCβ4, and chr19:6,969,643 to chr:6,954,299 for PLCβ3), and each of five VCF files were generated and merged for further analysis. After mpileup, VCF files were generated with vcfutils.pl varFilter using option –D100. Functional annotation of each variant was performed with ANNOVAR software.

SpliceAI Analysis.

The human gene, PLCγ1, was analyzed by SpliceAI to predict RNA splicing types on the genetic variations (4). SpliceAI, a residual convolutional neural network, predicts the probability of RNA splicing taking place in either donor or acceptor given DNA sequences. SpliceAI trained the complex patterns of RNA splicing from DNA sequences of 10 kbps around nucleotides of interests from the GENCODE (Genome ENCyclopedia Of DNA Elements) V24 databases (GRCh38/hg38) including more than 10,000 genes in chromosomes 2, 4, 6, 8, 10 to 22, X, and Y. We conducted the analyses using the pretrained SpliceAI, downloaded at https://github.com/Illumina/SpliceAI. We computed the probability of RNA splicing types (probabilities of donor, acceptor, or neither) for the alternate nucleotides on the starting locus of every exon of the gene, PLCγ1, based on the reference sequences (GRCh38/hg38). Then, the probabilities of the variant being splicing altering on acceptor and donor are computed by

where $$ a_{ref},\hspace{1em}\hspace{1em}{a}_{alt},\hspace{1em}\hspace{1em}{d}_{ref}$$, and $$ d_{alt}$$ denote splice acceptor and splice donor probabilities of the reference and alternate nucleotide, respectively. Each $$ \mathrm{\Delta }Score$$ is defined as the maximum of difference between the reference and alternate probabilities on splicing sites. To be specific, the gain score is denoted by the maximum disparity of alternate from the reference probability, which means the probability of corresponding splicing increases by the gain score. Contrary to the gain score, the loss score demonstrates the maximum difference of reference from alternate probability inferring the probability of the splicing decreased by the loss score. SpliceAI also provided the $$ \mathrm{\Delta }Postion$$ of acceptor gain, acceptor loss, donor gain, and donor loss, which indicate location where splicing changes with the probability relative to the position of interest. The positive delta position is downstream of the variant, whereas negative values are upstream.

cDNA Synthesis.

The cDNA synthesis for mature mRNA was performed with PrimeScript first-strand cDNA synthesis kit (Takara) according to the manufacturer’s protocol. To synthesize cDNA, we used 1 µgof total RNA from each sample, and added the reaction mixture containing the Oligo dT Primer (2.5 µM), dNTP mixture (0.5 mM), 5× PrimeScript buffer, RNase Inhibitor (20 U) and PrimeScript RTase (200 U). The synthesized cDNA was heat incubated at 42 °C for 60 min and 95 °C for 5 min, amplified with a thermal cycler dice. RNA templates were reverse transcribed into cDNA, which corresponds to pre-mRNA transcript using a High-Capacity reverse transcription kit (Applied Biosystems) according to the protocol described. Following the reference component suggested, each reaction was composed of 10× RT buffer, 25× dNTP Mix (4 mM), 10× RT Random Primers, MultiScribe Reverse Transcriptase (50 U), RNase inhibitor, and adjusted total volume with nuclease-free water. Then, template RNA sample with equal volume was added to set the total reaction volume. The reverse transcription was performed as suggested optimal thermal cycler condition, heat incubated at 25 °C for 10 min, 37 °C for 120 min, and 85 °C for 5 min, amplified with a thermal cycler dice.

RT-PCR.

PCR was performed using the Phusion-HF DNA polymerase kit (NEB). Primers employed were PLCγ1 Ex27-29 forward, 5′-GAT​TGG​CAC​AGA​ACG​TGC​TT-3′, reverse, 5′-CTC​CGT​CTT​TTG​CTT​GGT​GC-3′; PLCγ1 Ex27-28 forward, 5′-GAT​TGG​CAC​AGA​ACG​TGC​TT-3′, reverse, 5′-GCA​AAT​GAC​ACA​GGG​TTC​CA-3′; PLCγ1 Ex28-29 forward, 5′-GCA​GAT​GAA​CCA​GGC​CC-3′, reverse, 5′-CTC​CGT​CTT​TTG​CTT​GGT​GC-3′; PLCγ1 Ex26-30 forward, 5′-CTG​AGG​GGA​AGA​TGA​TGG​A-3′, reverse, 5′-CAG​GCC​TTT​TAC​TGG​GAA​AG-3′; GAPDH forward, 5′-CAT​GGC​CTT​CCG​TGT​TCC​TA-3′, reverse, 5′-GCG​GCA​CGT​CAG​ATC​CA-3′; PLCγ1 Ex27-32 forward, 5′-GAT​TGG​CAC​AGA​ACG​TGC​TT-3′, reverse 5′-CAA​ATG​GCT​GCT​GGT​ATC​TG-3′. Each reaction mixture contains 5× phusion HF buffer, dNTP (0.2 mM), primer forward and reverse (0.5 μM), Phusion DNA polymerase (1 U), and template cDNA, and nuclease-free water was added. The amplification was performed by the recommended condition. PCR products were detected through the ChemiDoc (Bio-Rad) imaging system.

RT-qPCR.

Total RNA was isolated from mouse whole blood using the TRIzol method. The cDNA synthesis for RT-qPCR was performed by High-Capacity reverse transcription kit (Applied Biosystems). Primers for RT-qPCR employed were PLCγ1 coding forward, 5′-GGT​GAC​CTC​AGT​CCT​TTC​AG-3′, reverse, 5′- GAA​ATC​TTC​AAA​TGG​CTG​CTG-3′. RT-qPCR amplification was performed for 10 min for initial denaturation, followed by 45 cycles with 10 s at 95 °C for denaturation, 30 s at 60 °C for annealing, and 1 min at 72 °C for extension. PCR products were analyzed by LightCycler 480 (Roche).

Western Blot and Antibodies.

Electrophoresis was performed with 8% polyacrylamide gels, and proteins were then electrotransferred to nitrocellulose membranes. Membranes were immunoblotted with mouse monoclonal anti-PLCγ1 antibody (B16-5), generated as described previously (49), and mouse monoclonal anti-PLCβ1 (sc-5291, Santa Cruz), rabbit polyclonal anti-PLCβ3 (sc-403, Santa Cruz), rabbit polyclonal anti-PLCβ4 (sc-404, Santa Cruz), and mouse monoclonal anti-β-actin (GT5512, GeneTex) antibodies in Tween 20/Tris-buffered saline containing 5% skim milk. After incubation with the appropriate peroxidase-conjugated secondary antibody, proteins were detected with the enhanced chemiluminescence system (ECL system, Amersham).

ChIP-Seq Data.

H3K27ac binding profiles throughout the mouse forebrain were adapted from the data generated by GSE52386 data. H3K27ac ChIP-seq data are available in GEO (https://www.ncbi.nlm.nih.gov/geo/) (50).