# Noninvasive Analysis Using Data-Independent Acquisition Mass Spectrometry: New Epidermal Proteins That Reveal Sex Differences in the Aging Process

Authors: Shirui Chen; Hui Zhang; Mengting Liu; Yaochi Wang; Cong Xin; Jing Ma; Xiaodong Zheng; Xuejun Zhang; Liangdan Sun; Sen Yang
Journal: Oxidative Medicine and Cellular Longevity
 (2021)

---

## Abstract

The development of mass spectrometry has provided a method with extremely high sensitivity and selectivity that can be used to identify protein biomarkers. Epidermal proteins, lipids, and cornified envelopes are involved in the formation of the skin epidermal barrier. The epidermal protein composition changes with age. Therefore, quantitative proteomic changes may be indicative of skin aging. We sought to utilize data-independent acquisition mass spectrometry for noninvasive analysis of epidermal proteins in healthy Chinese individuals of different age groups and sexes. In our study, we completed high-throughput protein detection, analyzed protein differences with MaxQuant software, and performed statistical analyses of the proteome. We obtained interesting findings regarding ceruloplasmin (CP), which exhibited significant differences and is involved in ferroptosis, a signaling pathway significantly associated with aging. There were also several proteins that differed between sexes in the younger group, but the sex differences disappeared with aging. These proteins, which were associated with both aging processes and sex differences, are involved in signaling pathways such as apoptosis, oxidative stress, and genomic stability and can serve as candidate biomarkers for sex differences during aging. Our approach for noninvasive detection of epidermal proteins and its application to accurately quantify protein expression can provide ideas for future epidermal proteomics studies.


---

## Body

## 1. Introduction

Since 1980, many studies have demonstrated that the epidermis is structurally and biochemically diverse in terms of its metabolic state while playing key roles in epidermal barrier function and skin senescence. Epidermal lipid and protein compositions and the number of stratum corneum (SC) layers affect these processes [1–3]. Skin aging is a complex physiological process in which several mechanisms simultaneously participate. Both inherent aging and photoaging affect changes in epidermal proteins such as collagen and other extracellular matrix proteins.In past studies, genomics and transcriptomics have been used to investigate several biological pathways of aging and the expression of human aging-related genes [4–7]. Proteomics research on skin aging is still scarce. Epidermal proteins are involved in the skin barrier [8], and proteomics studies provide more direct information about biological pathways than other types of studies. Epidermal proteins continue to change with age [9]. Therefore, research on epidermal proteins can provide insight into biological markers of aging [10]. Many proteins, such as members of the sirtuin family [11], are involved in aging-related signaling pathways. Among these pathways, the most well known are the mTOR, AMPK, WNT, MAPK, and p53 signaling pathways, which ultimately affect the cell cycle, genomic instability, apoptosis, and ultimately the outcome of aging [12, 13]. There are also differences in the rates of aging between the sexes that may be related to oxidative stress, hormone levels, immunity, and underlying mechanisms related to the sex chromosome [14, 15]. Men have thicker skin than women, but women have thicker subcutaneous tissue than men, and men’s skin is more susceptible to environmental stress and UV radiation than women’s skin [16, 17]. Women also generally live longer than men [18], and studies have shown that proteins that change with age also differ between sexes. Notably, aging does not proceed at a linear speed; there are obvious inflection points at different ages, and there are differences between men and women [19]. Many studies thus far have suggested that the speed of aging differs between men and women, so research on sex- and age-related intersections will be very valuable. Therefore, we evaluated the proteomes of different sexes at different ages and observed whether there were correlations between them.In this study, we analyzed the differences in the epidermal protein expression profiles of 20 healthy men and women in China with a data-independent acquisition (DIA) method [20, 21]. We selected the age categories with consideration of the significant changes that occur in the skin during menopause [22]. Noninvasive techniques were used to obtain epidermal proteins and to explore proteins with differences and associations at the levels of both age and sex. To better identify aging-sex correlates, we performed statistical analyses and linked information about these changes to aging-related signaling pathways. The highlighted proteins may be sex-specific markers of aging, and their identification may help elucidate the reasons for the difference in the rate of aging between men and women.


## 2. Materials and Methods

### 2.1. Study Participants

Skin samples were collected from 9 healthy Chinese males and 11 healthy Chinese females, none of whom had experienced excessive light exposure. The participants were divided into two groups according to age: the 20–30 y group (3 males and 7 females) and the 55–70 y group (6 males and 4 females). Individuals were excluded from the study if they had an allergy to tape, were long-term outdoor workers, or had a skin disease or other systemic condition involving skin problems. The participants did not use topical emollients or other cosmetic products for 24 hours prior to the experiment. This study was conducted in accordance with the recommendations of the Medical Ethics Committee of Anhui Medical University, and informed consent was obtained from all enrolled subjects.


#### 2.1.1. Materials

Sodium dodecyl sulfate (SDS), L3 lysate without SDS, trypsin, dithiothreitol (DTT), Coomassie Brilliant blue G-250, 40% ethanol, 10% acetic acid, ethylenediaminetetraacetic acid (EDTA), and acrylonitrile (ACN) were obtained from Wallis Technology, Beijing, China. Ultrafiltration membranes (10K MWCO, 1.5 ml, plate) were obtained from Pall Corp. (NY, USA). Pierce C18 pipette tips (10μl bed), Empore™ C18 47 mm extraction discs (model 2215), a ThermoMixer (MS-100), and a CentriVap vacuum concentrator and accessories were obtained from Thermo Fisher Scientific (Shanghai, China). The samples were analyzed using a Q Exactive high-frequency mass spectrometer and an UltiMate 3000 high-performance liquid chromatography system (Thermo Fisher Scientific, San Jose, CA, USA).


#### 2.1.2. Sample Preparation

The palm side of the forearm was gently wiped with a sterile cotton ball to cleanse it of epidermal contaminants. A3.0∗3.0cm piece of a 3M medical tape was prepared to fully and precisely cover the back of the slide. The slide was then moved back and forth quickly with even pressure for 2 minutes. Even pressure was applied so that the slide was covered. Once the sample was collected, another slide was placed over it to protect the skin sample. To reduce sampling variation, the same technician collected samples from all volunteers during the study.


#### 2.1.3. Protein Extraction and Enzymolysis

(1) Protein Extraction. Small pieces (0.5cm∗0.5cm) of tape/skin samples were cut out of the slides with a sterile blade and then placed into corresponding 1.5 ml centrifuge tubes according to the sample numbers. An appropriate amount of L3 lysis buffer without SDS was added, and cocktail containing EDTA was added to a final concentration of ×1. The samples were placed on ice for 5 min, and DTT was added to a final concentration of 10 mM. The samples were soaked overnight and then centrifuged at 25,000×g and 4°C for 15 min. The supernatant was obtained, and DTT was added to a final concentration of 10 mM. The samples were incubated in a 56°C water bath for 1 hour, after which iodoacetamide (IAM) was added to a final concentration of 55 mM. The samples were placed in a dark room for 45 min and then centrifuged at 25,000×g at 4°C for 15 min to obtain the supernatant, which contained the proteins.(2) Protein Extraction and Quality Control. The protein concentrations were measured using the Bradford method [23] with acceptable enzymatic efficiency, and it was determined that the extracted protein content was sufficient. For the Bradford quantitative assay, protein standards (0.2 μg/μl BSA; 0, 2, 4, 6, 8, 10, 12, 14, 16, and 18 μl) were added to wells A1 through A10 of 96-well MICROLON ELISA plates, and then, pure water (20, 18, 16, 14, 12, 10, 8, 6, 4, and 2 μl) was added. Then, 180 μl of Coomassie Brilliant blue G-250 quantitative work solution was added to each well. A linear standard curve was prepared based on the optical density at 595 nm (OD595) and the protein concentration. The protein sample solutions were diluted several times, 180 μl of the quantitative solution was added to 20 μl of protein solution, and the OD595 was read. The protein concentration of each sample was calculated according to the standard curve and the sample OD595. The purity of the extracted proteins was verified by SDS-PAGE and Coomassie Brilliant blue staining.(3) SDS-PAGE. For each sample, 30 μg of protein was added to an appropriate amount of loading buffer. After thorough mixing, the proteins in buffer were heated at 95°C for 5 min, centrifuged at 25,000×g for 5 min, and placed into the wells of a 12% SDS polyacrylamide gel. After electrophoresis, Coomassie Brilliant blue staining was conducted for 2 hours. An appropriate amount of decolorization solution (40% methanol and 10% acetic acid) was added to each sample, and the samples with decolorization solution were placed in a shaker 3 to 5 times for 30 min each.(4) Proteolysis. Trypsin (2.5 μg) was added to 100 μg of protein for each sample at a protein : enzyme ratio of 40 : 1, and the protein was hydrolyzed at 37°C for 4 hours. Trypsin was added once more according to the above ratio, and enzymolysis was continued for 8 hours at 37°C. The enzymatically digested peptides were desalted with a Strata-X column and extracted under vacuum. With this method, we extracted sufficient amounts of proteins with qualified enzymatic efficiency.


#### 2.1.4. High-pH Reverse-Phase Separation

A pooled sample was created from 10μg of each sample, and 200 μg of the pooled sample was mixed with 2 ml of mobile phase A (5% acetonitrile, ACN, pH 9.8) and diluted into a Shimadzu LC-20AB liquid chromatography system. The samples were separated in a liquid phase on a 5 μm 4.6×250mm Gemini C18 column. Elution was conducted at a flow rate of 1 ml/min with the following gradient: 5% mobile phase B (95% ACN, pH 9.8) for 10 min, 5% to 35% mobile phase B for 40 min, 35% to 95% mobile phase B for 1 min, mobile phase B which lasted for 3 min, and equilibration at 5% mobile phase B for 10 min. The elution peaks were monitored at 214 nm, and one fraction was collected every minute. The samples were combined with the chromatogram elution peaks to obtain 10 fractions, which were then extracted by freezing.


#### 2.1.5. High-Performance Liquid Chromatography

The extracted peptide samples were redissolved (centrifuged with mobile phase A (2% ACN and 0.1% formic acid (FA)) at20,000×g for 10 min), and the supernatant liquid was sampled. Separation was performed with a Thermo UltiMate 3000 UHPLC. Each sample was first injected into a trap column for enrichment and desalting before being run through a self-loading C18 column (150 μm I.D., 1.8 μm pore size, 25 cm column length) coupled to tandem mass spectrometers. The samples were separated at a flow rate of 500 nl/min through the following effective gradient: 0–5 min, 5% mobile phase B (98% ACN, 0.1% FA); 5–160 min, linear increase in mobile phase B from 5% to 35%; 160–170 min, 35% to 80% mobile phase B; 170–175 min, 80% mobile phase B; and 176–180 min, 5% mobile phase B. The output of the liquid chromatograph was connected directly to the mass spectrometer.


#### 2.1.6. Data-Dependent Acquisition (DDA) and Data-Independent Acquisition (DIA) Mass Spectrometry Analyses

After liquid chromatography, the peptides were separated by nanoelectrospray ionization (ESI) and detected in a DDA mode with a Q Exactive HF (Thermo Fisher Scientific, San Jose, CA) tandem mass spectrometer. The main parameter settings were as follows: the ion source voltage was 1.6 kV, the primary mass spectrometer range was 350~1,500 m/z, the resolution was 60,000, the MS2 starting m/z was fixed at 100, and the resolution of the second-stage mass spectrometer was 15,000. The screening criteria for the precursor ions in the MS2 scans were a charge of 2+ to 7+ with a peak intensity of over 10,000 and a ranking in the top 20. The main parameter settings for DIA mode detection were as follows: the ion source voltage was 1.6 kV, the MS1 scanning range was 350–1,500 m/z with a resolution of 120,000, and the 350–1,500 Da range was divided into 40 windows for fragmentation and signal acquisition. Higher-energy collisional dissociation (HCD) was used for the DDA and DIA modes of mass spectrometric detection of ion fragmentation patterns. The fragmented ions were detected in an Orbitrap. The dynamic exclusion time was set to 30 s. The automatic gain control (AGC) settings were as follows: level 1: 3E6 and level 2: 1E5.


### 2.2. Analysis of Epidermal Proteomics Data

The DDA data from the machine were identified using the Andromeda search engine integrated with MaxQuant [24], and Spectronaut was then used to build a spectral library with the results. The resulting data were reviewed using the UniProtKB/Swiss-Prot Homo sapiens proteome database. For large-scale DIA data, after constructing the spectral library information, convolutional extraction of the data was completed using Spectronaut, and data analysis and quality control were completed using the mProphet algorithm, resulting in reliable protein quantification results. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway functional annotations were performed. Based on the high-quality quantitative results, we searched for differentially expressed proteins (DEPs) between groups.MaxQuant was used to complete the identification of the DDA data in order to create the spectral library for subsequent DIA data analysis. During this procedure, the original offline data were used as the input data, the corresponding parameters and database were configured, and then, the identification and quantitative analysis were performed. The identification information that met the criterion of a false discovery rateFDR≤1% was used to build the final spectral library. The offline DIA data were analyzed using Spectronaut [25], and iRT peptides were used to correct the retention time. Spectronaut integrates the mProphet scoring algorithm, which can accurately reflect the degree of matching of isolated pairs. Then, based on the target-decoy model applied by SWATH-MS, the false-positive rate was controlled at 1% FDR, which yielded significant quantitative results.This procedure was used to preprocess the data according to the set comparison group, and a significance test was then performed based on the model. After that, the DEPs were screened with the criteria of a fold change (FC)value≥1.5 and a p value < 0.05 (the top 10 DEPs are listed in Table 1). The R package MSstats [26] from the Bioconductor resource library was used to complete differential analysis of the DEPs.Table 1
Top 10 differentially expressed proteins.
UpregulatedDownregulatedProtein IDlog2FCp valueProtein IDlog2FCp valueOld vs youngLG3BP0.7000.001ABRX2−0.7500.001HBA3.6060.001PLD3−0.6980.003RENR1.0140.001VTI1B−0.5870.004APOA23.0690.001VAT1−0.7830.004CLUS1.6670.001KRT35−1.6450.004SG1D22.6610.002BAF−2.4920.005SEM7A1.4140.002ECM1−0.6310.009B2MG1.5900.002K1C15−0.6000.010S10A40.7340.002MYH9−0.9070.011HBB3.4620.003PSB7−0.6780.012


### 2.3. Screening for DEPs

After the raw data were median-normalized by the software, the resulting data were selected for follow-up analysis. DEPs were identified through symmetrical scatter plot screening. The dots represent the relative mean expression values (normalized andlog2 transformed) of the proteins in the two groups, and the log2FC values were calculated using the R package MSstats. The dashed line in the graph represents the threshold line at ∣log2FC∣=1.5. Thus, the volcano plot of DEPs was produced with the average protein expression values, log2FC values, and p values. The color of the dot for each protein indicates the significance of the p value, with the blue to red gradient indicating p values ranging from nonsignificant to significant. The graph shows information about the gradient change. We performed protein-protein interaction network analysis on all the DEPs and screened the proteins with more than 10 nodes of neighboring proteins in the protein relationship network. The proteins with the most significant differences in log2FC values and p values were also selected (Table 2). A gene-level network analysis was conducted for the selected most promising candidate proteins, and the biological pathways in which they were enriched were investigated by using NetworkAnalyst 3.0 inline tool [27, 28], the value of the interaction relationship between genes retained the interaction relationship with experimental evidence, and the confidence score cutoff is set to 900. In addition, to observe the correlations between our differential expression analysis data and other aging-related proteomics results in external databases, we compared the DEP data with data from the Human Aging Gene Database (https://genomics.senescence.info/genes/index.html) (last updated on February 9, 2020) and from plasma proteomics studies related to aging changes in healthy populations [29].Table 2
Functional role of the 21 proteins.
Protein IDUniProtGene symbolFull nameFunctionHBBP68871HBBHemoglobin subunit betaThis protein includesiron ion binding and oxygen bindingHBAP69905HBA1, HBA2Hemoglobin subunit alphaInvolved in oxygen transport from the lung to the various peripheral tissuesAPOA2P02652APOA2Apolipoprotein A-IIMay stabilize the HDL (high-density lipoprotein) structure by its association with lipids and affect the HDL metabolismVAT1Q99536VAT1Synaptic vesicle membrane proteinPossesses ATPase activity (by similarity) and plays a part in calcium-regulated keratinocyte activation in epidermal repair mechanismsLG3BPQ08380LGALS3BPGalectin-3-binding proteinPromotes integrin-mediated cell adhesion and may stimulate host defense against viruses and tumor cellsKRT35Q92764KRT35Keratin, type I cuticular Ha5A member of the keratin gene familyABRX2Q15018ABRAXAS2/FAM175BBRISC complex subunit Abraxas 2Required for normal induction of p53/TP53 in response to DNA damageBAFO75531BANF1Barrier-to-autointegration factorPlays fundamental roles in nuclear assembly, chromatin organization, gene expression, and gonad developmentG3PP04406GAPDHGlyceraldehyde-3-phosphate dehydrogenaseThe product of this gene catalyzes an important energy-yielding step in carbohydrate metabolismAPOA1P02647APOA1Apolipoprotein A-IThis protein is associated with the transport of cholesterol and the major protein component of high-density lipoprotein (HDL)HPTP00738HPHaptoglobinThis protein produces haptoglobin, which encodes a protein with antimicrobial activity against bacteriaCLUSP10909CLUClusterinThis protein is involved in several basic biological events such as cell death, tumor progression, and neurodegenerative disordersA1AG1P02763ORM1Alpha-1-acid glycoprotein 1This protein is classified as an acute-phase reactant and may be involved in aspects of immunosuppressionA2MGP01023A2MAlpha-2-macroglobulinThis gene can inhibit inflammatory cytokines and thus disrupts inflammatory cascadesCYTCP01034CST3Cystatin-CThis protein was shown to have an antimicrobial functionFIBGP02679FGGFibrinogen gamma chainThe protein encoded by this gene is the gamma component of fibrinogenCERUP00450CPCeruloplasminA glycoprotein with ferroxidase activityFIBAP02671FGAFibrinogen alpha chainThis gene encodes the alpha subunit of the coagulation factor fibrinogenANT3P01008SERPINC1Antithrombin-IIIA member of the serpin superfamily and a plasma protease inhibitorFIBBP02675FGBFibrinogen beta chainThe protein encoded by this gene is the beta component of fibrinogenVTDBP02774GCVitamin D-binding proteinInvolved in vitamin D transport and storage and scavenging of extracellular G-actin


### 2.4. Independent Samplet-Tests for within-Age Group Differences

Information was compiled on the DEPs identified in the age groups in the DIA experiment. To reduce error caused by NA values in the experimental data, SPSS 25.0 software was used to perform independent samplet-tests for the different age groups. Only when a p<0.05 between age groups was obtained were independent sample t-tests performed between the male and female subgroups within the age groups. Both the p value and the average value of protein expression were recorded for the different sexes. GraphPad Prism 8.3.0 was used to create violin plots of the proteins with significant differences in the t-tests in order to better visualize the distribution of the data and the differences in expression between sexes. GraphPad was also used to generate receiver operating characteristic (ROC) curves and heatmaps for selected DEPs.


## 2.1. Study Participants

Skin samples were collected from 9 healthy Chinese males and 11 healthy Chinese females, none of whom had experienced excessive light exposure. The participants were divided into two groups according to age: the 20–30 y group (3 males and 7 females) and the 55–70 y group (6 males and 4 females). Individuals were excluded from the study if they had an allergy to tape, were long-term outdoor workers, or had a skin disease or other systemic condition involving skin problems. The participants did not use topical emollients or other cosmetic products for 24 hours prior to the experiment. This study was conducted in accordance with the recommendations of the Medical Ethics Committee of Anhui Medical University, and informed consent was obtained from all enrolled subjects.


### 2.1.1. Materials

Sodium dodecyl sulfate (SDS), L3 lysate without SDS, trypsin, dithiothreitol (DTT), Coomassie Brilliant blue G-250, 40% ethanol, 10% acetic acid, ethylenediaminetetraacetic acid (EDTA), and acrylonitrile (ACN) were obtained from Wallis Technology, Beijing, China. Ultrafiltration membranes (10K MWCO, 1.5 ml, plate) were obtained from Pall Corp. (NY, USA). Pierce C18 pipette tips (10μl bed), Empore™ C18 47 mm extraction discs (model 2215), a ThermoMixer (MS-100), and a CentriVap vacuum concentrator and accessories were obtained from Thermo Fisher Scientific (Shanghai, China). The samples were analyzed using a Q Exactive high-frequency mass spectrometer and an UltiMate 3000 high-performance liquid chromatography system (Thermo Fisher Scientific, San Jose, CA, USA).


### 2.1.2. Sample Preparation

The palm side of the forearm was gently wiped with a sterile cotton ball to cleanse it of epidermal contaminants. A3.0∗3.0cm piece of a 3M medical tape was prepared to fully and precisely cover the back of the slide. The slide was then moved back and forth quickly with even pressure for 2 minutes. Even pressure was applied so that the slide was covered. Once the sample was collected, another slide was placed over it to protect the skin sample. To reduce sampling variation, the same technician collected samples from all volunteers during the study.


### 2.1.3. Protein Extraction and Enzymolysis

(1) Protein Extraction. Small pieces (0.5cm∗0.5cm) of tape/skin samples were cut out of the slides with a sterile blade and then placed into corresponding 1.5 ml centrifuge tubes according to the sample numbers. An appropriate amount of L3 lysis buffer without SDS was added, and cocktail containing EDTA was added to a final concentration of ×1. The samples were placed on ice for 5 min, and DTT was added to a final concentration of 10 mM. The samples were soaked overnight and then centrifuged at 25,000×g and 4°C for 15 min. The supernatant was obtained, and DTT was added to a final concentration of 10 mM. The samples were incubated in a 56°C water bath for 1 hour, after which iodoacetamide (IAM) was added to a final concentration of 55 mM. The samples were placed in a dark room for 45 min and then centrifuged at 25,000×g at 4°C for 15 min to obtain the supernatant, which contained the proteins.(2) Protein Extraction and Quality Control. The protein concentrations were measured using the Bradford method [23] with acceptable enzymatic efficiency, and it was determined that the extracted protein content was sufficient. For the Bradford quantitative assay, protein standards (0.2 μg/μl BSA; 0, 2, 4, 6, 8, 10, 12, 14, 16, and 18 μl) were added to wells A1 through A10 of 96-well MICROLON ELISA plates, and then, pure water (20, 18, 16, 14, 12, 10, 8, 6, 4, and 2 μl) was added. Then, 180 μl of Coomassie Brilliant blue G-250 quantitative work solution was added to each well. A linear standard curve was prepared based on the optical density at 595 nm (OD595) and the protein concentration. The protein sample solutions were diluted several times, 180 μl of the quantitative solution was added to 20 μl of protein solution, and the OD595 was read. The protein concentration of each sample was calculated according to the standard curve and the sample OD595. The purity of the extracted proteins was verified by SDS-PAGE and Coomassie Brilliant blue staining.(3) SDS-PAGE. For each sample, 30 μg of protein was added to an appropriate amount of loading buffer. After thorough mixing, the proteins in buffer were heated at 95°C for 5 min, centrifuged at 25,000×g for 5 min, and placed into the wells of a 12% SDS polyacrylamide gel. After electrophoresis, Coomassie Brilliant blue staining was conducted for 2 hours. An appropriate amount of decolorization solution (40% methanol and 10% acetic acid) was added to each sample, and the samples with decolorization solution were placed in a shaker 3 to 5 times for 30 min each.(4) Proteolysis. Trypsin (2.5 μg) was added to 100 μg of protein for each sample at a protein : enzyme ratio of 40 : 1, and the protein was hydrolyzed at 37°C for 4 hours. Trypsin was added once more according to the above ratio, and enzymolysis was continued for 8 hours at 37°C. The enzymatically digested peptides were desalted with a Strata-X column and extracted under vacuum. With this method, we extracted sufficient amounts of proteins with qualified enzymatic efficiency.


### 2.1.4. High-pH Reverse-Phase Separation

A pooled sample was created from 10μg of each sample, and 200 μg of the pooled sample was mixed with 2 ml of mobile phase A (5% acetonitrile, ACN, pH 9.8) and diluted into a Shimadzu LC-20AB liquid chromatography system. The samples were separated in a liquid phase on a 5 μm 4.6×250mm Gemini C18 column. Elution was conducted at a flow rate of 1 ml/min with the following gradient: 5% mobile phase B (95% ACN, pH 9.8) for 10 min, 5% to 35% mobile phase B for 40 min, 35% to 95% mobile phase B for 1 min, mobile phase B which lasted for 3 min, and equilibration at 5% mobile phase B for 10 min. The elution peaks were monitored at 214 nm, and one fraction was collected every minute. The samples were combined with the chromatogram elution peaks to obtain 10 fractions, which were then extracted by freezing.


### 2.1.5. High-Performance Liquid Chromatography

The extracted peptide samples were redissolved (centrifuged with mobile phase A (2% ACN and 0.1% formic acid (FA)) at20,000×g for 10 min), and the supernatant liquid was sampled. Separation was performed with a Thermo UltiMate 3000 UHPLC. Each sample was first injected into a trap column for enrichment and desalting before being run through a self-loading C18 column (150 μm I.D., 1.8 μm pore size, 25 cm column length) coupled to tandem mass spectrometers. The samples were separated at a flow rate of 500 nl/min through the following effective gradient: 0–5 min, 5% mobile phase B (98% ACN, 0.1% FA); 5–160 min, linear increase in mobile phase B from 5% to 35%; 160–170 min, 35% to 80% mobile phase B; 170–175 min, 80% mobile phase B; and 176–180 min, 5% mobile phase B. The output of the liquid chromatograph was connected directly to the mass spectrometer.


### 2.1.6. Data-Dependent Acquisition (DDA) and Data-Independent Acquisition (DIA) Mass Spectrometry Analyses

After liquid chromatography, the peptides were separated by nanoelectrospray ionization (ESI) and detected in a DDA mode with a Q Exactive HF (Thermo Fisher Scientific, San Jose, CA) tandem mass spectrometer. The main parameter settings were as follows: the ion source voltage was 1.6 kV, the primary mass spectrometer range was 350~1,500 m/z, the resolution was 60,000, the MS2 starting m/z was fixed at 100, and the resolution of the second-stage mass spectrometer was 15,000. The screening criteria for the precursor ions in the MS2 scans were a charge of 2+ to 7+ with a peak intensity of over 10,000 and a ranking in the top 20. The main parameter settings for DIA mode detection were as follows: the ion source voltage was 1.6 kV, the MS1 scanning range was 350–1,500 m/z with a resolution of 120,000, and the 350–1,500 Da range was divided into 40 windows for fragmentation and signal acquisition. Higher-energy collisional dissociation (HCD) was used for the DDA and DIA modes of mass spectrometric detection of ion fragmentation patterns. The fragmented ions were detected in an Orbitrap. The dynamic exclusion time was set to 30 s. The automatic gain control (AGC) settings were as follows: level 1: 3E6 and level 2: 1E5.


## 2.1.1. Materials

Sodium dodecyl sulfate (SDS), L3 lysate without SDS, trypsin, dithiothreitol (DTT), Coomassie Brilliant blue G-250, 40% ethanol, 10% acetic acid, ethylenediaminetetraacetic acid (EDTA), and acrylonitrile (ACN) were obtained from Wallis Technology, Beijing, China. Ultrafiltration membranes (10K MWCO, 1.5 ml, plate) were obtained from Pall Corp. (NY, USA). Pierce C18 pipette tips (10μl bed), Empore™ C18 47 mm extraction discs (model 2215), a ThermoMixer (MS-100), and a CentriVap vacuum concentrator and accessories were obtained from Thermo Fisher Scientific (Shanghai, China). The samples were analyzed using a Q Exactive high-frequency mass spectrometer and an UltiMate 3000 high-performance liquid chromatography system (Thermo Fisher Scientific, San Jose, CA, USA).


## 2.1.2. Sample Preparation

The palm side of the forearm was gently wiped with a sterile cotton ball to cleanse it of epidermal contaminants. A3.0∗3.0cm piece of a 3M medical tape was prepared to fully and precisely cover the back of the slide. The slide was then moved back and forth quickly with even pressure for 2 minutes. Even pressure was applied so that the slide was covered. Once the sample was collected, another slide was placed over it to protect the skin sample. To reduce sampling variation, the same technician collected samples from all volunteers during the study.


## 2.1.3. Protein Extraction and Enzymolysis

(1) Protein Extraction. Small pieces (0.5cm∗0.5cm) of tape/skin samples were cut out of the slides with a sterile blade and then placed into corresponding 1.5 ml centrifuge tubes according to the sample numbers. An appropriate amount of L3 lysis buffer without SDS was added, and cocktail containing EDTA was added to a final concentration of ×1. The samples were placed on ice for 5 min, and DTT was added to a final concentration of 10 mM. The samples were soaked overnight and then centrifuged at 25,000×g and 4°C for 15 min. The supernatant was obtained, and DTT was added to a final concentration of 10 mM. The samples were incubated in a 56°C water bath for 1 hour, after which iodoacetamide (IAM) was added to a final concentration of 55 mM. The samples were placed in a dark room for 45 min and then centrifuged at 25,000×g at 4°C for 15 min to obtain the supernatant, which contained the proteins.(2) Protein Extraction and Quality Control. The protein concentrations were measured using the Bradford method [23] with acceptable enzymatic efficiency, and it was determined that the extracted protein content was sufficient. For the Bradford quantitative assay, protein standards (0.2 μg/μl BSA; 0, 2, 4, 6, 8, 10, 12, 14, 16, and 18 μl) were added to wells A1 through A10 of 96-well MICROLON ELISA plates, and then, pure water (20, 18, 16, 14, 12, 10, 8, 6, 4, and 2 μl) was added. Then, 180 μl of Coomassie Brilliant blue G-250 quantitative work solution was added to each well. A linear standard curve was prepared based on the optical density at 595 nm (OD595) and the protein concentration. The protein sample solutions were diluted several times, 180 μl of the quantitative solution was added to 20 μl of protein solution, and the OD595 was read. The protein concentration of each sample was calculated according to the standard curve and the sample OD595. The purity of the extracted proteins was verified by SDS-PAGE and Coomassie Brilliant blue staining.(3) SDS-PAGE. For each sample, 30 μg of protein was added to an appropriate amount of loading buffer. After thorough mixing, the proteins in buffer were heated at 95°C for 5 min, centrifuged at 25,000×g for 5 min, and placed into the wells of a 12% SDS polyacrylamide gel. After electrophoresis, Coomassie Brilliant blue staining was conducted for 2 hours. An appropriate amount of decolorization solution (40% methanol and 10% acetic acid) was added to each sample, and the samples with decolorization solution were placed in a shaker 3 to 5 times for 30 min each.(4) Proteolysis. Trypsin (2.5 μg) was added to 100 μg of protein for each sample at a protein : enzyme ratio of 40 : 1, and the protein was hydrolyzed at 37°C for 4 hours. Trypsin was added once more according to the above ratio, and enzymolysis was continued for 8 hours at 37°C. The enzymatically digested peptides were desalted with a Strata-X column and extracted under vacuum. With this method, we extracted sufficient amounts of proteins with qualified enzymatic efficiency.


## 2.1.4. High-pH Reverse-Phase Separation

A pooled sample was created from 10μg of each sample, and 200 μg of the pooled sample was mixed with 2 ml of mobile phase A (5% acetonitrile, ACN, pH 9.8) and diluted into a Shimadzu LC-20AB liquid chromatography system. The samples were separated in a liquid phase on a 5 μm 4.6×250mm Gemini C18 column. Elution was conducted at a flow rate of 1 ml/min with the following gradient: 5% mobile phase B (95% ACN, pH 9.8) for 10 min, 5% to 35% mobile phase B for 40 min, 35% to 95% mobile phase B for 1 min, mobile phase B which lasted for 3 min, and equilibration at 5% mobile phase B for 10 min. The elution peaks were monitored at 214 nm, and one fraction was collected every minute. The samples were combined with the chromatogram elution peaks to obtain 10 fractions, which were then extracted by freezing.


## 2.1.5. High-Performance Liquid Chromatography

The extracted peptide samples were redissolved (centrifuged with mobile phase A (2% ACN and 0.1% formic acid (FA)) at20,000×g for 10 min), and the supernatant liquid was sampled. Separation was performed with a Thermo UltiMate 3000 UHPLC. Each sample was first injected into a trap column for enrichment and desalting before being run through a self-loading C18 column (150 μm I.D., 1.8 μm pore size, 25 cm column length) coupled to tandem mass spectrometers. The samples were separated at a flow rate of 500 nl/min through the following effective gradient: 0–5 min, 5% mobile phase B (98% ACN, 0.1% FA); 5–160 min, linear increase in mobile phase B from 5% to 35%; 160–170 min, 35% to 80% mobile phase B; 170–175 min, 80% mobile phase B; and 176–180 min, 5% mobile phase B. The output of the liquid chromatograph was connected directly to the mass spectrometer.


## 2.1.6. Data-Dependent Acquisition (DDA) and Data-Independent Acquisition (DIA) Mass Spectrometry Analyses

After liquid chromatography, the peptides were separated by nanoelectrospray ionization (ESI) and detected in a DDA mode with a Q Exactive HF (Thermo Fisher Scientific, San Jose, CA) tandem mass spectrometer. The main parameter settings were as follows: the ion source voltage was 1.6 kV, the primary mass spectrometer range was 350~1,500 m/z, the resolution was 60,000, the MS2 starting m/z was fixed at 100, and the resolution of the second-stage mass spectrometer was 15,000. The screening criteria for the precursor ions in the MS2 scans were a charge of 2+ to 7+ with a peak intensity of over 10,000 and a ranking in the top 20. The main parameter settings for DIA mode detection were as follows: the ion source voltage was 1.6 kV, the MS1 scanning range was 350–1,500 m/z with a resolution of 120,000, and the 350–1,500 Da range was divided into 40 windows for fragmentation and signal acquisition. Higher-energy collisional dissociation (HCD) was used for the DDA and DIA modes of mass spectrometric detection of ion fragmentation patterns. The fragmented ions were detected in an Orbitrap. The dynamic exclusion time was set to 30 s. The automatic gain control (AGC) settings were as follows: level 1: 3E6 and level 2: 1E5.


## 2.2. Analysis of Epidermal Proteomics Data

The DDA data from the machine were identified using the Andromeda search engine integrated with MaxQuant [24], and Spectronaut was then used to build a spectral library with the results. The resulting data were reviewed using the UniProtKB/Swiss-Prot Homo sapiens proteome database. For large-scale DIA data, after constructing the spectral library information, convolutional extraction of the data was completed using Spectronaut, and data analysis and quality control were completed using the mProphet algorithm, resulting in reliable protein quantification results. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway functional annotations were performed. Based on the high-quality quantitative results, we searched for differentially expressed proteins (DEPs) between groups.MaxQuant was used to complete the identification of the DDA data in order to create the spectral library for subsequent DIA data analysis. During this procedure, the original offline data were used as the input data, the corresponding parameters and database were configured, and then, the identification and quantitative analysis were performed. The identification information that met the criterion of a false discovery rateFDR≤1% was used to build the final spectral library. The offline DIA data were analyzed using Spectronaut [25], and iRT peptides were used to correct the retention time. Spectronaut integrates the mProphet scoring algorithm, which can accurately reflect the degree of matching of isolated pairs. Then, based on the target-decoy model applied by SWATH-MS, the false-positive rate was controlled at 1% FDR, which yielded significant quantitative results.This procedure was used to preprocess the data according to the set comparison group, and a significance test was then performed based on the model. After that, the DEPs were screened with the criteria of a fold change (FC)value≥1.5 and a p value < 0.05 (the top 10 DEPs are listed in Table 1). The R package MSstats [26] from the Bioconductor resource library was used to complete differential analysis of the DEPs.Table 1
Top 10 differentially expressed proteins.
UpregulatedDownregulatedProtein IDlog2FCp valueProtein IDlog2FCp valueOld vs youngLG3BP0.7000.001ABRX2−0.7500.001HBA3.6060.001PLD3−0.6980.003RENR1.0140.001VTI1B−0.5870.004APOA23.0690.001VAT1−0.7830.004CLUS1.6670.001KRT35−1.6450.004SG1D22.6610.002BAF−2.4920.005SEM7A1.4140.002ECM1−0.6310.009B2MG1.5900.002K1C15−0.6000.010S10A40.7340.002MYH9−0.9070.011HBB3.4620.003PSB7−0.6780.012


## 2.3. Screening for DEPs

After the raw data were median-normalized by the software, the resulting data were selected for follow-up analysis. DEPs were identified through symmetrical scatter plot screening. The dots represent the relative mean expression values (normalized andlog2 transformed) of the proteins in the two groups, and the log2FC values were calculated using the R package MSstats. The dashed line in the graph represents the threshold line at ∣log2FC∣=1.5. Thus, the volcano plot of DEPs was produced with the average protein expression values, log2FC values, and p values. The color of the dot for each protein indicates the significance of the p value, with the blue to red gradient indicating p values ranging from nonsignificant to significant. The graph shows information about the gradient change. We performed protein-protein interaction network analysis on all the DEPs and screened the proteins with more than 10 nodes of neighboring proteins in the protein relationship network. The proteins with the most significant differences in log2FC values and p values were also selected (Table 2). A gene-level network analysis was conducted for the selected most promising candidate proteins, and the biological pathways in which they were enriched were investigated by using NetworkAnalyst 3.0 inline tool [27, 28], the value of the interaction relationship between genes retained the interaction relationship with experimental evidence, and the confidence score cutoff is set to 900. In addition, to observe the correlations between our differential expression analysis data and other aging-related proteomics results in external databases, we compared the DEP data with data from the Human Aging Gene Database (https://genomics.senescence.info/genes/index.html) (last updated on February 9, 2020) and from plasma proteomics studies related to aging changes in healthy populations [29].Table 2
Functional role of the 21 proteins.
Protein IDUniProtGene symbolFull nameFunctionHBBP68871HBBHemoglobin subunit betaThis protein includesiron ion binding and oxygen bindingHBAP69905HBA1, HBA2Hemoglobin subunit alphaInvolved in oxygen transport from the lung to the various peripheral tissuesAPOA2P02652APOA2Apolipoprotein A-IIMay stabilize the HDL (high-density lipoprotein) structure by its association with lipids and affect the HDL metabolismVAT1Q99536VAT1Synaptic vesicle membrane proteinPossesses ATPase activity (by similarity) and plays a part in calcium-regulated keratinocyte activation in epidermal repair mechanismsLG3BPQ08380LGALS3BPGalectin-3-binding proteinPromotes integrin-mediated cell adhesion and may stimulate host defense against viruses and tumor cellsKRT35Q92764KRT35Keratin, type I cuticular Ha5A member of the keratin gene familyABRX2Q15018ABRAXAS2/FAM175BBRISC complex subunit Abraxas 2Required for normal induction of p53/TP53 in response to DNA damageBAFO75531BANF1Barrier-to-autointegration factorPlays fundamental roles in nuclear assembly, chromatin organization, gene expression, and gonad developmentG3PP04406GAPDHGlyceraldehyde-3-phosphate dehydrogenaseThe product of this gene catalyzes an important energy-yielding step in carbohydrate metabolismAPOA1P02647APOA1Apolipoprotein A-IThis protein is associated with the transport of cholesterol and the major protein component of high-density lipoprotein (HDL)HPTP00738HPHaptoglobinThis protein produces haptoglobin, which encodes a protein with antimicrobial activity against bacteriaCLUSP10909CLUClusterinThis protein is involved in several basic biological events such as cell death, tumor progression, and neurodegenerative disordersA1AG1P02763ORM1Alpha-1-acid glycoprotein 1This protein is classified as an acute-phase reactant and may be involved in aspects of immunosuppressionA2MGP01023A2MAlpha-2-macroglobulinThis gene can inhibit inflammatory cytokines and thus disrupts inflammatory cascadesCYTCP01034CST3Cystatin-CThis protein was shown to have an antimicrobial functionFIBGP02679FGGFibrinogen gamma chainThe protein encoded by this gene is the gamma component of fibrinogenCERUP00450CPCeruloplasminA glycoprotein with ferroxidase activityFIBAP02671FGAFibrinogen alpha chainThis gene encodes the alpha subunit of the coagulation factor fibrinogenANT3P01008SERPINC1Antithrombin-IIIA member of the serpin superfamily and a plasma protease inhibitorFIBBP02675FGBFibrinogen beta chainThe protein encoded by this gene is the beta component of fibrinogenVTDBP02774GCVitamin D-binding proteinInvolved in vitamin D transport and storage and scavenging of extracellular G-actin


## 2.4. Independent Samplet-Tests for within-Age Group Differences

Information was compiled on the DEPs identified in the age groups in the DIA experiment. To reduce error caused by NA values in the experimental data, SPSS 25.0 software was used to perform independent samplet-tests for the different age groups. Only when a p<0.05 between age groups was obtained were independent sample t-tests performed between the male and female subgroups within the age groups. Both the p value and the average value of protein expression were recorded for the different sexes. GraphPad Prism 8.3.0 was used to create violin plots of the proteins with significant differences in the t-tests in order to better visualize the distribution of the data and the differences in expression between sexes. GraphPad was also used to generate receiver operating characteristic (ROC) curves and heatmaps for selected DEPs.


## 3. Results

### 3.1. Protein Quantification and Data Quality Control

In this project, mass spectrometry data were collected for 20 samples using a Q Exactive HF instrument in the DIA mode. Peptide quantification and protein quantification were completed using Spectronaut and MSstats software. The specific quantitative information for each sample is shown in Figures1(d) and 1(e). In total, we identified 1,318 proteins in the 20 samples, which were further analyzed with MSstats software, resulting in 1,270 proteins. The quality of the data was assessed by analysis of intragroup coefficients of variation, principal component analysis, and quantitative sample correlation analysis. We calculated Pearson’s correlation coefficients for the expression levels of all protein expressions between the two groups and displayed these coefficients as a heat map, as shown in Figure 2(b). We also acquired DDA mass spectrometry data from the samples and then used MaxQuant to complete library identification in order to obtain nonredundant, high-quality MS/MS spectral information for subsequent quantification in the DIA mode. The spectral library contained fragment ion intensities and retention times, which characterized the peptide spectral peaks. With regard to statistical information for the peptides and proteins in the spectral library, we obtained the peptide distribution, protein amount distribution, and protein coverage distribution (Figures 1(a)–1(c)).Figure 1
(a) Identification of DDA profiles for subsequent DIA quantification: unique peptide distribution with the horizontal axis being the only matching peptide number per protein and the vertical axis being the number of proteins; (b) protein mass distribution with the horizontal axis being the protein mass interval (Kilodalton) and the vertical axis being the corresponding protein number; (c) protein coverage distribution with the horizontal axis being the percentage protein coverage interval and the vertical axis being the number of proteins. (d, e) Unique proteins and peptides were counted in each sample collected by Q Exactive HF in the DIA mode.
(a)(b)(c)(d)(e)Figure 2
(a) Flowchart of protein extraction, enzymatic hydrolysis experiments, and information analysis. (b) Thex-axis and y-axis of both samples. Color represents the correlation coefficient (the bluer the color, the higher the correlation and the lighter the color, the lower the correlation).
(a)(b)


### 3.2. Identification of 95 DEPs in Aging via Proteomics Analysis

In order to obtain quantitative proteomic maps of the stratum corneum of keratinocytes derived from elderly and young populations, we used data extracted from MaxQuant and Spectronaut convolutions and applied bioinformatics analysis using 1% FDR analysis to identify 95 unique proteins that were differentially expressed between the age groups. The supplementary tables show important information about the identified proteins, including their UniProtKB accession numbers, IDs, protein and gene names, adjustedp values, and associated functions. The proteins that were significantly differentially expressed according to age included 38 downregulated proteins and 57 upregulated proteins (Supplementary Table 1 and Supplementary Table 2). In an interaction analysis of the 95 DEPs, we counted the proteins with more than 10 nodes of protein interactions (Figure 3(a)). We compared our DEPs with 307 human aging-related genes from the GenAge database and identified 6 genes that intersected with the DEPs, namely, CDC42, CLU, A2M, LMNA, VCP, and AIFM1 (Figure 4(c)). In a proteomics study on the senescence-associated secretory phenotype [30], serine protease inhibitors (SERPINs) were shown to be significantly associated with age changes in human cohorts and to be potential biomarkers. In the induced aging model in that study, the secretion levels of the proteins SERPINE1, SERPINF1, SERPINF2, SERPING1, and SERPINH1 were significantly increased. In contrast, SERPINA3, SERPINB1, SERPINB6, and SERPINC1 were significantly upregulated with age in our study. In an age-associated proteomics study in a healthy population [29], 217 age-associated proteins were present, 9 of which intersect with the DEPs in our study: ANXA2, S100A4, B2M, FGA, FGB, FGG, SERPINA3, CST3, and HP (Figure 4). In the previous study, SERPINA3 and SERPING1 were positively correlated with age, while SERPINF2 was negatively correlated with age, which is consistent with our results.Figure 3
(a) Proteins possessing more than 10 neighboring proteins in the protein-protein interaction network analysis. (b–e) Network analysis of gene expression was performed for the 21 selected genes, where the genes we selected were indicated by red circles. The pathway with the most genes involved in (b) is complement and coagulation cascades; the pathway with the most genes involved in (c) is cardiomyopathy; the pathway with the most genes involved in (d) is ferroptosis; and the pathway with the most genes involved in (e) is apoptosis and autophagy. (Supplementary Tables 5, 6, 7, and 8 for details).
(a)(b)(c)(d)(e)Figure 4
(a) Symmetric scatter plots compare the differential status of proteins in the age group. The horizontal and vertical coordinates represent the mean relative expression of young and old groups; the dashed line in the graph represents the threshold line at∣log2FC∣=1.5. The color of the scatter represents the p value. (b) Wayne plots show the distribution of protein expression between the young and old groups, with coverage areas representing differentially expressed proteins. (c) Wayne diagram of our proteomic results and 307 genes from the GenAge Human Aging Gene Database. (d) Venn diagram of our proteomic results and an aging-related plasma proteomics study in a healthy population.
(a)(b)(c)(d)


### 3.3. Functional Analysis of the DEPs

Through DIA mass spectrometry, we identified a total of 1,270 unique proteins, of which 95 were differentially expressed in the different age groups. The 95 identified proteins were analyzed using the GO and KEGG databases (the functional classifications of the 95 DEPs are shown in Supplementary Table3). The functional pathways into which the DEPs were classified are shown in Figure 5(a) including the cellular process, environmental information processing, genetic information processing, human disease, metabolism, and organismal system pathways. We analyzed the DEPs with the KEGG and found that the most enriched pathways included platelet activation, complement and coagulation cascades, and the VEGF signaling pathway. Among the DEPs, BANF1 was found to be responsible for the chromatin structure and dynamics, replication, recombination, and repair. Energy production and transfer are also entry points for the aging process, and LDHB and HBA1 were enriched in these pathways. The signaling pathways that we focused on in relation to aging included ferroptosis and apoptosis, both of which are cell growth- and death-related metabolic pathways. Among the six proteins intersecting with GenAge data, AIFM1 and LMNA are involved in apoptosis. In the gene-level network analysis of the 21 most promising proteins we screened (Figures 3(b)–3(e)), the significantly enriched pathways were the complement and coagulation cascades, platelet activation, ferroptosis, apoptosis, autophagy, and PPAR signaling pathways. Among them, the PPAR signaling pathway may be closely related to the aging process [31]. (For more detailed gene node information, please refer to Supplementary Tables 4, 5, 6, 7, and 8 and the following website: https://www.networkanalyst.ca/NetworkAnalyst/faces/Share?ID=_9qljssi8x.) A protein of particular interest was ceruloplasmin (CP), which had 16 neighboring proteins in the protein-protein interaction network. CP is closely associated with aging [32, 33], and the biological metabolic pathway in which it participates is the well-known ferroptosis pathway (Figure 6), a signaling pathway that is closely associated with aging. Our age-associated proteomics results are consistent with the results obtained in other age-associated proteomics studies with independent cohorts (Supplementary Tables 1 and 2 at 10.1038/s41591-019-0673-2).Figure 5
(a) Pathway classification of differentially expressed proteins, thex-axis represents the number of differentially expressed proteins, the y-axis represents the pathway-annotated entries. (b) Pathway relationships; red and blue dots indicate up- and downregulated differential proteins, respectively. Purple dots indicate the top ten enriched pathways, darker dots indicate significant enrichment, lighter dots indicate enrichment but not significant, and larger areas indicate higher enrichment. The different color lines indicate the different pathway categories: blue for environmental information processing, purple for human disease, orange for metabolism, and yellow for organismal systems.
(a)(b)Figure 6
Changes in proteins in the Ferroptosis signaling pathway with aging. CP is a protein that is upregulated with aging and also highlights important gender differences with aging, and it plays an important role in this well-known cellular aging and death pathway.


### 3.4. Potential Candidate Biomarkers of Sex Differences in the Aging Process

We compared the proteins that changed most significantly under the influence of two factors, age and sex, according to a dual statistical threshold defined by thep value and log2FC. The proteins that were most significantly upregulated in the elderly group were LGALS3BP, HBA1, and HBB, while the proteins that were most significantly upregulated in the young group were ABRAXAS2, BANF1, VAT1, and KRT35 (Figure 4(a)). The proteins in the graph are colored according to the significance of the p value, and the graph clearly shows that a smaller p value was associated with a greater difference in protein expression between the groups. Women generally live longer than men [18], and aging rates differ between men and women [34]. From the 95 proteins identified in our quantitative proteomics experiments and through statistical analysis via independent sample paired t-tests between the age and sex groups, we selected proteins as candidate biomarkers for distinguishing differences in the aging process between men and women (Figure 7). The ROC curves of several proteins of interest are shown in Figure 8, and the corresponding area under the curve (AUC) values and p values were calculated. ABRAXAS2, VAT1, and LGALS3BP were identified as promising candidate biomarkers of aging. t-tests for these three proteins revealed a p value < 0.05 between men and women in the younger group and a p value > 0.05 between men and women in the older group (Table 3). Thus, the three proteins were expressed differently between sex groups in the young population, but this difference disappeared with aging. Among these proteins, ABRAXAS2 and VAT1 were downregulated in the elderly group and LGALS3BP was upregulated in the elderly group. The results of tandem mass spectrometry of the three proteins are shown in Figures 9(a)–9(c).Figure 7
(a–c) Violin plots of the relative expression abundance of three proteins, ABRAXAS2, LGALS3BP, and VAT1, in the young-old group. The red solid line is the mean value and the blue dashed line is the quartile. The data in the vertical coordinates are protein expression abundance values taken aslog2. (d) Heatmap of six proteins of a particular value. (e) Mean values of the expression abundance of the three proteins in young and old males and females. The p values between the groups are marked. It can be seen that the three proteins differ between young men and women and do not differ between older men and women.
(a)(b)(c)(d)(e)Figure 8
(a–f) Receiver operating characteristic curve (ROC curve) with an AUC andp value obtained by comparing six proteins in young and older people.
(a)(b)(c)(d)(e)(f)Table 3
Different protein expressions at different ages in different genders. Thet-test showed significant differences among age groups (∗p<0.05) and between genders for ABRAXAS2, VAT1, and LGALS3BP (∗∗p<0.05). More details are described in Section 3. The case group is an older group (55–70 years) and the control group is a young group (20–30 years). M: mean protein expression of male; F: mean protein expression of female.
AccessionIDDescriptionGeneOld and young differencesOld group differences in gendersYoung group differences in gendersDownQ15018ABRX2_HUMANBRISC complex subunit Abraxas 2ABRAXAS2∗p=0.001p=0.302M=13.45761365F=13.10540960∗∗p=0.016M=13.6814853567F=14.2325095143Q99536VAT1_HUMANSynaptic vesicle membrane protein VAT-1 homologVAT1∗p=0.004p=0.729M=15.87932308F=15.78413334∗∗p=0.038M=17.2356184533F=16.3622206371UpQ08380LG3BP_HUMANGalectin-3-binding proteinLGALS3BP∗p=0.001p=0.549M=15.20773285F=15.01199494∗∗p=0.034M=14.8493189733F=14.2752645800Figure 9
(a) The sequence NH2-AIYQVYNALQEK-COOH allows the identification of ABRAXAS2. (b) The sequence NH2-AVDTWSWGER-COOH allows the identification of LGALS3BP. (c) The sequence NH2-TEAASDPQHPAASEGAAAAAASPPLLR-COOH allows the identification of VAT1.
(a)(b)(c)


## 3.1. Protein Quantification and Data Quality Control

In this project, mass spectrometry data were collected for 20 samples using a Q Exactive HF instrument in the DIA mode. Peptide quantification and protein quantification were completed using Spectronaut and MSstats software. The specific quantitative information for each sample is shown in Figures1(d) and 1(e). In total, we identified 1,318 proteins in the 20 samples, which were further analyzed with MSstats software, resulting in 1,270 proteins. The quality of the data was assessed by analysis of intragroup coefficients of variation, principal component analysis, and quantitative sample correlation analysis. We calculated Pearson’s correlation coefficients for the expression levels of all protein expressions between the two groups and displayed these coefficients as a heat map, as shown in Figure 2(b). We also acquired DDA mass spectrometry data from the samples and then used MaxQuant to complete library identification in order to obtain nonredundant, high-quality MS/MS spectral information for subsequent quantification in the DIA mode. The spectral library contained fragment ion intensities and retention times, which characterized the peptide spectral peaks. With regard to statistical information for the peptides and proteins in the spectral library, we obtained the peptide distribution, protein amount distribution, and protein coverage distribution (Figures 1(a)–1(c)).Figure 1
(a) Identification of DDA profiles for subsequent DIA quantification: unique peptide distribution with the horizontal axis being the only matching peptide number per protein and the vertical axis being the number of proteins; (b) protein mass distribution with the horizontal axis being the protein mass interval (Kilodalton) and the vertical axis being the corresponding protein number; (c) protein coverage distribution with the horizontal axis being the percentage protein coverage interval and the vertical axis being the number of proteins. (d, e) Unique proteins and peptides were counted in each sample collected by Q Exactive HF in the DIA mode.
(a)(b)(c)(d)(e)Figure 2
(a) Flowchart of protein extraction, enzymatic hydrolysis experiments, and information analysis. (b) Thex-axis and y-axis of both samples. Color represents the correlation coefficient (the bluer the color, the higher the correlation and the lighter the color, the lower the correlation).
(a)(b)


## 3.2. Identification of 95 DEPs in Aging via Proteomics Analysis

In order to obtain quantitative proteomic maps of the stratum corneum of keratinocytes derived from elderly and young populations, we used data extracted from MaxQuant and Spectronaut convolutions and applied bioinformatics analysis using 1% FDR analysis to identify 95 unique proteins that were differentially expressed between the age groups. The supplementary tables show important information about the identified proteins, including their UniProtKB accession numbers, IDs, protein and gene names, adjustedp values, and associated functions. The proteins that were significantly differentially expressed according to age included 38 downregulated proteins and 57 upregulated proteins (Supplementary Table 1 and Supplementary Table 2). In an interaction analysis of the 95 DEPs, we counted the proteins with more than 10 nodes of protein interactions (Figure 3(a)). We compared our DEPs with 307 human aging-related genes from the GenAge database and identified 6 genes that intersected with the DEPs, namely, CDC42, CLU, A2M, LMNA, VCP, and AIFM1 (Figure 4(c)). In a proteomics study on the senescence-associated secretory phenotype [30], serine protease inhibitors (SERPINs) were shown to be significantly associated with age changes in human cohorts and to be potential biomarkers. In the induced aging model in that study, the secretion levels of the proteins SERPINE1, SERPINF1, SERPINF2, SERPING1, and SERPINH1 were significantly increased. In contrast, SERPINA3, SERPINB1, SERPINB6, and SERPINC1 were significantly upregulated with age in our study. In an age-associated proteomics study in a healthy population [29], 217 age-associated proteins were present, 9 of which intersect with the DEPs in our study: ANXA2, S100A4, B2M, FGA, FGB, FGG, SERPINA3, CST3, and HP (Figure 4). In the previous study, SERPINA3 and SERPING1 were positively correlated with age, while SERPINF2 was negatively correlated with age, which is consistent with our results.Figure 3
(a) Proteins possessing more than 10 neighboring proteins in the protein-protein interaction network analysis. (b–e) Network analysis of gene expression was performed for the 21 selected genes, where the genes we selected were indicated by red circles. The pathway with the most genes involved in (b) is complement and coagulation cascades; the pathway with the most genes involved in (c) is cardiomyopathy; the pathway with the most genes involved in (d) is ferroptosis; and the pathway with the most genes involved in (e) is apoptosis and autophagy. (Supplementary Tables 5, 6, 7, and 8 for details).
(a)(b)(c)(d)(e)Figure 4
(a) Symmetric scatter plots compare the differential status of proteins in the age group. The horizontal and vertical coordinates represent the mean relative expression of young and old groups; the dashed line in the graph represents the threshold line at∣log2FC∣=1.5. The color of the scatter represents the p value. (b) Wayne plots show the distribution of protein expression between the young and old groups, with coverage areas representing differentially expressed proteins. (c) Wayne diagram of our proteomic results and 307 genes from the GenAge Human Aging Gene Database. (d) Venn diagram of our proteomic results and an aging-related plasma proteomics study in a healthy population.
(a)(b)(c)(d)


## 3.3. Functional Analysis of the DEPs

Through DIA mass spectrometry, we identified a total of 1,270 unique proteins, of which 95 were differentially expressed in the different age groups. The 95 identified proteins were analyzed using the GO and KEGG databases (the functional classifications of the 95 DEPs are shown in Supplementary Table3). The functional pathways into which the DEPs were classified are shown in Figure 5(a) including the cellular process, environmental information processing, genetic information processing, human disease, metabolism, and organismal system pathways. We analyzed the DEPs with the KEGG and found that the most enriched pathways included platelet activation, complement and coagulation cascades, and the VEGF signaling pathway. Among the DEPs, BANF1 was found to be responsible for the chromatin structure and dynamics, replication, recombination, and repair. Energy production and transfer are also entry points for the aging process, and LDHB and HBA1 were enriched in these pathways. The signaling pathways that we focused on in relation to aging included ferroptosis and apoptosis, both of which are cell growth- and death-related metabolic pathways. Among the six proteins intersecting with GenAge data, AIFM1 and LMNA are involved in apoptosis. In the gene-level network analysis of the 21 most promising proteins we screened (Figures 3(b)–3(e)), the significantly enriched pathways were the complement and coagulation cascades, platelet activation, ferroptosis, apoptosis, autophagy, and PPAR signaling pathways. Among them, the PPAR signaling pathway may be closely related to the aging process [31]. (For more detailed gene node information, please refer to Supplementary Tables 4, 5, 6, 7, and 8 and the following website: https://www.networkanalyst.ca/NetworkAnalyst/faces/Share?ID=_9qljssi8x.) A protein of particular interest was ceruloplasmin (CP), which had 16 neighboring proteins in the protein-protein interaction network. CP is closely associated with aging [32, 33], and the biological metabolic pathway in which it participates is the well-known ferroptosis pathway (Figure 6), a signaling pathway that is closely associated with aging. Our age-associated proteomics results are consistent with the results obtained in other age-associated proteomics studies with independent cohorts (Supplementary Tables 1 and 2 at 10.1038/s41591-019-0673-2).Figure 5
(a) Pathway classification of differentially expressed proteins, thex-axis represents the number of differentially expressed proteins, the y-axis represents the pathway-annotated entries. (b) Pathway relationships; red and blue dots indicate up- and downregulated differential proteins, respectively. Purple dots indicate the top ten enriched pathways, darker dots indicate significant enrichment, lighter dots indicate enrichment but not significant, and larger areas indicate higher enrichment. The different color lines indicate the different pathway categories: blue for environmental information processing, purple for human disease, orange for metabolism, and yellow for organismal systems.
(a)(b)Figure 6
Changes in proteins in the Ferroptosis signaling pathway with aging. CP is a protein that is upregulated with aging and also highlights important gender differences with aging, and it plays an important role in this well-known cellular aging and death pathway.


## 3.4. Potential Candidate Biomarkers of Sex Differences in the Aging Process

We compared the proteins that changed most significantly under the influence of two factors, age and sex, according to a dual statistical threshold defined by thep value and log2FC. The proteins that were most significantly upregulated in the elderly group were LGALS3BP, HBA1, and HBB, while the proteins that were most significantly upregulated in the young group were ABRAXAS2, BANF1, VAT1, and KRT35 (Figure 4(a)). The proteins in the graph are colored according to the significance of the p value, and the graph clearly shows that a smaller p value was associated with a greater difference in protein expression between the groups. Women generally live longer than men [18], and aging rates differ between men and women [34]. From the 95 proteins identified in our quantitative proteomics experiments and through statistical analysis via independent sample paired t-tests between the age and sex groups, we selected proteins as candidate biomarkers for distinguishing differences in the aging process between men and women (Figure 7). The ROC curves of several proteins of interest are shown in Figure 8, and the corresponding area under the curve (AUC) values and p values were calculated. ABRAXAS2, VAT1, and LGALS3BP were identified as promising candidate biomarkers of aging. t-tests for these three proteins revealed a p value < 0.05 between men and women in the younger group and a p value > 0.05 between men and women in the older group (Table 3). Thus, the three proteins were expressed differently between sex groups in the young population, but this difference disappeared with aging. Among these proteins, ABRAXAS2 and VAT1 were downregulated in the elderly group and LGALS3BP was upregulated in the elderly group. The results of tandem mass spectrometry of the three proteins are shown in Figures 9(a)–9(c).Figure 7
(a–c) Violin plots of the relative expression abundance of three proteins, ABRAXAS2, LGALS3BP, and VAT1, in the young-old group. The red solid line is the mean value and the blue dashed line is the quartile. The data in the vertical coordinates are protein expression abundance values taken aslog2. (d) Heatmap of six proteins of a particular value. (e) Mean values of the expression abundance of the three proteins in young and old males and females. The p values between the groups are marked. It can be seen that the three proteins differ between young men and women and do not differ between older men and women.
(a)(b)(c)(d)(e)Figure 8
(a–f) Receiver operating characteristic curve (ROC curve) with an AUC andp value obtained by comparing six proteins in young and older people.
(a)(b)(c)(d)(e)(f)Table 3
Different protein expressions at different ages in different genders. Thet-test showed significant differences among age groups (∗p<0.05) and between genders for ABRAXAS2, VAT1, and LGALS3BP (∗∗p<0.05). More details are described in Section 3. The case group is an older group (55–70 years) and the control group is a young group (20–30 years). M: mean protein expression of male; F: mean protein expression of female.
AccessionIDDescriptionGeneOld and young differencesOld group differences in gendersYoung group differences in gendersDownQ15018ABRX2_HUMANBRISC complex subunit Abraxas 2ABRAXAS2∗p=0.001p=0.302M=13.45761365F=13.10540960∗∗p=0.016M=13.6814853567F=14.2325095143Q99536VAT1_HUMANSynaptic vesicle membrane protein VAT-1 homologVAT1∗p=0.004p=0.729M=15.87932308F=15.78413334∗∗p=0.038M=17.2356184533F=16.3622206371UpQ08380LG3BP_HUMANGalectin-3-binding proteinLGALS3BP∗p=0.001p=0.549M=15.20773285F=15.01199494∗∗p=0.034M=14.8493189733F=14.2752645800Figure 9
(a) The sequence NH2-AIYQVYNALQEK-COOH allows the identification of ABRAXAS2. (b) The sequence NH2-AVDTWSWGER-COOH allows the identification of LGALS3BP. (c) The sequence NH2-TEAASDPQHPAASEGAAAAAASPPLLR-COOH allows the identification of VAT1.
(a)(b)(c)


## 4. Discussion

The process of skin aging is very complicated. In addition, due to differences in hormone levels and the immune system, there are also differences in the aging process between men and women [18, 35, 36]. The expression of epidermal proteins can be decrypted using proteomics techniques with DIA. This study is aimed at identifying and quantifying proteins that were differentially expressed at different ages while looking for key information regarding the differences between different sexes during aging. Reliable statistical analysis techniques revealed that 95 proteins were significantly differentially expressed according to age; of these, 57 were upregulated with aging and 38 were downregulated with aging.The results of this study are consistent with those of previous studies, which confirms the reliability of our proteomics experimental data. In our study, the results for SERPIN protein family members, FGA, FGB, FGG, and other proteins were consistent with other aging-related proteomics findings. Xu et al. found that low expression of SERPINC1 may suppress the PI3K/AKT/mTOR signaling pathway [37] and the suppression of the PI3K/AKT/mTOR signaling pathway can delay aging [38]. We also obtained some results that differ from previous findings, such as the results regarding the expression of GPC1 protein; in previous studies, GPC1 has been found to be significantly downregulated with aging [39], while it was upregulated with aging in our study. This difference between studies may have been related to sample size differences, laboratory errors, etc. The results of pathway enrichment analysis showed that the most enriched pathway was platelet activation [40]. There is evidence that oxidative stress contributes to the high reactivity of platelet activation during aging. FGA, FGB, FGG, VASP, and SNAP23 detected in the skin are all related to this pathway. In addition, another pathway, glycerophospholipid metabolism, was noteworthy, as metabolic dysfunction is a common feature of the aging mechanism [41]. The VEGF signaling pathway was also of interest. As age increases, the blood flow in the dermal blood vessels decreases, which leads to decreases in the proliferation of skin fibroblasts that may be related to impairment of the VEGF signaling pathway [42]. Finally, energy production and conversion are also related to the aging of cells [43]. One of the proteins involved in these processes is GPD2, which helps block the IMPM2L-GPD2 and IMMP2L-AIF signaling axes in order to drive the aging process instead of cell death under oxidative stress [44]. One study has shown that activation of respiratory chain-related GPD2 increases reactive oxygen species (ROS) release from mitochondria [45]. Although recent studies have demonstrated adaptive responses to oxidative stress, ROS remain indirect contributors to oxidative aging [46, 47]. This situation is also consistent with our finding that GPD2 expression is upregulated with increasing age. The very specific metabolic patterns summarized above are helpful for understanding the potential mechanisms leading to aging.Three of the 95 proteins that were differentially expressed in the different age groups were differentially expressed between men and women, but only among young people. This difference disappeared with aging. Studies have shown that aging does not proceed at a constant rate; there are inflection points during youth, middle age, and twilight age. Due to inconsistent proteomics results, Wyss-Coray et al. hypothesized that the aging rates of men and women also differ [19]. There are many pathways related to aging; the most studied of which are the insulin-like signaling pathway, the mTOR signaling pathway, the p53 signaling pathway, ferroptosis, and autophagy. Preliminary studies have shown that autophagy suppresses chromosomal instability by maintaining cellular metabolic homeostasis [48] and genomic instability is considered to be one of the nine characteristics of aging [49]. The occurrence of iron ion-dependent lipid peroxidation in cells is a marker of ferroptosis, which requires accumulation of iron ions, and blood iron metabolism influences healthy longevity [50]. Time-dependent excessive deposition of iron ions leads to genomic instability and DNA damage [51], and DNA damage can further promote aging. At the same time, excess iron ions inhibit expression of the protein p53, which can repair DNA damage, eventually initiating ferroptosis and accelerating aging [52]. Accumulation of excess iron ions in aging cells results in ferroptosis resistance, autophagy dysfunction, and evasion of the self-renewal mechanism and accelerates organismal aging [53]. In summary, ferroptosis contributes to aging in multiple ways by affecting the p53 signaling pathway, increasing genomic instability, affecting autophagy, and increasing ROS accumulation. In our study, CP, a primary copper-binding protein and an important antioxidant [54] that prevents the production of free radicals [55], was identified as a DEP. CP has been shown to be associated with oxidative modifications in age-related changes [56] and with brain aging [32]. It is also involved in ferroptosis, an important signaling pathway associated with aging.In our study, on the basis oft tests, only 3 of the 95 DEPs between age groups had statistically significant sex differences and the differential expression of the 3 proteins between sexes existed only in young people. This finding may be related to the difference in the aging rate between men and women. Among the DEPs, ABRAXAS2, also known as ABRO1, is a scaffold protein that recruits a variety of peptides to assemble BRCC36 isopeptide complex (BRISC) deubiquitinating enzyme (DUB) complexes. It has been shown that inhibitors of DUBs can induce autophagy and thereby inhibit cell death [57]. The BRISC enzyme has deubiquitination activity and consists of four known subunits: MERIT40; BRE, which is expressed in the brain and reproductive organs; BRCA1/BRCA2 complex, subunit 3 (BRCC3); and ABRO1 [58]. BRCC3, which was also among the 95 DEPs that we identified, exhibited a relationship with ABRO1 and PSMB7 in the protein-protein interaction network. Excessive oxidative stress leads to aging [59], especially in the skin, which is the first line of defense against various types of external stimuli and is subject to oxidative stress. ABRO1, a newly discovered p53-dependent regulator, plays important roles in tumor suppression and the DNA damage response [60]; consistent with these roles, studies have pointed out that ABRO1 is related to telomere replication and genome repair [61]. A blocked chromosome replication fork is essential for the maintenance of chromosome stability. ABRO1 protects the integrity of the blocked replication fork and thus maintains the stability of the genome [62]. The expression of ABRO1 was downregulated with age in our study, and cells’ ability to respond to oxidative stress decreases with age, which may imply that ABRO1 is an important indicator of aging. Second, VAT1 was expressed at high levels in the skin, with the overall expression being higher in males than females, which is consistent with results from the Human Protein Atlas (HPA) database (https://www.proteinatlas.org/ENSG00000108828-VAT1/tissue). VAT1 participates in calcium-regulated keratinocyte activation in the epidermal repair mechanism and has ATPase activity. During the EuKaryotic Orthologous Groups (KOG) analysis of our data, VAT1 was found to be enriched for the processes of energy production and conversion. Studies have shown that hVAT-1 plays a specific role in the physiology of keratinocytes, especially in the process of calcium regulation, and the gene encoding VAT1 is located on the BRCA1 locus [63, 64]. It may also form the BRISC enzyme together with the aforementioned ABRO1 protein, which is involved in deubiquitination activities, while ubiquitination is involved in DNA repair and stabilization of protein complexes [65]. Proteasome activity is reduced during aging, in which ubiquitination is the target event. Some of the most important biological processes associated with skin aging are ubiquitin-induced proteolysis and cellular metabolism [66]. Therefore, we speculate that the downregulation of VAT1 and ABRO1 with age is an important marker of ubiquitination and even aging. Third, LGALS3BP, also known as Mac-2 binding protein (M2BP), is involved mainly in biological processes that promote integrin-mediated cell adhesion and platelet degranulation and may stimulate the host to defend against viruses and tumor cells. M2BP can combine with collagens IV, V, and VI and with fibronectin [67], and M2BP is significantly related to reactive oxygen metabolites [68]. Increases in ROS, platelet activation, adiponectin, and oxidative stress are important factors leading to aging. The functional relationship between M2BP and aging needs further study. To date, no publications have reported the proteins ABRAXAS2, VAT1, and LGALS3BP as candidate biomarkers for aging. In addition, the differential expression of these proteins may suggest a difference in the rate of aging between men and women, but there are no data regarding the mechanisms by which imbalances in the three proteins play roles in skin aging.In our experiments, we identified 95 age-differentiated proteins in the forearm skin of 20 healthy Chinese volunteers and 3 of them were key proteins that differed between age groups and sexes. Comparison of our results with the results of other studies designed to identify skin aging biomarkers revealed some differences that may be attributable to differences in regional categories, genetic factors, sex makeup, and sample processing during the experiment. Results that were highly consistent with the results of other studies were also obtained, which indicates the validity of the method and results of the current experiment. However, in the present experiment, more correlations than causal relationships were found, so further experimental confirmation is needed. In quantitative proteomics studies, posttranslational modifications of proteins are also worthy of attention and both protein glycosylation and protein carbamylation are hallmarks of aging [69, 70]. Multiomics approaches will provide richer perspectives on aging than single proteomics techniques [71]. Multiomics analyses have been widely used for research on cancer [72] but have received very little application in the field of aging. Network-based integrated multiomics analyses may provide new insights into the nonlinear processes of aging and aid in the search for markers of inflection points.


## 5. Conclusion

We conducted an in-depth exploration of the epidermis based on a proteomics technology platform. We have demonstrated that our method is able to detect senescence- and sex-related epidermal proteins with significant differences, and we also found some candidate biomarkers that exhibit age-dependent sex differences. These results demonstrate the usefulness of our new proteomics assay for further elucidating biological processes, pathways, and biomarkers of skin aging. In conclusion, our findings indicate that the combination of noninvasive extraction of epidermal proteins and DIA mass spectrometry has great potential for future research in the field of skin health and disease.



---

*Source: 8849328-2021-04-09.xml*