Whole blood transcriptome analysis for age- and gender-specific gene expression profiling in Japanese individuals

Abstract

Whole blood transcriptome analysis is a valuable approachin medical research, primarily due to the ease of sample collection and the richness of the information obtained. Since the expression profile of individual genes in the analysis is influenced by medical traits and demographic attributes such as age and gender, there has been a growing demand for a comprehensive database for blood transcriptome analysis. Here, we performed whole blood RNA sequencing (RNA-seq) analysis on 576 participants stratified by age (20–30s and 60–70s) and gender from cohorts of the Tohoku Medical Megabank (TMM). A part of female segment included pregnant women. We did not exclude the globin gene family in our RNA-seq study, which enabled us to identify instances of hereditary persistence of fetal hemoglobin based on the HBG1 and HBG2 expression information. Comparing stratified populations allowed us to identify groups of genes associated with age-related changes and gender differences. We also found that the immune response status, particularly measured by neutrophil-to-lymphocyte ratio (NLR), strongly influences the diversity of individual gene expression profiles in whole blood transcriptome analysis. This stratification has resulted in a data set that will be highly beneficial for future whole blood transcriptome analysis in the Japanese population.

Graphical Abstract

age and gender difference, globin, HPFH, RNA-seq, whole blood transcriptome

Abbreviations

 
• AHSP

  alpha hemoglobin stabilizing protein

 
• BirThree

  birth and three-generation

 
• CommCohort

  community-based cohort

 
• DEG

  differentially expressed gene

 
• DV₂₀₀

  percentage of RNA greater than 200 nucleotides in length

 
• eQTL

  expression quantitative trait locus

 
• FDR

  false discovery rate

 
• GTEx

  Genotype-Tissue Expression

 
• HMBS

  hydroxymethylbilane synthase

 
• HPFH

  hereditary persistence of fetal hemoglobin

 
• IMM

  Iwate Tohoku Medical Megabank Organization

 
• IRB

  institutional review board

 
• NLR

  neutrophil-to-lymphocyte ratio

 
• PCA

  principal component analysis

 
• PC1

  first principal component

 
• PC2

  second principal component

 
• RIN

  RNA integrity number

 
• RNA-seq

  RNA sequencing

 
• rRNA

  ribosomal RNA

 
• TMM

  Tohoku Medical Megabank

 
• ToMMo

  Tohoku Medical Megabank Organization

 
• TPM

  transcripts per million

 
• WGS

  whole genome sequencing

In March 2011, a huge earthquake and tsunami devastated northeastern coast of Japan (1). In an effort to innovatively regenerate the tsunami-damaged regions, we established Tohoku Medical Megabank (TMM) Project (2, 3). This project aims to realize a next generation medicine through prospective genome cohort studies. To facilitate the TMM Project, Tohoku Medical Megabank Organization (ToMMo) was established at Tohoku University, and Iwate Tohoku Medical Megabank Organization (IMM) was established at Iwate Medical University. Both ToMMo and IMM have collaborated to create two types of cohorts (4). One is the community-based cohort (CommCohort Study), in which a total 84,073 participants were recruited (5). The other is the birth and three-generation cohort (BirThree Cohort Study), which included pregnant mothers, their partners, parents, newborn babies and their elder siblings (6). In the BirThree Cohort Study, a total of 73,500 participants were recruited from the obstetrics hospitals and clinics in the Miyagi Prefecture from 22,493 families.

Employing these two prospective cohorts strategically, we established TMM biobank, which houses biospecimens, health and clinical information, as well as genome-omics data (7). We carried out comprehensive health examination during the initial survey and have been conducting similar health examinations in the follow-up survey every five years. Importantly, we made up our mind to establish the TMM biobank as an integrated biobank. For this purpose, we have set up an analytical center within the TMM biobank that carries out large-scale genome and omics analyses (8). Through the analytical center, we have developed a Japanese whole-genome reference panel based on the whole genome sequence data (9) and a custom SNP array specially designed for the Japanese population (10).

In our analytical center, we are conducting a range of omics analyses, including transcriptomics, epigenomics, proteomics, metabolomics and metagenomics (11–14). Our goal is to elucidate gene–environment interactions within the Japanese population and gather robust evidence that will facilitate the advancement of future medicine in Japan (15). Given the large-scale nature of the cohort samples we handle, we have strived to establish high-throughput procedures for omics analyses. Among the ongoing omics analyses, RNA sequencing (RNA-seq) or transcriptome analysis is particularly compelling. This method offers the ability to broadly survey and/or scrutinize gene expression profiles in a manner that whole genome sequencing (WGS) cannot provide (16). However, an inherent challenge of cohort transcriptome analysis lies in the limitation of available RNA sources.

Prospective genome cohorts are typically composed of healthy volunteers, which restrict the types of biological materials available for study. Therefore, we chose to use whole blood samples, as these are comparatively easy to collect from cohort participants. We hypothesized that transcriptome analysis of whole blood samples would provide extensive information reflective of individual health conditions (17). Indeed, transcriptome analysis of whole blood samples has been conducted globally (18–21). For instance, the Genotype-Tissue Expression (GTEx) Project has determined gene expression profiles for various tissues, including whole blood (22). The whole blood transcriptome data from GTEx have been utilized to estimate gene expression levels in other tissues (23). Similarly, whole blood transcriptome analysis was carried out in a case-control study targeting COVID-19 in Japan (24).

However, whole blood gene expression analysis targeting healthy populations in prospective cohorts has yet to be conducted. This attempt is particularly important as uncertainty remains about the variations in gene expression profiles in relation to age and gender in whole blood transcriptomes. Additionally, there is no consensus methodologically on whether it would be more advantageous to remove globin mRNAs from the whole blood transcriptome (25). In an effort to address these gaps, we embarked on a large-scale transcriptome analysis using whole blood samples. These samples were collected from 576 female and male as well as young and elder participants of the TMM Project, with the aim of generating fundamental data about gene expression in the whole blood of the Japanese population.

Materials and Methods

Study design and sample selection

This project was performed as a part of the prospective cohort study at the TMM, with the approval of the institutional review board (IRB) at Tohoku University. The samples used here were from the cohort participants, all of whom gave TMM their written consent. In order to conduct whole blood transcriptome analysis, 576 participants' samples were selected randomly, with adjustments made for the gender and age. The 576 whole blood samples include 120 females in age of 20 to 30 years, 172 females in 60 to 70 years, 109 males in 20 to 30 years and 175 males in 60 to 70 years. Of note, 14 of the 120 females in the age of 20 to 30 years were during pregnancy. Cohort participants suffered from cancers, high blood pressure, diabetes mellitus and hyperlipidemia were excluded.

Total RNA extraction from PAXgene® blood RNA tubes

As a part of baseline analysis, we collected from 4337 participants approximately 2.5 mL whole blood in PAXgene® Blood RNA Tubes (#762165, BD Biosciences) containing 6.9 mL of the additive and stored them at −80°C. Total RNA was extracted using PAXgene® Blood RNA kit (#762174, BD Biosciences) according to the manufacture protocol, which included DNase I treatment to remove trace amounts of DNA. Concentrations of RNA samples were quantified by NanoDrop 2000 (Thermo Fisher Scientific). Quality of RNA samples was evaluated by means of RIN (RNA integrity number) and DV₂₀₀ (percentage of RNA greater than 200 nucleotides in length) values using RNA ScreenTape and Reagents (#5067–5576 and #5067–5577, Agilent Technologies) by TapeStation 2200 (Agilent Technologies).

RNA-seq analysis

The extracted RNA samples were diluted to 11 ng/μl concentration in 96-well plates (20 μl or 220 ng/well). After the RNA was subjected to the removal of ribosomal RNA (rRNA) using rRNA Depletion Kit (#1000005953, MGI Tech), RNA libraries were prepared using a MGISP-960 high-throughput automated sample preparation system (MGI Tech). Through the RNA libraries, DNA libraries were constructed with MGIEasy RNA Directional Library Prep Set (#1000006386, MGI Tech) according to the manufacturer's protocol. The DNA libraries were sequenced with DNBSEQ-G400RS High-Throughput Sequencing Set (#1000016952 with FCL PE150 or #1000016995 with App-A FCL PE150, MGI Tech) using a DNBSEQ-G400 sequencer (MGI Tech).

Data analysis

Raw sequence reads were quality-controlled by removing low-quality bases using Trimmomatic (26) with default parameters. Trimmed reads were aligned to the human genome GRCh38 using STAR v2.5.3a (27) based on GENCODE release 28 gene annotations (28). Read count and expression level of individual genes estimated as TPM (Transcript per Million) (29) were calculated using RSEM v1.3.0 (30) with default parameters.

Expression levels of individual genes excluding the effect of globin genes were calculated based on the read count matrix, which removes components corresponding to 12 genes of the globin gene family (ENSG00000206172 HBA1, ENSG00000188536 HBA2, ENSG00000244734 HBB, ENSG00000229988 HBBP1, ENSG00000223609 HBD, ENSG00000213931 HBE1, ENSG00000213934 HBG1, ENSG00000196565 HBG2, ENSG00000206177 HBM, ENSG00000086506 HBQ1, ENSG00000130656 HBZ, ENSG00000206178 HBZP1), using in-house Python script (contents are available upon request).

Principal component analysis (PCA) was performed using expression profiles of the genes that were found to be expressed in all samples by using the prcomp function in the R environment. Differentially expressed genes (DEGs) were determined using the DESeq2 package (31) in the R environment with a false discovery rate (FDR) cutoff of 1%.

Fig. 1

Whole blood transcriptome analysis toward the establishment of gene expression reference panel for Japanese population. (A) Overview of the study design. (B) Schematic diagram RNA-seq analysis flow.

Data availability

The gene expression statistics are available at the jMorp (https://jmorp.megabank.tohoku.ac.jp/) website via a web interface (15). The raw sequence data are under controlled access, as the sequence data contain information under restriction of the participants’ consent. The sequence data are available upon request after approval of IRB and Materials and Information Distribution Review Committee.

Results

Whole blood transcriptome analysis toward establishment of gene expression reference panel for Japanese population

To elucidate an age and gender-specific gene expression profile of the Japanese population using whole blood, we performed RNA-seq analysis of 576 blood samples, including 120 females in age of 20–30s, 172 females in 60–70s, 109 males in 20–30s and 175 males in 60–70s derived from the TMM biobank (Fig. 1A, left panel). Of note, 14 females out of the 120 females in age of 20–30s were during pregnancy. This is due to the nature of our cohort-biobank in which we have recruited pregnant women into the BirThree Cohort Study. In this study, we selected samples randomly from whole available participants.

Total RNA samples of sufficient quantity and quality for RNA-seq experiments were obtained from all 576 blood samples taken by the PAXgene® Blood RNA Tubes, which were cryopreserved in the TMM biobank for an extended period (Supplementary Fig. S1-S2, Table 1). It should be noted that in whole blood RNA-seq cases, both rRNA and globin mRNAs are generally depleted experimentally before sequencing to improve the quality of transcriptome data. However, in this study, we decided to exclude only rRNAs, but to include globin mRNAs in the transcriptome experiment to gain insight into the population-level expression characteristics of the globin gene family (Fig. 1A, middle panel).

As it has been reported that in silico removal of globin mRNA genes from sequencing results improves the quality of the transcriptome data, we processed the transcriptome data and prepared gene expression profiles of non-globin genes by excluding the reads assigned to the globin gene family (Fig. 1A, right panel). We subsequently evaluated the gene expression data set separately for globin and non-globin genes.

Schematic diagram of RNA-seq data analysis and approach of excluding globin mRNAs during the data analysis is shown in Figure 1B. In this protocol, raw sequence reads were first quality controlled by removing low-quality bases using a software Trimmomatic, and trimmed reads were aligned to the human genome GRCh38. The read count and expression level of individual genes estimated as TPM were calculated using RSEM v1.3.0, and the expression levels were calculated based on the non-globin genes only read count matrix using in-house python script. PCA was performed based on the expression profile of genes found to be expressed in all samples.

Characteristics of globin gene expression

Taking advantage of our approach for the whole blood transcriptome analysis, we first examined expression characteristics of the globin family genes. We found significant level of both HBA1/HBA2 (α1/α2 globin) and HBB (β globin) gene expressions, which encode the subunits of adult-type hemoglobin, and this was confirmed in all 576 samples (Table 2). In this transcriptome analysis, globin gene content varies sample to sample, but occupies 18.9% to 57.1%, showing very good agreement with that reported 20.0% to 62.6% (32). When we assessed relationship between the expression levels of each globin gene by TPM, we found a strong correlation between expression levels of the HBA1 and HBB genes (Fig. 2A). This indicates that transcript levels of HBA1 and HBB genes are well coordinated in adult peripheral blood. This observation is also consistent with the global gene co-expression data provided by COXPRESdb (https://coxpresdb.jp/locus/?gene_id=3043).

Table 1

Cohort RNA samples extracted from PAXgene® Blood RNA Tubes

	Ave.		SD	Min		Max
RNA (ng/μL)	33.5	±	19.6	4.2	–	160.0
RNA (μg)	2.7	±	1.6	0.2	–	12.8
RIN	7.9	±	0.6	4.9	–	9.1
DV200	89.3	±	6.6	35.1	–	97.7

RIN, RNA Integrity Number; DV₂₀₀, percentage of RNA greater than 200 nt in length

Table 2

Expression statistics of globin genes

	Gene	Detection Number	Detection Rate [%]	Minimum TPM	Median TPM	Maximum TPM
α-Globin Locus	HBZ	0	0.00	0.00	5.76	46.83
	HBM	0	0.00	17.18	85.71	699.21
	HBA1	576	100	2.31 × 1005	6.34 × 1005	1.13 × 1006
	HBA2	576	100	2.59 × 1004	8.49 × 1004	1.69 × 1005
	HBQ	0	0.00	8.43	31.27	85.35
β-Globin Locus	HBE	0	0.00	0.00	0.09	5.12
	HBG2	7	1.22	42.69	3.49 × 1003	5.28 × 1005
	HBG1	2	0.35	11.70	1.16 × 1003	1.31 × 1005
	HBD	0	0.00	5.26	37.34	574.46
	HBB	576	100	5.08 × 1004	1.60 × 1005	3.56 × 1005

Column ‘Detection Number’ means the number of samples in which read count composition of the gene in question exceeds 1%.

We found that TPM for HBA1 was approximately one order of magnitude higher than that of HBA2, which were 6.34 x10⁰⁵ and 8.49 x10⁰⁴ (medium TPM), respectively, indicating that HBA1 gene is mainly transcribed among a globin gene (Table 1). An intriguing observation here is that TPM for HBA1 was approximately 4-fold higher than that for HBB (Table 1 and Fig. 2A). Expression profiles of these globin genes showed almost normal distributions in 576 individuals (Fig. 2A). Consistent with these finding, similar expression profiles for these globin transcripts have been reported (33). Hematopoietic indices of these 576 participants showed no sign of anemia caused by chain imbalance of α and β globin chains (Fig. 2A). The reason(s) for the transcript level difference in peripheral blood transcriptome analysis remains to be clarified.

We also found expressions of HBG1 and HBG2 genes, which encode fetal β-type globin γ^A andγ^G, respectively. In this analysis, we assigned HBG1 and HBG2 positive if read count composition of these genes exceeds 1% of total read count composition. According to this criterion, HBG1 and HBG2 were found to be expressed in 2 (0.35%) and 7 (1.22%) of 576 participants in adults, respectively, indicating that our approach for whole blood transcriptome analysis successfully identified fetal γ globin genes in rare adults (Table 1 and Fig. 2B).

We found that two HBG1 positive cases showed high-level expressions. Especially, one case showed very high level of HBG1 and low level of HBB gene expression, while the other showed high level of HBG1 and moderate level of HBB gene expression (Fig. 2C). Expression of HBG2 in these two cases showed almost similar expression profiles, while rest of the five HBG1 positive cases showed moderate level of HBG1 TPM increase, but HBG2 level was not increased (Fig. 2D).

These globin gene expression profiles are collectively shown in Figure 2E. These results demonstrate that there are hereditary persistence of fetal hemoglobin (HPFH) cases in Japanese population and that our approach of RNA-seq including globin genes has a potential to detect these HPFH cases.

Whole blood transcriptome profile reflects cell contents in peripheral blood

In order to examine blood cell contents in whole blood sample, we next investigated gene expression profile of non-globin genes. Our hypothesis is that we can assess each white blood cell content through the expression profiles of lineage specific genes.

We examined overview of the gene expression profile based on PCA. To assess influence of the globin transcript removal, we first performed the PCA in the presence of globin gene transcript. We then recalculated the gene expression data after excluding the read information assigned to the globin gene group. The results are shown side-by-side, Figure 3A-C with globin transcript and Figure 3D–F without globin gene transcript, respectively.

When PCA was performed using data with globin transcripts, clusters highly related to age or gender group could not be observed in both transcript clusters (Fig. 3A and3D). We then examined functions of gene set, which showed strong associations to the first principal component (PC1) axis that showed high-level contribution rate of 48.8% (with globin transcript; Fig. 3B) and 28.6% (without globin transcript; Fig. 3E). Removal of globin transcript reduced the PC1 axis, perhaps due to the decrease of globin transcript contribution.

PCA with globin transcript detected functions related to ‘red blood cell function’ and ‘iron metabolism’, suggesting that the expression level information on the globin genes, which is in general excluded experimentally as a noise, has a strong influence on the entire data set (Fig. 3B). Indeed, a strong positive correlation (SCC = 0.87) was observed between HBA1 gene expression levels and the PC1 scores (Fig. 3C).

In PCA for data without the globin genes, functions detected by the PC1 were significantly different from those detected by the data with globin transcript (Fig. 3E). Comparison of PC1 data with and without globin genes showed that top 5 functions in the with-globin analysis were all disappeared in the without-globin gene analysis, which included ‘Erythrocytes take up carbon dioxide and release oxygen’, ‘Binding and uptake of ligands by scavenger receptors’ and ‘Heme signaling’ (Fig. 3B). We surmise that in the assignment of these functions, globin gene transcripts may contribute substantially.

Of note, functions related to the ‘neutrophil degranulation phenomenon’ was detected as a positively related function for PC1, while ‘Processing of capped intron-containing pre-mRNA’ ‘tRNA processing’ and ‘rRNA processing’ were detected as negatively related functions for PC1 (Fig. 3E). We selected C5AR1 (34) as a candidate transcript expressed specifically in neutrophils and examined the association between the expression level of C5AR1 and PC1 score. We found a positive correlation (SCC = 0.69) between the expression level of C5AR1 and the PC1 score (Fig. 3F), indicating that neutrophil-related genes are contributing to this association.

As the samples for this analysis were from the TMM biobank that has been collecting the data of blood cell components, we capitalized the data and conducted association studies between the PC1 score and number of various types of blood cells. As expected, we found a strong positive correlation (SCC = 0.80) between the neutrophils number and the PC1 score (Fig. 3G). In contrast, we found a strong negative correlation (SCC = -0.78) between the lymphocytes number and the PC1 score observed (Fig. 3H). This correlation was reproducible in the PC2 score of the with-globin gene transcript analysis (Table 2B), indicating the importance of globin gene transcript removal for blood cell component assessment. Notably, more strong correlation (SCC = 0.80) was also observed between the neutrophil-to-lymphocyte ratio (NLR), which is utilized as one of the indicators of immune function, and the PC1 score. These results support our contention that whole blood transcriptome results reflect contents of blood cells, which can be used to predict the proportion of the blood cell components in peripheral blood.

Substantial numbers of differentially expressed genes are found by age and gender

To understand the influence of age and gender to the gene expression profile in whole blood transcriptome, we characterized expressions of various genes in a population stratified by age and gender. To this end, we have used the transcript data without-globin-gene from 106 females in age of 20s to 30s (excluding 14 pregnant women), 172 females in age of 60–70s, 109 males in age of 20s to 30s and 175 males in 60s to 70s (Fig. 4A). In this figure, arrow A to B indicates higher level gene expression in B.

We first assessed influence of age on the whole blood transcriptome and identified 1309 genes for women (Fig. 4A, #1 plus #2) and 366 genes for men (#3 plus #4). Expression levels of these two groups of genes changed significantly with age (Fig. 4B). Of the 1309 genes identified in female, expressions of 687 genes were down-regulated, while 642 genes were up-regulated in the aged group. Of the 366 genes identified in male, expressions of 274 genes were down-regulated, while 92 genes were up-regulated in the aged group. Thus, the most salient finding here is that females have more significant gene expression changes with age than male do. In addition, males have three-times more down-regulated genes than up-regulated genes along with age, while females have almost similar numbers of down-regulated and up-regulated genes.

We then analyzed genes which expressions differed in gender and further characterized these genes in combination with age. The number of genes whose expression differed between females and males was 428 in their 20s to 30s (Fig. 4A, #5 plus #6) and 1129 in their 60s to 70s (#7 plus #8), indicating that the gender difference in gene expression in whole blood transcriptome increases with age. Of the 1129 genes, 609 showed higher level expression in males than in females, while 520 genes showed higher level expressions in females than in males (Fig. 4C). These results clearly demonstrate that number of genes with differential expression by gender increases approximately three times in the aged groups.

Functional characteristics of DEGs by age and gender

It seems of interest to examine functions of genes that resided in each category set in the previous section. We first examined genes influenced by age. As there are approximately 4 times more genes in females than in males in this category (Fig. 4B), we characterized age-influenced genes in female at first. In women, the hallmark gene sets showing age-dependent up-regulation include ‘UV RESPONSE DN’ and ‘HEME METABOLISM’ (Fig. 4D, #1), while genes down-regulated with age has no significantly enriched functions. The genes in each category are shown in Supplementary Table S1. In men, the hallmark gene sets showing age-dependent down-regulation include ‘G2M CHECKPOINT’, ‘E2F TARGETS’ and ‘MYC TARGETS V1’ (Fig. 4D, #4), while genes up-regulated with age has no significantly enriched functions.

We then conducted functional annotation of gene expressions that are changed by gender in combination with age. In the 20s to 30s, ‘HEME_METABOLISM’ was detected as functions of gene groups whose expression was higher in males than in females (Fig. 4E, #5), and ‘INTERFERON GAMMA RESPONSE’ and ‘INTERFERON ALPHA RESPONSE’ were detected as a function of genes whose expression levels were higher in females than in males (Fig. 4E, #6). It seems interesting to note that in the 60s to 70s, ‘REACTIVE OXYGEN SPECIES PATHWAY’ was detected as functions of genes whose expression levels were higher in males than in females (Fig. 4E, #7). Similarly, ‘HEME_METABOLISM’ was detected whose gene expression levels were higher in males than in females (Fig. 4E, #5), whereas ‘INTERFERON GAMMA RESPONSE’ and ‘INTERFERON ALPHA RESPONSE’ were detected as a function of genes whose expression levels were higher in females than in males (Fig. 4E, #6).

Fig. 2

RNA sequencing of globin gene family has a potential to detect HPFH patients. (A) Relationship between HBA1 and HBB gene expression levels. TPM (Transcripts Per Million) values are used in the analysis. (B) Relationship between HBG1 and HBG2 gene expression levels. (C) Relationship between HBG1 and HBB gene expression levels. (D) Relationship between HBG2 and HBB gene expression levels. (E) Gene expression composition of globin family genes.

Annotation of genes changed by age and gender

Given the significant differences in gene expression profiles we discovered by age and gender in the human whole blood transcriptome analysis, our next step was to scrutinize the expression profiles of individual genes within each category, stratified by age and gender. We aimed to compare overlaps or redundancies among DEGs. To achieve this, we refined the criteria for the DEG estimation procedure, applying FDR cutoff of 5% in this analysis. Consequently, we identified more genes than those analyzed in the previous section.

We began by comparing the differences in expressions of age-dependent genes between female and male groups. To accomplish this, we initially stratified gene expressions by age, further subdividing them by gender. This process resulted in nine groups (Fig. 5A). For instance, Group F_Y > E M_Y > E comprises genes with higher expressions in young generations (20–30s) compared to the elderly (60–70s), irrespective of gender.

Figure 3

Continues

Fig. 3

Neutrophil and lymphocyte composition mainly affects the divergence of whole blood gene expression profile. (A) PCA of gene expression profile without in silico removal of globin gene counts. (B) Pathway enrichment analysis of strongly PC1-associated genes in A (|Factor Loading| > 0.8). (C) Correlation of HBA1 gene expression level with PC1 score in A. (D) PCA of gene expression profile after in silico removal of globin gene counts. (E) Pathway enrichment analysis of strongly PC1-associated genes in D (|Factor Loading| > 0.8). (F) Correlation of C5AR1 gene expression level with PC1 score in D. C5AR1 is selected as a neutrophil marker. (G) Correlation of neutrophil content with PC1 score in D. (H) Correlation of lymphocyte content with PC1 score in D. Note the inverse relationship of neutrophil and lymphocyte content to PC1 score. (I) Correlation of neutrophil-to-lymphocyte ratio (NLR) with PC1 score in D.

Fig. 4

Gene expression alteration by age and gender in human whole blood. (A) Schematic summary of the gene expression alternation by age and gender. Note that dot and arrow represent the higher expression of the gene. (B) Number of age-associated genes stratified by gender. Numbers correspond to those in panel (A). (C) Numbers of gender-associated genes stratified by age group. (D) Gene set enrichment analysis against MSigDB hallmark gene sets for significantly age-related gene set as in (A). (E) Gene set enrichment analysis against MSigDB hallmark gene sets for of significantly gender-biased expression gene set as in (A).

Fig. 5

Interplay of age and gender on whole blood gene expression. (A) Schematic diagram of gene classification based on age-related expression level alteration stratified by gender. F and M stand for female and male, respectively. Y and E stand for young (20–30s) and elder (60–70s), respectively. (B) Schematic diagram of gene classification based on gender-related expression level alteration stratified by age group. (C) Distribution of the mean normalized expression level per sample group for gender-biased age-related DEGs. (D) Gene set enrichment analysis against MSigDB hallmark gene sets for each category in C (FDR < 0.01). (E) Distribution of the mean normalized expression level per sample group for Age-biased Gender-related DEGs. (F) Gene set enrichment analysis against MSigDB hallmark gene sets for each category in (E) (FDR < 0.01).

We then compared how expressions of genes exhibiting gender-dependent variations differ between young and elderly groups. For this analysis, we initially stratified gene expressions by gender and then further subdivided them based on age. This method also yielded nine distinct groups (Fig. 5B). For example, Group Y_F > M E_F > M includes genes where expressions are higher in females than males, regardless of age.

How age-dependent gene expressions differ in female and male groups

In this analysis, we first examined gene expression variations with age and then listed them separately in female and male groups (Fig. 5C). Regardless of the gender, expression of 99 genes decreased with age (Group F_Y > E M_Y > E), while expression of 137 genes increased with age (Group F_Y < E M_Y < E). Thus, expressions of total 236 genes were changed either positively or negatively with age irrelevant to the gender.

Fig. 6

Effect of age and gender on the expression of oxidative and inflammatory stresses related genes. (A) Expression pattern for the REACTIVE_OXYGEN_SPECIES_PATHWAY genes. Asterisk (*) indicates genes that are categorized into Group Y_F = M E_F < M in Figure 5D. Expression pattern of 4 genes in the pathway are shown; (B) NFE2L2 (NRF2), (C) GLRX (Glutaredoxin), (D) TXN (Thioredoxin) and (E) PDLIM1 (PDZ and Lim domain protein 1). Data from pregnant females in B–E also show pregnancy-induced changes of these gene expressions. P, Pregnant, NP, Non-Pregnant.

In contrast, if we restricted the gene expression changes in either female or male, much wider range of changes could be seen. In case of females, expressions of 1270 genes were found to decrease (Group F_Y > E M_Y = E), while expressions of 1282 genes were increased with age (Group F_Y < E M_Y = E). In case of males, numbers of genes in this category were markedly smaller than those in females; expressions of 433 genes (Group F_Y = E M_Y > E) decreased and expressions of 88 genes (Group F_Y = E M_Y < E) increased with age. Only one gene each for Groups F_Y > E M_Y < E and F_Y < E M_Y > E was found in this analysis, which correspond to genes whose expression increase in male but decrease in female and whose expression decrease in male but increase in female, respectively. Therefore, we did not show these two groups in the Figure 5C.

We then explored enriched functions and found nine functions in Groups F_Y > E M_Y = E, F_Y < E M_Y = E and F_Y = E M_Y > E (Fig. 5D). In the Group F_Y > E M_Y = E in which age-related decrease of gene expression occurred most drastically in females, ‘IL6 JAK STAT3 SIGNALING’, ‘PROTEIN SECRETION’, ‘INFLAMMATORY RESPONSE’, ‘COMPLEMENT’ and ‘MTORC1 SIGNALING’ were identified (Fig. 5D). In the Group F_Y < E M_Y = E in which age-related increase of gene expression occurred most drastically in females, ‘HEME METABOLISM’ was detected (Fig. 5D). This shows very good coincidence with that of #1 group in Figure 4D. In the Group F_Y = E M_Y > E in which age-related decrease of gene expression occurred most drastically in males, ‘G2M CHECKPOINT’, ‘E2F TARGETS’ and ‘MYC TARGETS V2’ were identified (Fig. 5D). Again, this shows very good coincidence with that of #4 group in Figure 4D.

How gender-dependent gene expressions differ in young and elderly groups

We next examined gene expression variations with gender in young and elderly generation groups (Fig. 5E). In this analysis, we first examined gene expression variations with gender and then listed them separately in age groups. Regardless of age, 74 genes were highly expressed in female (Group Y_F > M E_F > M), while 120 genes were highly expressed in male (Group Y_F < M E_F < M). In contrast, 216 genes were found to be expressed more in females than males in the young generation (Group Y_F > M E_F = M), and 458 genes were found to be expressed more in males than females in the young generation (Group Y_F < M E_F = M), indicating that expressions of these two groups of genes differ in gender only in the young generation.

In the 60s, 1140 genes showed higher level expression in females than in males (Group Y_F = M E_F > M) and 1142 genes showed higher level expression in males than in females (Group Y_F = M E_F < M). These results indicate that a number of genes show difference in their expression in both genders. This finding coincides well with the finding of Figure 4C.

When the enriched functions of each gene group were investigated, ‘INTERFERON GAMMA RESPONSE’, ‘INTERFERON ALPHA RESPONSE’ and ‘INFLAMMATORY RESPONSE’ were found in the Group Y_F > M E_F = M (Fig. 5F). For the Group Y_F < M E_F = M, ‘HEME METABOLISM’ and ‘EPITHELIAL MESENCHYMAL TRANSITION’ functions were significantly enriched (Fig. 5F). For the Group Y_F = M E_F > M, ‘ALLOGRAFT REJECTION’, ‘ESTROGEN RESPONSE EARLY’ and ‘MYC TARGETS V2’ were detected as enriched functions (Fig. 5F). For the Group Y_F = M E_F < M, ‘MTORC1 SIGNALING’, ‘OXIDATIVE PHOSPHORYLATION’, ‘REACTIVE OXYGEN SPECIES PATHWAY’, ‘PROTEIN SECRETION’ and ‘COMPLEMENT’ were detected as enriched functions (Fig. 5F).

While in this study we explored gene expression variability with aging and sex differences at a population level, it is also important to address gene expression variability at the individual level. As a preliminary approach, we examined gene–gene relationships within the UV RESPONSE DN and HEME METABOLISM gene sets and found that not all gene pairs from these functional groups exhibit a high correlation in expression levels. Due to the complexity of presenting scatter plots for numerous gene pairs, we calculated correlation coefficients as a summary statistic to assess gene expression similarities at the individual level (Supplementary Table S1). The expression level of one gene related to HEME METABOLISM shows a positive correlation with multiple genes in the UV RESPONSE DN set, suggesting its potential role in regulating gene expression across various biological functions. In contrast, the other gene showed a lesser association. These insights underscore the significance of individual gene expression patterns. Individual gene expression data of this study are available in the jMorp database [https://jmorp.megabank.tohoku.ac.jp].

Effect of age and gender on the expression of oxidative and inflammatory stresses related genes

In order gain more insight into the influence of age and gender on the whole blood transcriptome profiles, we next tried to focus on specific pathway of genes. For this purpose, we selected expression profiles of the REACTIVE_OXYGEN_SPECIES_PATHWAY gene set in MSigDB, which include 49 genes (Fig. 6A).

In this analysis, we first noticed that expression of NFE2L2 gene is expressed significantly higher in young females, and the expression is down-regulated with age (Fig. 6B). The NFE2L2 gene encodes NRF2, a transcription factor regulating inducible expression of antioxidant and phase II detoxification-related genes. Expression levels of the NFE2L2 gene in young and elderly males are low, which are comparable with elderly females.

To gain further insight, we compared expression levels of NFE2L2 gene in 14 pregnant females (see  Fig. 1B) with those of 106 non-pregnant females. As a result, we discovered that the NFE2L2 gene expression is significantly higher in pregnant females than in young females (Fig. 6B). As the young females show much higher level of NFE2L2 gene expression than that in elderly females and in both young and elderly males, these results unequivocally demonstrate that the NFE2L2 expression levels in whole blood transcriptome analysis increases significantly with the pregnancy.

We also found that approximately a quarter of the genes (15 out of 49 genes) in this pathway exhibit similar expression level changes to that of NFE2L2 gene, which shows higher level expression only in young females (Fig. 6A). We surmise that this group of genes may include NRF2 target genes. As expected, we found several typical NRF2 target genes including GLRX and TXN genes. Therefore, we performed similar analysis, including the pregnant females. We found that expressions of the GLRX (Fig. 6C) and TXN (Fig. 6D) genes are reproducibly and significantly higher in pregnant young females than in non-pregnant females. Again, both expressions in young females are higher than those in elderly females and both young and elderly males.

To approach to a more comprehensive understanding of our findings related to NRF2, we calculated the correlation coefficients for the gene expression levels between the NRF2 gene and its major downstream genes, including both GLRX and TXN (Supplementary Table S2). As a result, we found that detoxification and antioxidant genes show positive correlation in general, while proinflammatory genes, like TNF-α and IL-6, showed negative correlation. These results show very good coincidence to the previous analyses, as the expressions of the former genes appear to be activated by NRF2, while expressions of the latter genes are shown to be repressed by NRF2 (35).

We found several immune-related genes, such as PDLIM1, show higher level expression in elderly females than in young females, which is an opposite expression change from NFE2L2 (Fig. 6A). PDLIM1 plays an essential role in inflammatory responses, which inhibits NF-κB nuclear translocation through binding to NF-κB. Intriguingly, PDLIM1 gene expression is significantly lower in pregnant females than in non-pregnant females and in elderly females, elderly males and young males (Fig. 6E). These results suggest that anti-inflammation levels might be activated in pregnant females.

These results of the REACTIVE_OXYGEN_SPECIES_PATHWAY gene set thus revealed that expressions of oxidative and inflammatory stress-related response genes in whole blood transcriptome are under elaborate and meticulous regulations (Supplementary Fig. S3), which may have important clinical relevance.

Discussion

In this study, we aimed to prepare a foundational platform of whole blood transcriptome of the Japanese population. For this endeavor, we chose 576 individuals from the participants of the TMM Cohorts and Biobank. The selected group consists of 120 young females and 109 young males, as well as 172 elderly females and 175 elderly males. Notably, our young female group includes 14 pregnant individuals. We conducted whole blood transcriptome analysis on all these stratified groups. By combining RNA sequencing of the globin-including library and in silico removal of globin gene reads, we efficiently captured both the expression dynamics of the globin gene family and accurate expression profiles of all transcripts in whole blood samples. Our results show that a significant number of genes are differentially expressed in both the female and male groups as well as the young and elderly groups. This includes genes found in the reactive oxygen species pathway and other pathways. While recent studies have reported whole blood transcriptome analyses on Japanese subjects (24, 36, 37), to the best of our knowledge, our study is the most extensive and comprehensive whole blood transcriptome data set of the Japanese population to date.

Traditionally, whole blood transcriptome analysis involves the removal of globin gene transcripts along with rRNAs (38). However, we decided to challenge this practice in our study by analyzing the expression of globin gene families. This approach allowed us to gain fundamental insights into the expression levels of each globin gene in adults. Indeed, we found that 7 of the 576 adult subjects analyzed displayed significant levels of expression of the fetal HBG1 or HBG2 genes (encoding Hb ^Aγ and ^Gγ, respectively), indicating that the incidence of HPFH in the Japanese population is approximately 1.22%. While the mechanistic basis for HPFH is yet to be clarified, it seems interesting to note that in our long-read sequence analysis of the Japanese population, we found a roughly 5-kb deletion in between the HBG1 and HBG2 genes in 2 of 222 individuals (39). The region has been reported to contain important negative regulatory region for β-type globin gene expression (40, 41).

It is well known that gene expression profiles in human blood samples are influenced by both the age and gender of subjects (18–21, 42). Consistent with the findings, our study revealed that a substantial number of genes are differentially expressed between the female and male groups, as well as the young and elderly groups. Intriguingly, we observed that females exhibited more significant age-related changes in gene expression profile (642 genes increased and 687 decreased with age) than males do (increase and decrease with age are 92 and 274, respectively). We hypothesize that these gender difference may impact the regulation of hematopoietic cell differentiation, which becomes more prevalent with age. This contention aligns with a recent study showing that females retain a lower percentage of natural killer cells and a higher percentage of plasma cells in peripheral blood compared to males, but these gender differences diminish with age (43). It has also been reported that some of the hematopoietic transcription factors are influenced by sex hormones. For instance, GATA1 is negatively influenced by estrogen (44).

It should be noted that variations in individual whole blood transcriptome profiles seem to be influenced substantially by the numbers of leukocytes and lymphocytes present in each sample. This observation suggests that the immune system activity of individual appears to be one of the major determinants of gene expression profiles in human whole blood transcriptome analysis. Conversely, this finding also points to the potential of utilizing whole blood transcriptome results for estimating peripheral blood cell composition or deconvoluting peripheral blood cells. With this in mind, we believe that the outcomes of whole blood transcriptome analysis may have extensive applications in the future.

We identified significant enrichments of biological functions or pathways for DEGs associated with age or gender in the whole blood transcriptome analysis. For instance, heme metabolism was noted as a function of genes whose expression was lower in young females in their 20s to 30s. Consistent with this result, concentrations of heme and its related metabolites in serum were found to be higher in males than females (45). In contrast, immune response-related genes exhibited significantly higher expression levels in young females, demonstrating that the immune response activity is more pronounced in young females than in young males. This finding is consistent with previous studies comparing gene expression profiles between females and males (46). Among the genes whose expression levels decline with age solely in males, we identified cell cycle regulation-related functions, such as ‘G2M CHECKPOINT’, ‘E2F TARGETS’ and ‘MYC TARGETS V1’. This presents potentially crucial fundamental information about gender-related differences in cellular senescence and/or cancer cell growth (47). Furthermore, ‘REACTIVE OXYGEN SPECIES PATHWAY’ was detected as an enriched function of genes with higher expression levels in elderly males compared to females. This finding may link the distinct response ability of men and women to oxidative stress with the variations in their rates to suffer from cardiovascular diseases (48).

Our data set includes whole blood transcriptome data taken during pregnancy, providing us with the opportunity to examine the impact of pregnancy on gene expression in whole blood. We observed significant increases in the expression of genes, such as glutaredoxin and thioredoxin, which are associated with the oxidative stress response (49), during pregnancy. Upon detailed scrutiny of the pregnancy-induced changes in gene expression, we identified functions related to immune response and heme metabolism as the functions of genes whose expression increased during pregnancy (Supplementary Fig. S4). Recent studies have shown that the proportion of neutrophils in the blood rises during pregnancy (50). Consequently, our observations may be indicative of the higher expression of genes related to immune function as the neutrophil ratio increases.

It has been shown that both lactotransferrin (LTF) and matrix metallopeptidase 8 (MMP8) are essential for achieving various immune function (51, 52). Our study found that the expression of LTF and MMP8 is much higher in pregnant women in their 30s than in non-pregnant women in their 20s to 30s. Furthermore, we observed an increase in the expression of genes essential for heme and heme protein biosynthesis, such as hydroxymethylbilane synthase (HMBS) and alpha hemoglobin stabilizing protein (AHSP), during pregnancy. Since the genetic polymorphisms of the HMBS gene have been reported to be involved in acute intermittent porphyria (53), our data on the expression levels of these heme biosynthesis-related genes could potentially contribute to a better understanding of the mechanisms underlying iron deficiency anemia during pregnancy.

In summary, our findings offer novel insights into individual gene expression variation within the whole blood transcriptome in the Japanese population. To further promote personalized and precision medicine, it will be essential to clarify the relationship between gene expression variation and genetic variation through eQTL analysis (54).

Supplementary Data

Supplementary data are available at JB Online.

Funding

This research was supported (in part) by the Japan Agency for Medical Research and Development, AMED under grant JP21tm0424601.

Conflict of Interest

None declared.

Author contributions

Y.A. contributed to data analysis and interpretation, drafted the original manuscript. K.T. contributed to blood RNA samples, drafted the original manuscript. H.A., I.M., Ke.K. supervised data analysis. J.K., A.O., L.B., T.S., F.K. supervised interpretation in RNA-seq analysis. N.I., A.H., Ka.K. provided blood samples. K.O. supervised interpretation in hematology. M.Y. wrote the manuscript. All authors reviewed the manuscript and approved the final version of the manuscript to be published.

Acknowledgements

The authors would like to thank Takumi Fukushi, Nanae Osanai, Noriko Takahashi, Keiko Tateno and HaploPharma Inc. (Sendai, Miyagi, Japan) for technical assistance in DNBSEQ-G400 sequencer.

Author notes

^‡co-first author