Age-associated bimodal transcriptional drift reduces intergenic disparities in transcription

Measurement of gene expression levels in young and old samples using SiNG-PCRseq

For gene expression analysis, we used a set of primer pairs that can amplify 285 exonic amplicon sequences of 175 epi-driver genes [7]. Using SiNG-PCRseq, we determined the amount of target transcript levels in the skeletal muscle, brain (striatum), and splenic CD4+ T cells collected from wild type young (wy, 2 months old) and old (wo, 16–19 months old) mice and age-matched HD mice (hy and ho, respectively). In running multiplexed SiNG-PCR, rat genomic DNA (gDNA) was used as a spike-in for accurate quantitation of mouse cDNA sequences relative to matched rat gDNA sequences [7]. Each of the primer pairs was designed to carry one (or >2 less frequently) nucleotide variation between the DNAs of different origin species to best achieve their equal amplification and to set distinguishable flags between the cDNA-derived reads and the rat gDNA reads after deep sequencing. Splenic CD4+ T cells were pooled from three or four mice and split into four replicates before use in multiplex PCR. Brain and muscle samples were analyzed from individual animals. We obtained 6.6 × 10⁶ (muscle), 8.5 × 10⁶ (brain), and 29.6 × 10⁶ (T cells) reads on average per index after sequencing about 80% of which were mapped to our amplicon-built reference (Supplementary Data File 1).

The expression level of each exonic sequence was measured by calculating the M/R ratio, the ratio of the mouse read counts relative to the rat’s counts (Supplementary Data File 1). Principal component analysis (PCA; Supplementary Fig. S1) showed that the muscle samples did not differ much, whereas brain tissue showed a moderate distinction between ages and also between the genotypes. The young-old difference was the most prominent in T cells. The young T cells of wild-type and HD mice grouped together but with age the two genotypes diverged greatly and became distinct from each other.

Age-associated transcriptional alteration in splenic T lymphocytes

We used a scatter plot to compare the expression levels between young versus old samples. Each of the tissue samples showed a high similarity between wild-type young versus HD young (wy-hy; Fig. 1A-C). Therefore, we took wy as the reference for the age-related changes occurring in both wild-type and HD samples. In muscles, the young to old comparisons (wy-wo and wy-ho) showed no substantial differences, as shown by their strong correlations (R² > 0.976). In the brain, however, a relatively large proportion of genes deviated to a degree from the regression curve for old HD (wy-ho, R² = 0.933). These alterations of epi-driver gene expression levels are potentially related to pathogenic events in the old HD brain and appeared to occur randomly as judged by the almost undisturbed plot slope (m = 0.9879; Fig. 1B).

Figure 1. Age-associated changes in quantitative transcriptome structure of splenic T cells. (A-C) Scatter plots comparing the epi-driver gene expression levels between wild-type young (wy) and HD young (hy) samples (wy-hy), between wild-type young (wy) versus wild-type old (wo) samples (wy-wo), and between wild-type young versus HD old (ho) samples (wy-ho) in skeletal muscle (A), brain tissue (B), and splenic CD4+ T cells (C). Blue and gray lines denote a regression curve and the reference slope (m = 1.0), respectively. (D) Weighted root mean square deviation (wRMSD) analysis for assessment of the extent of difference in gene expression level between young and old cells. (E-F) Comparison of the correlations (coefficient of determination, R²; (D) and the slopes (E) of the scatter plots. Mu, muscle; Br, brain; Tc, T cells. Dotted green lines indicate trend lines for T cell samples.

The age-related changes to epi-driver gene expression were distinctive in the splenic T cells, as shown by the markedly lower plot slopes compared to muscle and brain (Fig. 1C). The young-old difference in expression levels was quantified by calculating the weighted root mean square deviation (wRMSD) of old samples from young samples [14–16]. The deviation was the highest in old HD T cells (wy-ho) and next from the old wild-type T cells (wy-wo), with the wRMSD values 79.3 and 15.4 fold larger, respectively, compared with those of young muscles (Fig. 1D). In the T cell samples, the plot correlations and slopes also decreased in the order hy, wo, and ho (Fig. 1E and 1F, respectively). While the wRMSD, correlation coefficients, and slopes could all be used to quantify gene expression differences, the former two cannot express the change of direction (increase or decrease) in gene expression; thus, we chose the slope (m) as the tool to express the age-associated alteration of gene expression, A₃GE, which we defined as

‘$\sqrt{{(1-m)}^2}$‘.

Alteration in the quantitative transcriptome structure as a function of age in T cells

Figure 1 demonstrated that A₃GE is clearly identifiable in the CD4+ T cells. The exacerbated A₃GE in old HD T cells might result from accelerated aging in HD mice. To test whether aging progresses more rapidly in HD T cells, we assessed the biological ages of old T cells by measuring their telomere lengths [17,18]. The copy number ratio (T/S ratio) of the telomeric sequence relative to a single-copy gene sequence (36B4) [19,20] was first determined in T cells from 2, 20, 25, and 28 month old wild-type mice. The T/S ratio was found to decrease with age (Fig. 2A) at a rate of 1.83% per month, and there was a strong inverse correlation with age (R² = 0.9841; Fig. 2B).

Figure 2. A model for the alteration of gene expression levels as a function of age in T cells. (A) Determination of the copy-number ratio (T/S ratio) of the telomeric sequence to a single-copy gene sequence (S6B4) in splenic CD4+ T cells of wild-type mice with different ages (2, 20, 25, and 28 months). Error bars, standard deviation. (B) Determination of the biological age of the ho T cells. Regression curve (dotted blue line) was derived by the T/S ratios of wild-type T cells in A. The T/S ratio of the 19 month old wo sample corresponds to approximately 0.68 (left red triangle). The relative T/S ratio of the ho T cells was 69% of the wo T/S ratio; the T/S ratio of the ho is 0.47 (0.68 x 0.69), which corresponds to about 30 months of age in wild-type mice (right red triangle). (C) Relative T/S ratio of the old HD (ho) T cells to the old wild-type (wo) T cells. In A and B, quantitative real-time PCR was performed to obtain C_t values for the telomere repeats and S6B4 sequence before calculating the copy-number difference. Error bars, standard deviation. (D) A regression model of the age-associated alteration of gene expression (A₃GE defined as $\sqrt{{\left(1-m\right)}^2}$, where m is a slope) in splenic T cells. A₃GE increases at a rate of approximately 2.87% per month on average in wild-type mice.

Comparing old HD (16-19 months) and old wild-type (19 months) T cells, the T/S value of the former was 69% of the latter (Fig. 2C). This result provides evidence for accelerated aging in HD T cells. As shown in Fig. 2B, the relative T/S value of 19 month old, wild-type T cells corresponded to 0.68 on the linear regression if the T/S ratio of two month old T cells was set to 1.0. Hence, the T/S value of old HD T cells would be 0.47 (0.68 x 69%), which implies that they were comparable to 30 month old wild-type T cells. Hence, in terms of telomere length, the 16–19 month old HD T cells aged approximately 1.5-fold more than the age-matched wild-type T cells. Meanwhile, seven month old wild-type (wm) T cells had gene expression levels comparable to the young T cells, exhibiting no noticeable sign of A₃GE (Supplementary Fig. S2). This suggests that A₃GE does not start at the beginning of life but is somehow initiated during middle of life. Fig. 2D shows that A₃GE values of wild-type T cells of different ages and age-adjusted old HD T cells exhibit a linear relationship with age (R² = 0.9676). In this regression model, the regression coefficient was 0.0287, implying that A₃GE increases at a rate of 2.87% per month on average and reaches 50% at approximately 21 months of age in wild-type mice.

Categorical analysis of epi-driver genes for age-associated alteration of gene expression

We classified the epi-driver genes into eight different categories according to their roles in the modification of chromatin and compared their expression levels between young and old T cells. In all categories, a decrease of the plot slope was consistently observed in the old wild-type T cells, and it was even more severe in old HD T cells (Fig. 3A). When the A₃GE in individual epi-driver categories was compared a relatively uniform A₃GE was observed across the categories (Fig. 3B). An exception was the acetylation category in which the HD genotype dependent A₃GE intensification was undetectable. This probably indicates transcriptional stability of acetylation category genes with age. In the muscle and brain samples, however, A₃GE was not observed (see also Supplementary Fig. S3A and S3B). The slope and A₃GE values for individual categories are summarized in the Supplementary Fig. S4. The extent of A₃GE was so marked that the correlation among the genes was completely lost in some categories including the arginine-methylation (R² = 0.172), PRG (0.019), and ubiquitination categories (0.013; Fig. 3A). wRMSD calculations showed that the acetylation category was least affected whereas the PRC-regulated gene (PRG) category deviated with age the most (Fig. 3C), which is similar to the result of A₃GE measurements in Fig. 3B.

Figure 3. Comparison of A₃GE in different epi-driver gene categories between old wild-type and HD T cells. (A) Difference in the extent of A₃GE between wo and ho T cells. Two scatter plots which compare expression levels between young (wy-hy) and old (wy-wo; blue) and between young and HD old (wy-ho; orange) are merged. The slope (m) of the regression curve and the correlation value (R², coefficient of determination) are shown. Each slope represents the extent of A₃GE in the corresponding T cells. PRC, Polycomb-repressive complex; PRG, PRC-regulated genes. (B-C) Comparisons of the levels of A₃GE (B) and weighted root mean square deviation (wRSMD) (C) in individual gene expression categories in young (wy and hy) and old (wo and ho) T cells. A₃GE was minimal in muscle and brain samples which are displayed in the smaller plots. A₃GE, age-associated alteration of gene expression.

In addition, we examined a further gene category that consisted of aging-related genes [7]. It showed a similar A₃GE in the T cells (Supplementary Fig. S5A), indicating that A₃GE was not restricted to the epi-driver genes. Another group of genes, which were selected from the literature for their relatively constant expression in various tissues, showed only a weak A₃GE in T cells (Supplementary Fig. S5B). These observations suggest that the extent of A₃GE is determined by the intrinsic stringency of gene expression regulation. We also examined MACS-purified splenic macrophages from old wild-type mice. SiNG-PCRseq results showed that epi-driver gene expression levels did not noticeably change with age in these cells (m = 0.988 and R² = 0.907; Supplementary Fig. S6), which refutes the possibility that A₃GE affects all blood cells similarly.

Transcriptional drift − Age-related reduction of gene-to-gene disparity of transcriptional activity

Another interesting feature of A₃GE was that, as shown by the lowered plot slopes in Fig. 3A, those genes that were weakly expressed in young samples tended to be increasingly transcribed in old samples, whereas those genes which were relatively highly expressed in young cells tended to decline in old cells. Representative trends are presented for the acetylation and PRC category plots (Fig. 4A). The median values of relative (M/R) gene expression levels for all of the epi-driver genes ranged from 1.97 to 2.55 and did not greatly differ among the four T cell sample groups (Fig. 4B). However, interestingly, the variance of gene expression levels of them decreased to a half in the old wild-type T cells and to a quarter in old HD T cells compared to the young cells (Fig. 4C). In contrast, the muscle and brain samples showed no apparent change in variances. The reduced variance in old T cells demonstrates that the gene expression system in old T cells tends to narrow the differences in expression level between strongly and weakly transcribed genes.

Figure 4. Opposite transcriptional drift of genes with weak and strong transcriptional strengths. (A) Representative pattern of an age-related expression changes in the acetylation (upper) and PRC (Polycomb repressive complex; lower) category genes. The expression levels of individual epi-driver genes in young T cells (blue) are paired with corresponding genes in old T cells (orange). The thick arrows denote the direction of expression-level change in low expressing genes (left) and high expressing genes (right) with age. (B-C) Measurement of the median (B) and variance (C) of gene expression levels in the brain (B) and T cells (T) relative to those in the muscle samples. (D-E) Opposite directions of age-associated alteration of gene expression (A₃GE) between the lowly (LEGs, n = 50; D) and highly expressed genes (HEGs, n = 50); (E). Two scatter plots compare expression levels between young (wy-hy) and old (wy-wo; blue) and between young and HD old (wy-ho; orange) are merged. Box plots on the right show the opposing directions of change and also the statistical significance of the A₃GE differences between the samples. (F) A₃GE in the transcriptomes of CD4+ T cells. RNA-seq was performed using mRNAs of MACS-purified splenic CD4+ T cells pooled from three different mice of two or 20 months of age. Five hundred top and 500 bottom ranked genes were selected, HEGs (FPKM values ranging from 144.1 to 7,444.7) and LEGs (from 0.0041 - 0.035) respectively, after transcriptomes were sorted by expression level for FPKM levels of young T cells. Expression levels of old wild-type T cells relative to young were calculated using a log₂ scale. P-values in D-F, paired sample t-test.

In order to validate the general nature of this convergence of expression, we examined A₃GE in two groups of genes, top ranked (highly expressed genes, HEGs; n=50) and bottom ranked (lowly expressed genes, LEGs; n = 50) determined after sorting the genes of the young wild-type T cells by expression level. We observed that the LEGs mostly showed increased levels in the old T cells and, conversely, the HEGs showed decreased levels (Fig. 4D and 4E, respectively). This tendency to plateau was far more obvious for the LEGs; the fold change in expression level was 2.7 (wo) to 5.3 (ho) on a log₂ scale in the LEGs while it was 0.5 (wo) to 0.9 (ho) in the HEGs. Box plots demonstrate the opposing directions of expression change and the significant differences between LEGs and HEGs. The increasing transcription of LEGs in old T cells was not due to RNA polymerase II activity was decreased in old T cells (Supplementary Fig. S7), agreeing with a previous report of the downregulation of RNA splicing and processing genes in old cells [21].

We further examined the whole transcripts for A₃GE in old wild-type mice. RNA-seq was performed using mRNAs from splenic CD4+ T cells and peripheral blood mononuclear cells (PBMCs). Fifty-four million reads on average per sample (82.5%) were uniquely mapped onto the mm10 reference (Supplementary Data File 1). Comparison of the expression levels of LEGs and HEGs (n = 500 each) revealed a similar trend of transcriptional drift between two and 20 month old mice and displaying a significant difference between the HEGs and LEGs (p = 3.4 x 10^-60, paired sample t-test; Fig. 4F). In contrast to the CD4+ T cells, PBMCs showed no such a change for HEGs and only a slight increase in expression level for LEGs (Supplementary Fig. S8). These data indicate that this pattern of transcriptional drift, of all the cell types tested, is unique to CD4+ T cells.

Ethics statement

This study was carried out in strict accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Livestock Research Institute of Korea. The protocol was approved by the Committee on the Ethics of Animal Experiments of the Korea Research Institute of Bioscience and Biotechnology (KRIBB) and by Hanyang Institutional Animal Care and Use Committee (HY-IACUC-09-017).

HD model mice and tissue preparation

YAC128 mice (FVB-Tg(YAC128)53Hay/J; Jackson Laboratory) which contains 128 CAG repeats of the human HTT gene [47] and wild-type littermates were kept in a 12 hour light-dark cycle with unrestricted access to food and water. Transgene in YAC128 was detected by PCR genotyping using primers 5’-CCGCTCAGGTTCTGCTTTTA-3’ and 5’-TGGAAGGACTTGAGGGACTC-3’. We used a total 7 of female wild-type littermate mice at 16-19 months (n=4) and 2 months (n=3) of age, and a total 7 of female YAC128 HD mice at 16-19 months (n=3) and 2 months (n=4).

The method for isolation of splenic CD4+ T cells and macrophages was described in detail elsewhere [7]. Brains were dissected to obtain striatal tissues and immediately after isolation, samples were frozen in dry ice and stored at -80°C deep-freezer. For skeletal muscles, tissues were obtained from hind legs and homogenized using Biomasher II (Biomasher) prior to total RNA extraction. cDNAs were generated from 1 μg of total RNA by incubating with Super-Script III (Invitrogen), oligo(dT)₂₀ (Invitrogen), and random-hexamers (Invitrogen) at 50°C for 1 hour followed by enzyme inactivation at 75°C for 10 min.

Multiplex PCR

Detection of nucleotide inter-species variation (ISVs) present in target genes between mouse mRNA and rat genomic DNA sequences through sequence alignment via BLAST, and primer design for the ISV-containing homologous sequence stretches are described in detail elsewhere [7]. In total, 285 primer pairs targeting 175 genes were split into 10 groups composed of 28~29 primer pairs and used in multiplexed PCR with 15 ng of mouse cDNA and 2 ng of rat gDNA as templates in following conditions; 95°C for 15 min, then at 95°C/20 sec, at 57°C /40 sec, at 72°C/1 min for 45 cycles and 72°C for 5 min. All experiments were duplicated for each sample.

Amplicon library production and next generation sequencing

For construction of amplicon sequencing libraries for Illumina sequencing platform, PCR products of each sample generated by 10 groups of primer sets were mixed together and purified using Expin purification kit (GeneAll). The amplicon ends were phosphorylated with T4 polynucleotide kinase (NEB) by incubating at 37°C for 30 min, and were ligated with NGS adapters using T4 DNA ligase (Solgent) at 37°C for 2 hour. Finally, index PCR was conducted to complete Illumina sequencing libraries. The index PCR conditions are as follows; after 98°C for 15 min, proceeded to 20 cycles of 98°C/20 sec, 68°C/30 sec, and 72°C/1 min, followed by a final extension at 72°C for 5 min. The amplified libraries were purified and sequenced using HiSeq 2500. The barcoded sequencing libraries were pooled with equivalent amount and subjected to a multiple parallel sequence using Illumina HiSeq 2500 platform.

wRMSD analysis

In order to quantify young-old difference in expression levels, we calculated the root mean square deviation (RMSD) of old samples from young samples. To minimize the bias that could result from a measurement error of gene expression profile of a group with a low coefficient of variation (CV, the deviation of each gene expression level from the mean) was weighted with the CV of the gene in the group [14–16]. The wRMSD was obtained with the equation,

$\mathrm{w}\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{D}=\sqrt{\sum {w_i\bullet ({Em}_i-E_i)}^2}$

where w$i$, Em$i$, and E$i$ are the values obtained from the $i$th gene of interest for the weight of the mean square deviation of the gene expression, the reference expression level (e.g., the mean expression level in the group), and the expression level of the gene, respectively. The weight (w$i$) is the proportion of the CV for the expression level of the $i$th gene to the sum of CV for those of all genes in the group and was obtained with the equation,

$w_i=\frac{{CV}_i}{\sum {CV}_i}$.

Telomere length assay by qPCR

Genomic DNA was extracted from MACS-purified splenic T cells using DNeasy Blood & Tissue Kits (QIAGEN). For quantitative real-time PCR-mediated Telomere length assay [20], acidic ribosomal phosphoprotein P0 (36B4) single-copy gene was used as a standard. Primers used were 5′-CGGTTTGTTTGGGTTTGGGTTTGGGTTTGGGTTTGGGTT-3′ and 5′-GGCTTGCCTTACCCTTACCCTTACCCTTACCCTTACCCT-3’ for telomeric sequence, and 5'-ACTGGTCTAGGACCCGAGAAG-3' and 5'-TCAATGGTGCCTCTGGAGATT-3' for 36B4 gene sequence. Quantitative real-time PCR was performed using 10 ng of genomic DNA and 2X Power SYBR Green PCR Master Mix (ABI) in QuantStudio 3 Real-Time PCR System (ABI). Real-time PCR condition was set as 95°C /15 sec, 55°C /1 min and 72°C /1 min. Telomere to single copy gene (T/S) ratio is calculated as

${{[2}^{\mathrm{C}\mathrm{t}\left(\mathrm{t}\mathrm{e}\mathrm{l}\mathrm{o}\mathrm{m}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{s}\right)}/{\mathrm{ }2}^{\mathrm{C}\mathrm{t}\left(36\mathrm{B}4\right)}]}^{-1}=2^{-\mathrm{Ä}\mathrm{C}\mathrm{t}}$.

Target gene expression estimation

Prior to the aligning of sequence reads, a custom reference genome was generated by gathering mouse mRNA sequences (mm10) of 260 target genes, and the reference was indexed by ‘bowtie2-build’ [60]. The raw sequencing reads were groomed to remove unwanted sequences and low quality bases using ‘Trim_galore (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/’. The trimmed reads were mapped on the custom genome using ‘bowtie2’ with default parameters.

For quantification of mouse and rat read counts, sequence variation analyses were performed according to ‘The genome analysis toolkit (GATK)’ pipeline [61] with some modifications. The mapped reads were sorted by ‘Picard’ (https://broadinstitute.github.io/picard/), and sequence variations (SNVs) were called using ‘UnifiedGenotyper (GATK)’ with the down-sampling option disabled. Due to the nature of amplicon sequencing, duplicated reads were preserved throughout following analyses. Since INDEL was not included for targets, the INDEL realignment step was skipped. The resulting VCF files were filtered to obtain read count information only at pre-defined the SNVs sites. The number of reference (mouse) and alternative (rat) sequences was counted at the filtered SNV sites using home-brew bash scripts. The count data were corrected for deviations in the fractional quantities as previously described [8]. Finally, gene expression levels were estimated by calculating relative read counts of mouse to rat (M/R ratio). To draw PCA plots, raw count data for all samples were combined and normalized by ‘DESeq’ [62]. Using ‘plotPCA’ function, the plots were generated.

RNA-seq and data analysis

Total RNAs were extracted from T and blood cells of each age group using RNeasy Plus Mini Kit (Qiagen). Dynabeads mRNA DIRECT (Thermo) was used to isolate poly-A tailed RNAs from 2 μg of total RNAs. Each step was conducted according to the manufacturer's instruction. Isolated poly-A tailed RNAs were used to construct an RNAseq library using NEBNext Ultra RNA library prep kit for illumina (NEB). Isolated poly-A tailed RNAs were incubated at 94°C for 15 min for fragmentation. First strand cDNA was synthesized with fragmented RNAs and ProtoScript II Reverse Transcriptase. Second strand cDNA was synthesized with Second Strand Synthesis Enzyme Mix in the kit and purified. After end repair of the synthesized double-stranded DNAs using NEBNext End Prep Enzyme Mix in the condition of 20°C for 30 min and 65°C for 30 min, NEBNext Adaptor was added into the 3’ end adenylated DNA with Blunt/TA Ligase Master Mix and the mixture was incubated at 20°C for 15 min. DNA samples were enriched by 15 cycles of PCR reaction consisting of the purified ligates, NEBNext Q5 Hot Start HiFi PCR Master Mix 2X, and index primer. Each purification was conducted using Ampure XP beads (Beckman). Enriched DNA fragments were purified and sequenced using NextSeq 500.

Raw sequencing reads were trimmed to remove Illumina adapter sequences and low quality bases using ‘Trim_galore (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/’, and the trimmed reads were aligned on the reference transcriptome (mm10) using STAR: ultrafast universal RNA-seq aligner [63]. For estimation of gene expression levels, we followed ‘cufflinks pipeline’ [64]. Transcriptome was assembled by ‘cufflinks’ and ‘cuffmerge’, and gene expression levels were estimated and by ‘cuffquant’, and finally, a matrix of gene expression levels in RPKM was generated by ‘cuffnorm’.