The general law of plasma proteome alterations occurring in the lifetime of Chinese individuals reveals the importance of immunity

Baseline characteristic

Between 2020–2021, we recruited 253 healthy individuals, aged 0–100 years, with no adverse outcomes (such as, cardiovascular diseases, neurodegenerative diseases, and cancers) from the Beijing Hospital. We collected plasma from ethylenediaminetetraacetic acid (EDTA)-treated blood obtained via vein puncture. Samples with poor plasma quality and incomplete basic information were excluded from the study. Next, we screened samples from each age group, which consisted of a compilation of 10 years (i.e., Age group 1: 0-10; Age group 2: 11-20; and so on). Finally, we selected 97 samples. Among them, 37 were males (38.14%) and 87 (89.69%) were Han Chinese (Table 1). These subjects were well matched in sex and ethnicity (Supplementary Table 1).

Characterization of age-dependent protein expression

Mass spectrometry (MS)-based quantitative proteomics is a key technical method for protein identification and determination of the absolute or relative abundance of individual proteins [9]. In this study, we conducted MS quantitative proteomics based on Data-independent Acquisition (DIA). To quantify the proteomics data, peak strength of overall sample was calculated and normalized. Based on our analysis, the quantitative values distribution of the samples’ peptide segments achieved basically the same response strength (Supplementary Figure 2A). We examined protein quantification of all samples via LC-MS/MS to further validate our quantification data. Based on our analysis, DIA was very repetitive for large samples (Supplementary Figure 2B), suggesting that our proteomics technique was accurate and precise. In total, we identified 2682 age-related proteins (89.28%), out of the 3004 total proteins in plasma (false discovery rate [FDR]-corrected p≤0.01) (Supplementary Figure 2C). Following appropriate sample quality control (QC) and normalization procedures, we employed the Pearson correlation coefficients (PCC) to assess correlations among the 97 samples, based on proteomics data (2682 proteins). Apart from the A20-2 sample belonging to the Zhuang ethnic group, all PCCs between the technical replicates were relatively high, reaching mean of 0.97 (Supplementary Figure 2D).

Construction of weighted gene co-expression network (WGCNA)

WGCNA was conducted identify co-expression modules and generate co-expression networks. To conduct proteomics analysis, 96 samples were divided into 10 groups (B1-B10), with each group representing a compilation of 10 years old, following the removal of sample A20-2. Next, based on the expression level of all proteins, 10 sets of samples were clustered hierarchically. We observed no strong outliers (Figure 1A). Next, we employed adequate soft-threshold validation, with convergence toward scale-free topology and mean connectivity, to select 6 common values as power values (Figure 1B). The Dynamic Tree Cut algorithm provided the first set of modules and this was combined with associated modules (r > 0.7; Merged dynamic). Each module contained at least 50 genes (Figure 1C and Supplementary Figures 3, 4). 14 modules were identified according to age (Figure 1D).

Figure 1. Network construction of the weighted co-expressed genes and their associations with age. (A) The Group hierarchical clustering tree. (B) Determination of soft-thresholding power. (C) The heatmap of genes correlation in modules. (D) Hierarchical clustering dendrograms (HCD) of modules. HCD using dynamic tree cut, depicting 14 plasma protein modules. Each branch of the HCD presents a single protein, and the colored bar below refers to its corresponding protein module. The HCD height is the distance between proteins.

Correlation between modules and aging

The gene co-expression in each module was summarized by Eigengene (i.e., the first component gene expression of each module). Next, we computed the relationship between each Eigengene and age group. As presented in Figure 2A and Supplementary Figure 5, the strongly correlated modules were different for each age group; for instance, the brown module represented 0-10 years old, the black module represented 11-20 years old, the purple module represented 21-30 years old, the pink module represented 31-40 years old, the midnight blue module represented 41-50 years old, the green module represented 51-60 years old, the lightcyan module represented 61-70 years old, the blue module represented 71-80 years old, the red module represented 81-90 years old, and the turquoise module represented 91-100 years old. In the meantime, the intramodular connectivity (K.in) and module correlation degree (MM) corresponding to each gene were computed and genes with high connectivity were assigned as potential hub genes with important functions. These results are presented in Supplementary Figures 6, 7.

Figure 2. Heatmaps of age-related modules and the clustering of some critical module eigengene. (A) Heatmaps of the correlation between Eigengene and clinical traits of age. Each row corresponds to a module eigengene, and each column corresponds to an age group. Each cell contains the corresponding correlation and p-value. (B) The clustering of four module eigengenes (black, green, red and turquoise).

Functional enrichment analysis of critical modules

Next, we performed Gene Ontology (GO)-biological process (BP) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses to elucidate the potential functional implications of proteins identified in the strongly correlated modules of each age group.

The top five BP-related GO terms from each age group demonstrated that the proteins were involved in complement activation, protein activation cascade, and immune response (GO:0006956, GO:0072376, GO:0002455) in the 0-10 years old (module brown); leukocyte mediated immunity (GO:0002443), immune response (GO:0006955), and cell activation (GO:0001775) in the 11-20 years old (module black); nucleosome or chromatin assembly (GO:0006334, GO:0031497) in the 21-30 years old (module purple); keratinization and epidermis development (GO:0031424, GO:0008544) in the 31-40 years old (module pink); regulation of receptor binding (GO:1900121, GO:1900120), heart trabecula morphogenesis (GO:0061384), negative regulation of molecular function or high-density lipoprotein particle clearance (GO:0044092, GO:0010987) in the 41-50 years old (module midnight blue); exocytosis (GO:0045055, GO:0006887) and vesicle-mediated transport (GO:0016192) in the 51-60 years old (module green); keratinization (GO:0031424, GO:0030216) in the 61-70 years old (module lightcyan); platelet degranulation or activation (GO:0002576, GO:0030168) in the 71-80 years old (module blue); protein activation cascade, complement activation, and immune response (GO:0072376, GO:0006956, GO:0002455, GO:0002250, GO:0006958) in the 81-90 years old (module red); chromatin silencing, epigenetic, and immune system process (GO:0006342, GO:0045814, GO:0000183, GO:0002376) in the 91-100 years old (module turquoise) (Figure 3A and Supplementary Table 3).

Figure 3. The top enrichment of Gene Ontology (GO)-biological progress (BP) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway in each age groups. (A) The leading five enrichment data of GO-BP analysis. Differential genes (Genes) are distributed across multiple GO biological progress. Each row corresponds to the top five enriched pathways, and each column corresponds to an age group. The color represents P value and the smaller the P value, the color becomes darker (red). no enrichment to no color (White). There is no color (grey) without enrichment. (B) Top ten enrichment results of the KEGG pathway. Genes are involved in a wide range of KEGG. Each row corresponds to top ten enriched pathways, and each column corresponds to an age group. The color represents P value and the smaller the P value, the more obvious enrichment, the darker color (red), and no enrichment to no color (White). There is no color (grey) without enrichment. (C) The pathway crosstalk analysis of critical modules. The color of nodes denotes modules.

We also identified enriched pathways using KEGG enrichment analysis. Based on our analysis, infectious diseases (ko05150) were the most significantly enriched KEGG pathway, followed by transport and catabolism (ko04145) in the 0-10 years old (module brown). With the development of the body, signaling molecules and interaction (ko04512) with environmental information began to play key roles among the 11-20 years old (module black). The purple module, corresponding to the 21-30 years old, showed significant enrichment of pathways involved in human diseases, such as, immune diseases, cancers and substance dependence (ko05322, ko05203, ko05034). Subsequently, proteins were enriched in development and immune system (ko04380, ko04662) in the 31-40 years old (module pink), followed by infectious diseases (ko05146, ko05150) and digestive system (ko04976) in the 41-50 years old (module midnight blue). In the middle-aged and old stages [10] (51-60 years old, module green and 61-70 years old, module lightcyan), the proteins were significantly enriched in carbohydrate metabolism (ko00030, ko00010), transport, catabolism (ko04145, ko04144), and immune system (ko04666). These alterations implied that the cellular processes and metabolism were mobilized throughout the body. As the aging process continued, the proteins significantly enriched in cell motility and cellular community of cellular processes (ko04810, ko04510, ko04520), infectious diseases and cardiovascular diseases of human diseases (ko05131, ko05414) in the 71-80 years old (module blue). Next, proteins were enriched in immune and infectious diseases (ko05150, ko05143, ko05322, ko05146, ko05310) in the 81-90 years old (module red). Finally, in individuals over 90 years old, proteins were enriched in immune diseases, substance dependence, and cancers (ko05322, ko05034, ko05202; module turquoise; Figure 3B and Supplementary Table 4). Notably, the >90 age group was enriched in the same pathways as the 21-30 age group, however, the genes involved were different.

The pathway crosstalk analysis of critical modules

Based on the overall analysis over an average lifespan, we performed pathway crosstalk analysis among the critical module in each age group. Notably, these modules were not independent, as they shared certain common pathways, for example, phagosome (ko04145), pentose phosphate pathway (ko00030), staphylococcus aureus infection (ko05150). All pathways in these modules formed two major clusters (Figure 3C). One cluster consisted of the pentose phosphate pathway and glycolysis/gluconeogenesis and involved the module green and lightcyan. The other cluster was more complex and formed connections with all critical modules. Interestingly, most of the signaling pathways in these clusters were related to immunity. Viral carcinogenesis, which was ranked as the topmost node (20 nodes) and most interactive, can serve critical roles in the progression of life. These results demonstrate that immunity is important throughout life.

Identifying hub genes in the critical module

Furthermore, we generated a protein-protein interactions (PPI) network of proteins in the critical module of each age group (Table 2 and Supplementary Figure 8) to recognize chief hub genes associated with aging.

The all-connective genes in the ten modules were individually uploaded to STRING for PPI evaluation (Supplementary Figure 8). In all, the node scores ranged from 0.400 to 0.999. Using the PPI data, we were also able to establish the co-expression coefficient of different proteins. Based on the high connectivity of genes in each module, the elevated combined scores (CS; value > 0.9) and elevated co-expression coefficients (CEC; value > 0.8) suggested that protein-protein interactions dictate spatial protein positioning and gene expression within cells. We further analyzed the genes in relation to our previously identified hub genes in each age group. The molecular functions of all candidate hub genes are listed in Table 2. And, elevated CS (value > 0.9), but reduced CEC (value = 0) suggests that the protein-protein interaction not only regulates gene expression, but also signaling cascade (Supplementary Table 5).

Based on our results, there were forty hub genes in the module brown, black, green, blue, red, and turquoise with elevated CS combined scores (value > 0.9) and elevated CEC (value > 0.8) (Table 3). The proteins in these modules represented strong correlation with the corresponding age groups (Figure 4A). Interestingly, hub genes from age groups of 0-20 years old and 71-100 years old, namely, ARPC2, UBB, FGB, and FGG, are mostly involved in infectious diseases and the immune system. In contrast, the hub genes from the age groups of 51-60 years old, namely, ADH1A, PGD, and TALDO1, are mainly associated with carbohydrate metabolism (Table 3 and Supplementary Table 6).

Figure 4. The expression of hub genes and transcription factor (TF) regulatory network in the critical module. (A) The expression heatmap of hub gene. (B) The interaction regulatory network of TF and target genes. The core circles are two transcription factors (CDC5L and HMGB2). The surrounding circles represent their target genes, the color represents the source module, the circle size represents the expression, and the yellow circle represents the hub gene.

Construction of the transcription factor (TF) regulatory network

We further investigated the molecular mechanisms that regulate hub genes. We employed the Network Analyst tool to predict TFs that regulate expression of the forty hub genes. Based on our analysis, the TFs CDC5L and HMGB2 modulated the expressions of ADH1A, ARPC4, PGD, HSPA5, APOA1, ITIH2, APOA2, TPT1, FGB, F2, FGG, ACTC1, and ACTBL2 genes (Figure 4B) in the 51-60 years old and 71-100 age groups (modules blue, green, red, and turquoise). In addition, PCC analysis revealed that CDC5L was positively correlated with ADH1A, ACTC1, APOA1, ITIH2, and FGG expressions, whereas, HMGB2 was positively correlated with ARPC4, PGD, HSPA5, APOA1, ITIH2, APOA2, TPT1, FGB, F2, and ACTBL2 expressions.

Subjects

Samples in this study were selected from 253 healthy individuals, aged 0–100 years, without infectious diseases and adverse outcomes like cardiovascular diseases, neurodegenerative diseases, and cancers. We collected plasma from EDTA-treated blood obtained via vein puncture. According to the inclusion / exclusion criteria (inclusion criteria: 0-100 years old, healthy individuals with complete basic information; exclusion criteria: 0-100 years old, with infection or major chronic diseases, or with incomplete basic information), 153 individuals were screened out. Finally, we selected 97 samples (37 males and 60 females, 87 Han Chinese and 10 other Ethnic Groups) according to the principle of 10-year-old group (i.e., Age group 1: 0-10; Age group 2: 11-20; and so on), equal distribution of sample number and gender (Table 1). Subjects in the ten groups were well matched in sex and ethnic by chi-square test (Supplementary Table 1). We received informed consent from all subjects and their guardians, retrieved detailed information about this research and any risks associated with it before subject recruitment. Our protocols followed the Declaration of Helsinki and received approval from the ethical board at Beijing Hospital, Ministry of Health (2019BJYYEC-118-02).

Sample pretreatment for plasma proteome profiling using DIA MS

The plasma sample (20 μL for serum sample) was diluted with 1% SDS in 8 M urea buffer 40 times and allowed to lyse on ice for 30 min, while vortexing the tube every 10 min. The soluble protein lysates were then centrifuged at 12000 g for 20 min at 4° C and the accumulated protein quantified with a BCA protein assay (Supplementary Figure 1 and Supplementary Table 2). Next, 100 μg of protein was placed in a fresh Eppendorf tube with 8 M urea at a total volume to 100 μL. 2 μL of 0.5 M TCEP was then introduced and the sample maintained at 37° C for 1 h, before addition of 4 μL of 1 M iodoacetamide and maintenance in the dark at room temperature (RT) for 49 min. Following this, 5 volumes of -20° C pre-chilled acetone was employed for overnight (ON) protein precipitation at -20° C. The precipitated protein was then rinsed with pre-chilled 90% acetone aqueous solution two times before re-suspension in 100 μL 100 mM TEAB. This was followed by the addition of Sequence grade modified trypsin (Promega, Madison, WI; ratio 1:50, enzyme: protein, weight: weight) for ON protein digestion at 37° C. The peptide mixture was then desalted using C18 ZipTip, measured with Pierce™ Quantitative Colorimetric Peptide Assay (23275), and lyophilized using SpeedVac.

High PH reverse phase separation

The peptide mixture was re-suspended in 20mM ammonium formate in water, pH brought up to 10.0 with ammonium hydroxide. Next, it underwent fractionation via high pH separation using Ultimate 3000 system (Thermo Fisher Scientific, MA, USA) attached to a reverse phase column (XBridge C18 column, 4.6mm x 250 mm, 5μm, (Waters Corporation, MA, USA). It employed a linear gradient from 5-45% in 40 min of 20mM ammonium formate in 80% ACN, pH brought up to 10.0 with ammonium hydroxide. Column re-equilibration was done for 15 min to bring it back to the starting condition. The column flow rate (CFR) and temperature (CT) remained at 1mL/min and 30° C, respectively. We collected a total of ten fractions, each of which was vacuum-concentrated for subsequent analysis (Guangzhou Gene Denovo Bio-Tech Co., Ltd., China).

DDA: Nano-HPLC-MS/MS analysis

The peptides were re-suspended in 30 μL 0.1% formic acid in water and assessed via on-line nanospray LC-MS/MS on an Orbitrap Fusion Lumos coupled to EASY-nLC 1200 system (Thermo Fisher Scientific, MA, USA). 3 μL peptide was introduced to the analytical column (Acclaim PepMap C18, 75 μm x 25 cm) and peptide separation was done with a 120-min gradient spanning 5-35% concentration of 0.1% formic acid in CAN. The CFR was adjusted to 200 nL/min and the CT to 40° C. The electrospray voltage of 2 kV versus the MS inlet was employed.

MS was carried out in data dependent acquisition mode, which automatically changed between MS and MS/MS modes. The adjustment criteria were as follows: (1) MS: scan range (m/z) =350–1200; resolution=120,000; AGC target=400000; maximum injection time=50 ms; Filter Dynamic Exclusion: exclusion duration=30s; (2) HCD-MS/MS: resolution=15,000; AGC target=50000; maximum injection time=35 ms; collision energy=32.

Database search

Raw DDA data processing and analysis was done with Spectronaut X (Biognosys AG, Switzerland), using default settings, to develop the initial search list. The search criteria included homo sapiens and trypsin as the digestion enzyme. Carbamidomethyl (C) was assigned as the fixed alteration. Oxidation (M) was assigned as the variable alteration. The Q value (FDR) threshold for precursor and protein levels was set to 1%.

DIA: Nano-HPLC-MS/MS analysis

The peptides were re-suspended in 30 μL of 0.1% formic acid in water and assessed as described above. The parameter adjustments were similar to the ones described above, except for a few alterations as is listed below: (1) MS: AGC target=1e6; (2) HCD-MS/MS: resolution=30,000; AGC target=1e6; collision energy=32; stepped CE=5%. (3) DIA was conducted with variable Isolation window, and each window overlapped 1 m/z, and the window quantity was 60.

Data analysis

Data preparation and analysis were done as described above. We chose dynamic iRT as the retention time estimation type. Data retrieval via Spectronaut Pulsar [28] employed a large-scale mass calibration that used iRT calibration and gradient stability to generate an ideal dynamic extraction window. Precursor and protein q value (FDR) threshold was set to 1%. We also selected “mutated” decoy generation, which applied an arbitrary quantity of random numbers of AA position swamps (min=2, max=length/2). Subsequently, we quantified all identified precursors that passed the adjusted filters. The mean top 3 filtered peptides which passed the 1% Q value threshold were entered into the calculation of the key group quantities. Missing values are not filtered or filled in our study. The specific results of quantitative peptides and expression differential analysis, please refer to Supplementary Figure 2.

Relationship analysis of samples

Correlation analysis (CA) of replicas

CA was conducted via R. This analysis produces a correlation coefficient (CC), which enables users to evaluate the reliability and operational stability of two parallel experiments. Experiment repeatability enhances the closer CC is to 1.

WGCNA generation and age-related modules identification

Network construction

Co-expression networks were generated with the WGCNA (v1.47) package in R [29]. Following the elimination of irrelevant genes, the DEG values were entered into WGCNA to generate co-expression modules, based on the blockwise automatic network generation tool. Modules were thus generated with mostly default settings, except for a few adjustments as follows: power: 6; TOM Type: unsigned; merge Cut Height: 0.7; and min Module Size: 50. Subsequently, genes placed in clusters in 14 associated modules.

Module and gene selection

To identify clinically relevant modules, we employed module eigengenes to compute the CC of samples or sample traits. Intramodular connectivity (K.in) and module correlation degree (MM) of each gene were computed with the R package of WGCNA and genes with relatively high connectivity were identified as potential hub genes carrying clinically relevant functions. PCC between each gene and corresponding trait information was computed for all modules and gene significance value (GS) was measured for all genes. The generated networks were visualized in Cytoscape_3.3.0 [30].

Functional analysis of module genes

We also performed GO and KEGG pathway enrichment analysis on genes from all modules to discern the biological activity of the relevant genes in each module.

Functional analysis of module genes

Gene Ontology (GO) [31] is an internationally recognized gene function stratification system. It provides the most recent definitions and properties of any genes and their products in an organism. GO examines three major categories: molecular function (MF), cellular component (CeCom), and biological process (BP). Its smallest unit is GO-term, which refers to a specific type of ontology. GO enrichment analysis enables us to determine which GO terms are markedly enriched in a given module, relative to the genome background. Using this information, one can screen for module genes that represent certain BP. We first mapped all module genes to the GO terms in the Gene Ontology database (http://www.geneontology.org/). Next, we calculated gene numbers for every term, before screening for markedly enriched GO terms, relative to background genome, as evidenced by the hypergeometric test. P-value was determined as follows

$$P=1-{\displaystyle \sum _{i=0}^{m-1}\frac{\left(\begin{array}{c}M\\ i\end{array}\right)\left(\begin{array}{c}N-M\\ n-i\end{array}\right)}{\left(\begin{array}{c}N\\ n\end{array}\right)}}$$

Here, N represents GO annotation-associated gene quantity; n denotes gene quantity in the module N; M represents the total quantity of specific GO term-related genes; m denotes gene quantity in the module M. The computed p-value underwent FDR Correction, considering FDR ≤ 0.05 as the cut-off point. GO terms meeting these criteria were assigned as markedly enriched GO terms in module genes.

Pathway enrichment analysis (PEA)

Protein-protein interactions are crucial regulators of biological functions. PEA is an excellent tool for elucidating the biological roles of genes in cells. KEGG is a commonly used public pathway-associated database [32]. Using PEA, we recognized essential signaling pathways associated with the Genes in the ten examined modules, relative to a genome background. The formula used to perform PEA was the same as GO terms, only instead of entering genes associated with GO terms, we entered genes associated with the KEGG annotation.

Protein-protein interaction (PPI) network generation using select modules and identified hub genes

PPI network was generated via STRING (version 10.0; available online: http://string-db.org) [33], which plotted genes as nodes and interaction as lines, and it was visualized via the Cytoscape software [30], which displayed a core and hub gene biological interaction. The generated PPI network represented a database of known and predicted protein interactions, and only interactions with a combined score >0.9 were marked as statistically significant.

Protein advanced functional annotation

The protein domain analysis (Pfam annotation)

Pfam (Protein families database of alignments and hidden Markov models) is a database that includes protein structural information, based on the widely used multiple sequence alignment and hidden markov model (HMM). Protein domain prediction employed the Pfam_scan [34] program developed by sanger. The protein sequence was compared with the Pfam database to obtain relevant annotation information of protein structures.

TF analysis

TFs regulate gene expression and are critical for development, cellular function, and disease response. The predicted protein sequences were compared with the TF database (plantTFDB [35]/animalTFDB [36]) using hmmscan.

Subcellular localization analysis

Protein activities are highly reliant upon their subcellular localizations. Hence, predicting protein localization directly from protein sequences can be highly beneficial in inferring protein functions. The software WoLFPSort [37] was used to predict the subcellular location of proteins and study their functions.

Availability of data and material

The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.