Ancient proteins identify various Denisovan remains from Southwest China

Archaeological site and samples

Bianfu Cave (26° 27′ 43.65′′ N, 100° 09′ 55.45′′ E) is in Heqing County, western Yunnan Province in southwestern China, approximately 2 km south of the Yinhe River, a secondary tributary of the Yangtze River (Fig. 1a). In 2019, the Yunnan Institute of Cultural Relics and Archaeology and other cooperative institutes conducted excavations on this site. The stratigraphy comprises Holocene deposit, artefact-bearing cultural deposit and underlying cave-fissure deposit. The cultural deposit is further divided into upper (layers 3–6) and lower (layers 7–12) cultural layers (Fig. 1b). A Bayesian age model combining optically stimulated luminescence and U-series dates yielded a time range spanning approximately 190 ka to 70 ka¹, making it the earliest Palaeolithic site discovered in western Yunnan. All cultural layers yielded lithic artefacts, totalling 1,366 stone artefacts, including both core and flake tools¹. More than 64,000 mammal fossils were recovered from the cultural layers, among which more than 60,000 were fragmentary bones. On the basis of the limited number of morphologically identifiable mammal fossils, more than 100 individuals of deer and bovids were discovered, together with small numbers of rhinoceros, hyenas and other animals¹.

Four hominin teeth were unearthed in the middle of layer 7 (between 148–139 ka and 142–134 ka)¹, suggesting a high likelihood of finding more hominin fossils in this site. The four hominin teeth included a lower left lateral incisor (I₂), a lower left fourth premolar (P₄) and two lower right second molars (M₂). Morphological analysis showed that the lower M₂ teeth exhibited close morphological affinity with those of Xiahe and Laos Denisovans, as well as the East Asian late Middle Pleistocene archaic Homo from Hualongdong¹. Two of these teeth, a lower fourth premolar (YHB3518) and a lower second molar (YHB3075) (Fig. 1c,d), were selected for proteomic analysis to further confirm their taxonomic attribution. An acid etch approach was used for analysis of the enamel, whereas dentine powder was drilled for protein extraction, resulting in one enamel fraction and three dentine fractions for liquid chromatography coupled with tandem mass spectrometry (LC–MS/MS) analysis.

After manual inspection and traditional morphological measurement of more than 60,000 bone fragments on the basis of bone size, cortical thickness, articular surface morphology, muscle attachment marks on long bone surfaces, vascular impressions on cranial internal surfaces and so on, 22 fragments were considered as potential hominin remains, five from the upper layers and 17 from the lower layers. Before conducting protein extraction on these 22 remains, we randomly selected 36 bone fragments from the upper layers and 18 from the lower layers for preliminary ZooMS analysis to assess protein preservation of the site. Thus, a total of 76 bone fragments were analysed by ZooMS in 6 batches, with each batch including an extraction blank and a positive control (either modern camel bone or dentine powder) (Supplementary Data 1). After ZooMS identification, 3 hominin remains and 14 mammalian samples from different layers were submitted to LC–MS/MS analysis, including the corresponding extraction blanks. To obtain a more comprehensive proteome of the hominin remains at this site, we performed further proteomic extraction on the parietal BFD767, yielding another three fractions for LC–MS/MS analysis. We also performed U-series dating to estimate the minimum age of three hominin bone samples.

U-series dating

U-series dating was conducted on three hominin bone samples and their attached carbonates. For each sample, a series of 5–11 spot pairs were analysed using a laser ablation system coupled to a multicollector inductively coupled plasma mass spectrometer (LA–MC-ICPMS). The LA system (RESOlution-LR, Applied Spectra) was equipped with a Coherent COMPex Pro102 ArF Excimer laser source and an S155 large-format sample pool. The Neptune double-focusing MC-ICPMS (Thermo Fisher) had nine Faraday cups and a secondary electron multiplier, allowing simultaneous measurement of several isotopes over a relative mass range of 17%. The Neptune interface was upgraded with a Jet-sample cone, an X-skimmer cone and a high-efficiency dry pump to improve the ion transport efficiency. By tuning the analytical parameters (carrier gas flow, torch position and zoom optics) of our LA–MC-ICPMS, we could routinely obtain ²³⁸U signals of approximately 0.8–1.0 V and ²³²Th/²³⁸U ratios of approximately 90–95% on the international standard NIST 612, with a spot size of 173 μm, a pulse repetition of 10 Hz, an energy fluence of 5 J cm⁻² and a scan speed of 3 μm s⁻¹. We conducted a preablation step to remove impurities from the sample surface. Then, each spot was ablated for 120 s at 10 Hz pulse repetition and 5 J cm⁻² fluence. With signal intensities of ²³⁸U and ²³⁰Th, the spot size ranged from 50 μm to 180 μm. These laser ablations produced a crater approximately 150 μm deep. We measured isotopes ²³⁰Th and ²³⁴U using a secondary electron multiplier and simultaneously measured ²³²Th, ²³⁵U and ²³⁸U in Faraday cups^47,48,49,50. For each spot, around 150 cycles were performed to acquire isotopic data, with the first 40–50 cycles collected before sample ablation and used as the gas blank. An in-house fossil bone standard (RM-B1) and a carbonate standard (RM-C1) were used as the matrix-matched standards for calibration of laser-induced elemental and isotopic fractionation for the hominin bones and carbonates, respectively.

ZooMS screening

Around 100–200 mg of the sample was used for protein extraction with the conventional acid-insoluble approach⁵¹. After demineralization in 0.6 M HCl, the pellet was incubated with 50 mM ammonium bicarbonate (ABC) at 65 °C for 3 h. The supernatant was then collected and digested with trypsin. The digested peptide mixture was desalted using Pierce C18 Pipette Tips (Thermo Scientific), and then subjected to matrix-assisted laser desorption/ionization time-of-flight mass spectrometry (MALDI-TOF MS) analysis.

Palaeoproteomic extraction

For the identified hominin parietal BFD767, a second round of protein extraction was conducted to obtain more proteins. The protocol was modified from that reported in a previous publication⁵². Another subsample (around 300 mg) was demineralized in 0.6 M HCl, and the acid supernatant was collected and ultrafiltered through an Amicon Ultra-4 (3 kDa) centrifugal filter unit. After the filtration unit had been washed with buffer (50 mM ABC), we obtained the ‘acid-soluble fraction 3KD’ by dissolving the proteins retained on the filter in 50 mM ABC. The pellet was further incubated in 50 mM ABC at 65 °C for 3 h. The supernatant was reduced with Tris(2-carboxyethyl) phosphine HCl (final concentration 0.1 M) at 56 °C for 30 min and alkylated using iodoacetamide (final concentration 0.1 M) at room temperature in the dark for 30 min. The supernatant was divided into two aliquots, one for trypsin digestion and the other for elastase digestion. The trypsin aliquot was ultrafiltered through an Amicon Ultra-4 (3 kDa) centrifugal filter unit, washed with buffer (ABC) and dissolved in 50 mM ABC; this was referred to as the ‘acid-insoluble fraction trypsin’. The elastase fraction was ultrafiltered through an Amicon Ultra-4 (3 kDa) centrifugal filter unit, washed with buffer (Tris-HCl solution) and dissolved in 50 mM Tris-HCl; this fraction was referred to as the ‘acid-insoluble fraction elastase’.

The ‘acid-soluble fraction 3KD’ and ‘acid-insoluble fraction trypsin’ were digested with porcine trypsin (Promega, 2 μl, 0.5 μg μl⁻¹) overnight at 37 °C, after which trifluoroacetic acid (TFA) was added to the digestion (final concentration 0.1%) to stop the reaction. The ‘acid-soluble fraction elastase’ was digested with elastase (Promega, 6 μl, 0.2 μg μl⁻¹) overnight at 37 °C; then, TFA was added to the digestion (final concentration 0.5%) to stop the reaction. After a desalting step using Pierce C18 Pipette Tips, these peptide mixtures were dried for LC–MS/MS analysis.

For two hominin teeth, YHB3518 and YHB3075, enamel proteins were extracted using an acid etch method modified from that described in ref. ⁵³. After the initial 2-min etch using 5% (v/v) HCl, the solution was discarded. Then, two rounds of 15-min etch were carried out, and the solutions were combined to form the ‘acid etch fraction’. No reduction, alkylation or digestion was performed on this fraction. For dentine protein extraction, around 30 mg of powder (35 mg for YHB3518; 32 mg for YHB3075) was demineralized in 0.6 M HCl. The supernatant was collected and ultrafiltered through an Amicon Ultra-4 (3 kDa) centrifugal filter unit. After washing, the proteins were dissolved in 50 mM ABC and formed the ‘acid-soluble fraction 3KD’. The acid-insoluble pellet was washed and incubated with 50 mM ABC at 65°C for 3 h. The supernatant was collected as the ‘acid-insoluble fraction trypsin 1’. The remaining pellet was washed and incubated again. The supernatant with pellet was denoted the ‘acid-insoluble fraction trypsin 2’. No reduction or alkylation was conducted, and all three fractions were digested by trypsin. After desalting with Pierce C18 Pipette Tips, these peptide mixtures were analysed by MALDI-TOF MS for assessment of the collagen signal. As the ‘acid-soluble fraction 3KD’ of YHB3075 did not show any obvious collagen signal, it was combined with the ‘acid-insoluble fraction trypsin 1’ for further LC–MS/MS analysis. The other fractions were dried separately for LC–MS/MS analysis.

A blank was extracted along with the samples for estimation of modern contamination during the experiment. All the sample treatments were performed in the dedicated clean room for ancient protein analysis in the Molecular Paleontology Laboratory, Institute of Vertebrate Paleontology and Paleoanthropology (IVPP), Chinese Academy of Sciences (CAS).

MALDI-TOF MS

One microlitre of the peptide mixture was spotted on to an MTP384 Bruker ground-steel MALDI target plate. Then, 1 µl of α-cyano-4-hydroxycinnamic acid matrix solution (1% in 50% acetonitrile/0.1% TFA (v/v/v)) was added on top and mixed with the peptide mixture. Each sample was analysed on a Bruker autoflex maX MALDI-TOF mass spectrometer in triplicate. The acquired mass ranges were slightly different across several batches and included m/z 600 to 3,500, m/z 600 to 4,000 and m/z 700 to 3,500. Text files were converted from the raw files and processed using mMass v.5.5.0 (ref. ⁵⁴).

Liquid chromatography coupled with tandem mass spectrometry

For protein extraction from the three identified hominin bone samples (including the initial ZooMS extraction and the further protein extraction on the parietal BFD767), we performed LC–MS/MS analysis using two instruments from Central Laboratory at Capital Medical University: an Orbitrap Fusion Lumos mass spectrometer and an Orbitrap Exploris 480 mass spectrometer. Both machines were interfaced with an EASY-nLC 1200 system (Thermo Fisher Scientific). The peptides were loaded on to a trap column (100 μm i.d. × 2 cm, C18), followed by separation on an analytical column (150 μm i.d. × 25 cm, C18). Mobile phase A was 0.1% formic acid in water, and mobile phase B consisted of 80% acetonitrile and 0.1% formic acid. For Orbitrap Fusion Lumos, the peptides were eluted at a flow rate of 500 nl min⁻¹ using a 120-min linear gradient program: 0–8 min, 7–11% B; 8–96 min, 11–28% B; 96–108 min, 28–40% B; 108–113 min, 40–90% B; 113–120 min, 90% B. Full MS data were acquired across an m/z range of 375–1,400 with a resolution of 120,000, an AGC target of 250% and a maximum injection time of 50 ms. The MS/MS scan was conducted at a resolution of 15,000 with an AGC target of 100%, a maximum injection time of 22 ms and a normalized collision energy of 35%. For the Orbitrap Exploris 480, peptides were eluted at a flow rate of 350 nl min⁻¹ using a 125-min linear gradient program: 0–1 min, 1–9% B; 1–101 min, 9–28% B; 101–113 min, 28–42% B; 113–117 min, 42–99% B; 117–125 min, 99% B. Full MS data were acquired across an m/z range of 350–1,500 with a resolution of 60,000, an AGC target of 300% and a maximum injection time of 50 ms. The MS/MS scan was collected at a resolution of 15,000 with an AGC target of 75%, a maximum injection time of 22 ms and a normalized collision energy of 30%. Extraction blanks from different sample batches were analysed alongside the samples. Injection blanks were included before and after each sample analysis to minimize possible carryover between runs.

For the protein extraction from two hominin teeth, we performed LC–MS/MS analyses using two instruments: an Orbitrap Fusion Lumos mass spectrometer (Thermo Fisher Scientific) from Central Laboratory at Capital Medical University; and an Orbitrap Exploris 480 mass spectrometer (Thermo Fisher Scientific) from State Key Laboratory of Genetics and Development of Complex Phenotypes at Fudan University. Both devices were coupled to an EASY-nLC 1200 HPLC system (Thermo Fisher Scientific). For the Orbitrap Fusion Lumos at Capital Medical University, the gradient program and mass spectrometry parameters were the same as those used for the hominin bone sample. For the Orbitrap Exploris 480 at Fudan University, the peptides were loaded on a 75 μm i.d. × 25 cm analytical column, which was packed in-house using reversed-phase silica of 1.9 μm (Reprosil-Pur C18 AQ, Dr. Maisch GmbH). Mobile phases A and B were as described previously. The peptides were separated using an 80-min gradient: 5–8% B, 2 min, at a flow rate of 200 nl min⁻¹; 8–44% B, 38 min, 200 nl min⁻¹; 44–70% B, 8 min, 200 nl min⁻¹; 70–100% B, 2 min, 200 nl min⁻¹; 100% B, 10 min, 200 nl min⁻¹; 100–5% B, 2 min, 200 nl min⁻¹; 5% B, 2 min, 300 nl min⁻¹; 5–100% B, 6 min, 300 nl min⁻¹; 100% B, 10 min, 300 nl min⁻¹. Full MS scans were collected during the first 65 min, after which the column was washed and re-equilibrated for 15 min without data collection. The full MS data acquisition was conducted across an m/z range of 350–1,600 with a resolution of 60,000, an AGC target of ‘Standard’ and the maximum injection time mode set to ‘Auto’. The MS/MS spectra were acquired with a resolution of 15,000, an AGC target of ‘Standard’, a maximum injection time of 30 ms and a normalized collision energy of 30%.

To provide a comparative reference for the extent of fossil degradation across different layers of Bianfu Cave, we also selected 14 animal samples for further LC–MS/MS analysis. These samples were processed on an Orbitrap Fusion Lumos mass spectrometer from State Key Laboratory of Genetics and Development of Complex Phenotypes at Fudan University, coupled to an EASY-nLC 1200 HPLC system (Thermo Scientific). Mobile phases A and B were as described previously. The peptides were separated on an analytical column (75 μm i.d. × 20 cm, C18) using an 80 min gradient: 2–5% B, 3 min, at a flow rate of 200 nl min⁻¹; 5–35% B, 40 min, 200 nl min⁻¹; 35–44% B, 5 min, 200 nl min⁻¹; 44–100% B, 2 min, 200 nl min⁻¹; 100% B, 10 min, 200 nl min⁻¹; 100–5% B, 2 min, 200 nl min⁻¹; 5% B, 2 min, 300 nl min⁻¹; 5–100% B, 6 min, 300 nl min⁻¹; 100% B, 10 min, 300 nl min⁻¹. The full MS data acquisition was conducted over the first 65 min, with no data collected over the subsequent 15 min during column wash and re-equilibration. Full scans ranging from m/z 350 to m/z 1,600 were measured, with a resolution of 60,000, an AGC target of 100% and a maximum injection time of 50 ms. The MS/MS spectra were acquired at a resolution of 15,000 with an AGC target of 100%, a maximum injection time of 30 ms and a normalized collision energy of 30%.

Proteomic data analysis

To confirm the taxonomic identification of these bone fragments, we conducted two analyses on the MS/MS data file from each LC–MS/MS run (not including the second round of protein extraction for the hominin parietal BFD767). First, the raw data file was searched against a custom mammal type I collagen database through PEAKS Online v.12 or v.11.5 (ref. ⁵⁵). This database comprised two predominant bone proteins (COL1A1 and COL1A2), the sequences of which were collected from GenBank, UniProt and previous publications^{56,57,58,59,60,61,62,63,64,65,66}. The translated protein sequences from high-coverage genomes of Denisova 3 (ref. ³¹), two Neanderthals (Altai Neanderthal, Vindija19)^34,39 and Ust’-Ishim⁴⁰ were also included. PEAKS searches, including peptide de novo, PEAKS DB, PTM and SPIDER, were performed with the following parameters: parent ion mass tolerance of 10 ppm, fragment ion mass tolerance of 0.05 Da, enzyme specificity set to SemiTryptic, maximal missed cleavages of 2, maximal modification number of 6 for each peptide, and variable modifications including deamidation (NQ), oxidation (M), hydroxylation (P), pyro-Glu from E and pyro-Glu from Q. The peptide length range was set to 8–45. The false discovery rate (FDR) was set to 1% (peptide level), and protein scores were filtered with −10log₁₀P ≥ 20 and average local confidence (%) ≥ 50 (de novo only). A common contaminant database was also included in each database search, as in a previous study⁵, and we analysed the corresponding extraction blanks using the same strategy as for the samples to monitor possible modern contaminants during the experiments. Peptides detected in the negative control sample or assigned to contaminant proteins were excluded. The remaining peptides from each sample were used to calculate the extent of glutamine and asparagine deamidation using a modified Python script based on peptide intensity⁶⁷.

To validate the taxonomic identification results, we used another set of pipelines based on Mascot search and the ClassiCOL pipeline²⁵. Raw data files from the bone samples were searched against a published collagen database containing 45 collagen homologues from 221 mammalian species²⁵ using Mascot v.2.6 (Matrix Science). The parameters were set as follows: parent ion mass tolerance of 10 ppm, fragment ion mass tolerance of 0.05 Da, enzyme specificity SemiTryptic, maximal missed cleavages of 2, and variable modifications including deamidation (NQ), oxidation (M) and hydroxylation (P). After being filtered with a significance threshold of 0.01, the results were submitted to the ClassiCOL pipeline with the taxonomy restricted to Mammalia (python ClassiCOL.py -d {path_to_the_script} -l {path_to_search_results} -s MASCOT -t Mammalia -c 20).

Following taxonomic identification, we performed a database search against the Hominidae bone database using raw data files from all three identified hominin bone samples to retrieve their proteomes. This database comprised the entire human proteome (from UniProt, one protein sequence per gene; 20,575 entries), supplemented with common bone and dentine protein sequences translated from 4 ancient hominin individuals (Denisova 3, Altai Neanderthal, Vindija19 Neanderthal and Ust’-Ishim)^31,34,39,40. Published protein sequences for several other archaic individuals^2,3,4,5,41 were also included in the database. The main search was conducted in PEAKS Online⁵⁵ using the same parameters as for the search against the mammal type I collagen database with the following minor modifications. For the raw data files from the second round of protein extraction for the hominin parietal BFD767, carbamidomethylation (C) was added as a variable modification in the ‘acid-insoluble fraction trypsin’ and ‘acid-insoluble fraction elastase’, and enzyme specificity was set to ‘unspecific’ for the ‘acid-insoluble fraction elastase’. The possible mutations identified in the SPIDER search were added to the original database for iterative data search and verification. Peptides matching multiple genes were considered to be non-unique and excluded from subsequent analysis. Given the high similarity between different collagen chains, the permutation and Unipept confirmation filtering approach was used to validate the identification of low-abundance collagens (the pept2prot function from Unipept was used; Supplementary Text 2 and Supplementary Fig. 18). After this filtering and validation, the extent of glutamine and asparagine deamidation for each identified protein was calculated from peptide intensities using a Python script modified from a previous study⁶⁷. For a reliable calculation, each protein was required to have at least two Q/N-containing peptides with valid intensities. Endogenous proteins were confirmed on the basis of elevated deamidation values.

For the raw data files from three dentine fractions of two hominin teeth, a database search against the Hominidae bone database was carried out using the same parameters as those used for the search against the mammal type I collagen database for bone samples. The data processing strategy was the same as that used for the hominin bone samples. For the raw data files from the acid etch fraction, an initial search was conducted against the entire human proteome. The parameters were set as follows: precursor ion mass tolerance of 10 ppm; fragment ion mass tolerance of 0.05 Da; enzyme specificity of ‘unspecific’; variable modifications including deamidation (NQ), oxidation (M), hydroxylation (P), phosphorylation (STY), pyro-Glu from E and pyro-Glu from Q; maximum modifications per peptide of 3; and peptide length range of 6–45. FDR was set to 1% at the spectral level, and protein scores were filtered with −10log₁₀P ≥ 20 and average local confidence (%) ≥ 50 (de novo only). In addition to the main enamel proteins, dentine collagens were identified in YHB3518 with relatively abundant peptides and elevated deamidation rates, probably because the acid etch process accessed the dentine region of this tooth and extracted related proteins. Two serine protease inhibitors, SERPINC1 and SERPINA1, were identified in YHB3075; these proteins exhibited elevated deamidation rates. This was not unexpected, as these proteins have also been reported in archaic and modern enamel^32,68,69. Therefore, for each tooth, a specific Hominidae enamel database was built for a second round of search. In addition to the 12 commonly used proteins (AHSG, ALB, AMBN, AMELX, AMELY, AMTN, COL17A1, ENAM, KLK4, MMP20, ODAM and TUFT1), three main dentine collagens (COL1A1, COL1A2 and COL2A1) were included in the database for YHB3518, whereas two serine protease inhibitors (SERPINC1 and SERPINA1) were included in the database for YHB3075. Only species belonging to Hominidae were included; these comprised archaic hominins such as Denisovans and Neanderthals. The protein sequences, including both canonical and isoform sequences, were obtained from UniProt and the translated ‘Hominid Palaeoproteomic Reference Dataset’ available at Zenodo⁷⁰. In addition, previously published proteomes from extinct archaic hominins (H. antecessor and H. erectus)^32,41, Paranthropus⁷¹ and Gigantopithecus⁷² were incorporated. The search parameters were the same as those used in the initial search. Peptides detected in the negative control sample or assigned to contaminant proteins were excluded. The remaining peptides from each acid etch fraction were then reported and used to calculate deamidation profiles and validate their endogeneity⁶⁷ (Supplementary Text 2 and Supplementary Fig. 19).

Phylogenetically informative SAPs validation

Genus- and population-specific variants were screened using a previously established pipeline and database, excluding variants that had single-peptide support or were indistinguishable owing to PTM⁵. Variants within the Homo genus that met initial thresholds (PSM count greater than or equal to 2, PSM ratio greater than or equal to 10%, and intensity ratio greater than or equal to 10%) were further validated using the strict filter criteria described in previous publications to confirm their reliability^4,5,71. The peptides and PSMs supporting these sites were manually checked using the following criteria.

1. (1)

   Basic Local Alignment Search Tool (BLAST): all the supporting peptides were submitted to BLAST search against the NCBInr database to determine their specificity and endogeneity. Peptides matching any archaea, bacteria or fungi sequences were removed from the analysis. In addition, if a peptide matched any proteins from other genes, it was considered to be non-unique and excluded.

2. (2)

   Depth: after BLAST filtering, variants covered by at least two peptides were retained for the following evaluations.

3. (3)

   MS2 support: the MS/MS spectrum for each supporting PSM was manually inspected to check the fragmentary y and b ions surrounding the SAP. Two levels of filtering were applied: the MS2 support criterion required one y ion or one b ion to be associated with the SAP, whereas the strict MS2 support criterion required at least two y ions or two b ions flanking the SAP. When the SAP was situated at the terminus of the peptide, PSMs with only one associated y or b ion were also considered to meet the strict MS2 support criterion. The number of peptides with MS2 support PSM and strict MS2 support PSM was calculated.

4. (4)

   Deamidation rate: if at least one supporting PSM contained Q or N residues, deamidation values of this SAP were calculated on the basis of spectrum count and ion intensity. The quantities of PSMs and covered N and Q sites were also calculated. An elevated deamidation rate of Q or N residues in the supporting peptides indicates a probable endogenous origin and thus supports the reliability of the identified SAP.

5. (5)

   Codon degeneracy: if the codon of the derived amino acid variant differed from the codon of the ancestral allele by two or more bases, it was assumed to be an artefact, as such mutations are rare, and the variant was excluded.

6. (6)

   Collagen-specific features: for SAPs that originated from collagens, the supporting peptides were checked to confirm that appropriate collagen-specific features were present. If a SAP is in the triple-helix region of the collagen sequence, the sequence usually consists of GXY repeats, with glycine appearing every three amino acids and the hydroxylation typically occurring at the Y position of proline within the GXY repeats. These two features are important for the stability of the triple-helix structure of collagen^73,74,75,76. If the amino acid variant disrupted the original GXY repeats or the hydroxylation pattern was inconsistent with known biological specificity, the SAP was considered to be less reliable.

7. (7)

   Validation using multiple software tools: raw data files were searched against the same database used for the PEAKS search using two other software tools, MaxQuant (v.2.6.0.0)⁷⁷ and pFind (v.3.2.1)⁷⁸. One exception was the data search for ‘acid-insoluble fraction elastase’ from the parietal BFD767 in MaxQuant. As using ‘unspecific’ enzyme specificity in MaxQuant was time-consuming and the recovered endogenous proteome for this sample comprised only collagens, collagens were extracted from the Hominidae bone database to construct a smaller and more suitable database for this data search. For the acid etch fraction, the search parameters were the same as those used in PEAKS. For the dentine and bone fractions, minor modifications of the search parameters were used: in the pFind search, the peptide length range was set to 8–60 with a peptide mass range of 440–10,000 Da, and open search was enabled. In MaxQuant, the ion mass tolerance and maximum number of modifications per peptide were left as the default values. Pyro-Glu from E and pyro-Glu from Q were removed from the variable modification list, and oxidation K (hydroxylation K) was added. In both pFind and MaxQuant, the FDR was set to 1% at spectra level and 10% at protein level.

On the basis of these criteria, 12 positions were evaluated in the 5 hominin samples, and 6 positions were confidently identified (Supplementary Text 3 and Supplementary Data 4).

Construction of the protein consensus sequences and phylogenetic analysis

Consensus protein sequences were generated for the endogenous proteins identified in the five hominin samples (Supplementary Data 6). Heterozygous positions were determined on the basis of PSM count (greater than or equal to 2), PSM ratio (greater than or equal to 10%) and intensity ratio (greater than or equal to 10%)^5,71. The reliability of these positions was further confirmed using the SAP validation strategy described above, and alleles considered to be unreliable were excluded from further analysis. At reliable heterozygous positions, we randomly selected one variant to build the consensus sequence. Positions with single-peptide coverage or unreliable SAPs were masked as ‘X’ in the consensus sequences. One protein from the parietal BFD769 was removed because each position was supported by only a single peptide, resulting in five consensus protein sequences for this sample. For YHB3518, two proteins (COL1A1 and COL1A2) were identified in both enamel and dentine. To increase the coverage, the peptides for each protein were combined to generate its consensus sequence.

To determine the phylogenetic placement of the five hominin samples, we constructed a dataset for each sample on the basis of its endogenous proteome. We used 13, 16, 6, 5 and 5 proteins for tooth YHB3518, tooth YHB3075, parietal BFD767, parietal BFD769 and radius BFD771, respectively. In addition to the Bianfu Cave sequences, the reference dataset included three archaic individuals (Denisova 3 and Neanderthals from Altai and Vindija19), one chimpanzee (Pan troglodytes), one pygmy chimpanzee (Pan paniscus) one western gorilla (Gorilla gorilla) and one Sumatran orangutan (Pongo abelii, as an outgroup), whose protein sequences were chosen from the Hominidae bone database⁵. The protein sequences for H. sapiens (modern human) were obtained from SwissProt. These protein sequences were aligned and concatenated, resulting in datasets with lengths of 13,737, 13,701, 8,194, 7,654 and 7,654 amino acid positions for tooth YHB3518, tooth YHB3075, parietal BFD767, parietal BFD769 and radius BFD771, respectively. Leucines were converted to isoleucines in the datasets because they are isobaric amino acids and cannot be distinguished by tandem mass spectrometry. The best partition schemes and substitution models were chosen using Partition Finder (v.2.1.1)⁷⁹; these were then used in MrBayes (v.3.2.6)⁸⁰. The other parameters for the Bayesian analysis were as follows: two Markov chain Monte Carlo runs, 10,000,000 generations sampling every 5,000 generations and 25% burn-in. Convergence was evaluated using the average standard deviation of split frequencies (less than or equal to 0.01) and effective sample sizes (greater than or equal to 200). A 50% majority-rule consensus tree was generated, and clade frequencies were used to represent posterior probabilities.

Micro-CT scanning

The two parietal bones and one radius fragment were scanned using a Phoenix v|tome|x m microfocus CT system housed at IVPP, CAS. The scanning parameters were as follows: 160 kV, 130 µA. The isometric voxel size was 34.96 µm. Tomographic slices in TIFF format were imported into Mimics v.17.0 (Materialise) for virtual reconstruction of the specimens.

Cranial vault thickness and distribution pattern of the parietals

The linear thickness of the parietal fragments from Bianfu Cave was measured approximately at the centre using VGSTUDIO MAX 2023.1 (Volume Graphics). The distribution pattern of parietal thickness was visualized as a colourmap generated using the Surface Distance function in Avizo 9.2 (Thermo Fisher Scientific). Comparative samples included Asian late H. erectus, East Asian late Middle Pleistocene archaic Homo (Maba, Xuchang, Dali, Xujiayao, Jinniushan), H. heidelbergensis, Neanderthals and modern humans. The linear thickness data for comparative specimens were sourced from the literature (Supplementary Table 3), and the three-dimensional (3D) models of comparative specimens used for the colourmap analysis were generated by surface scanning of high-resolution models or micro-CT scanning of the specimens housed at IVPP, CAS.

CSG analysis of the radial mid-neck and diffeomorphic surface matching analysis of the proximal radius

Radial models were aligned to a standard orientation in Avizo 9.2 following the protocol established by Ruff⁸¹ to ensure precise anatomical correspondence before CSG analysis. We used the morphomap R package⁸² to extract the subperiosteal surface and calculate CSG properties at the radial mid-neck, defined as the midpoint between the distal margin of the radial head and the proximal limit of the radial tuberosity. CSG properties were calculated exclusively from the subperiosteal surface.

We used diffeomorphic surface matching to quantify morphological variation in the proximal radius. To ensure precise anatomical correspondence between BFD771 and the comparative sample, we truncated the diaphysis distal to the most prominent projection of the radial tuberosity for all specimens. The external surfaces were subsequently aligned through rigid and uniform scale transformations in Avizo 9.2 to eliminate variations in orientation and size. A reference template representing the mean shape of all specimens and a group of transformation parameters were generated from the aligned surfaces using Deformetrica v.4.3 (ref. ⁸³). In this framework, the transformation is parameterized by a grid of control points distributed in the ambient space, with geometric deformations governed by momenta vectors associated with these points^84,85. Principal component analysis of deformation-based shape residuals generated from control points and moment was conducted using the RToolsForDeformetrica package v.0.1 in R v.4.1.0. The extreme conformations of the first two principal component axes (PC1 and PC2) were generated using Deformetrica v.4.3 (ref. ⁸⁴).

The comparative framework for the radial morphological analyses comprised a diverse assemblage of hominins and extant great apes. The modern human sample consisted of cemetery-derived individuals (n = 14, pooled sex) housed at IVPP, CAS. Fossil hominin specimens include Neanderthals (Feldhofer 1, La Ferrassie 1, Krapina 189R. 1, 191R. 3, Spy 6), H. cf. erectus (SK 18b), A. afarensis (A.L. 288-1), A. sediba (UW88-85) and P. boisei (KNM-ER 1500). Finally, the extant great ape sample included pooled-sex representatives of Pan (n = 8), Gorilla (n = 9) and Pongo (n = 7).

Virtual data for the modern human sample were acquired using a Phoenix v|tome|x m microfocus CT system housed at IVPP, CAS. CT data for the bilateral radii of La Ferrassie 1 were provided by the Muséum national d’Histoire naturelle following a formal institutional request. Digital surface models for the radii of A. sediba (MH 2) and all extant great apes were obtained from the MorphoSource digital repository. Essential specimen information is summarized in Supplementary Table 4, and all citations and acknowledgments conform strictly to the MorphoSource Standard Usage Agreement^86,87. Digital surface models of the remaining fossil hominin specimens were generated by 3D scanning of the permanent fossil cast collections housed at IVPP, CAS.

Ethics statement

The protein analyses of five hominin specimens were approved by the Yunnan Provincial Institute of Cultural Relics and Archaeology, and the researcher responsible is a co-author.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.