Agta demography

Ethnographic data collection took place over two seasons in April–June 2013 and February–October 2014. We censused 915 Agta individuals (54.7% male) across 20 camps. The four selected camps and two multi-camps in the current study were the ones providing social interaction data and saliva samples from almost all camp members. Accurate ages were estimated following relative aging protocols (Diekmann et al., Reference Diekmann, Smith, Gerbault, Dyble, Page, Chaudhary and Thomas2017). Relatedness (biological and affinal) was based on household genealogies. To resolve inconsistencies, we took either the genealogy from the most knowledgeable individual (i.e. mother over aunt) or the genealogy that reduced other inconsistencies (i.e. discarding 6 month interbirth intervals). Genealogies contained 2953 living and dead Agta. We used the R packages pedigree, kinship2 and igraph to measure consanguineous relatedness (r) (Dyble et al., Reference Dyble, Salali, Chaudhary, Page, Smith, Thompson and Migliano2015; Page et al., Reference Page, Myers, Dyble and Migliano2019). For comparative purposes, we obtained 14 saliva samples from neighbouring Palanan farmers, making sure that individuals were unrelated by directly asking.

BaYaka demography

Ethnographic data collection took place over two seasons in April–June 2013 and February–October 2014. We collected saliva samples from 21 individuals for microbiome analyses.

Ethics

This study was approved by UCL Ethics Committee (UCL Ethics code 3086/003) and carried out with permission from local government and community members. Informed consent was obtained from all participants, after group and individual explanation of research objectives in the indigenous language. A small compensation (usually a thermal bottle or cooking utensils) was given to each participant. The National Commission for Indigenous Peoples (NCIP) advised us that the process of Free Prior Informed Consent should be obtained from the community leaders, youth and elders under the supervision and validation of NCIP. This was done in 2017 with the presence of all community leaders, elders and youth representatives at the NCIP regional office, with the mediation of the regional officer and the NCIP Attorney. The validation process was approved unanimously by the community leaders and the NCIP, and validated the full 5 years of data collection.

Microbial DNA extraction and 16S rRNA gene sequencing.

A total of 190 saliva samples were selected from Agta hunter–gatherers (n = 155) and Palanan farmers in the Philippines (n = 14), and BaYaka hunter–gatherers from the Congo (n = 21). Microbial DNA was extracted following the protocol for manual purification of DNA for Oragene⋅DNA/saliva samples. The 16S rRNA gene V3–V4 region was amplified by PCR with primers containing Illumina adapter overhang nucleotide sequences. All PCR products were validated through an agarose gel and purified with magnetic beads. Index PCR was then performed to create the final library also validated through an agarose gel. All samples were pooled together at equimolar proportions and the final pool was qPCR-quantified before MiSeq loading. Raw Illumina pair-end sequence data were demultiplexed and quality-filtered with QIIME 2 2019.1 (Bolyen et al., Reference Bolyen, Rideout, Dillon, Bokulich, Abnet, Al-Ghalith and Caporaso2019) and DADA2 (Callahan et al., Reference Callahan, McMurdie, Rosen, Han, Johnson and Holmes2016), which generates single nucleotide exact amplicon sequence variants (ASV or ESV). ASVs are biologically meaningful as they identify a specific sequence and allow for higher resolution than operational taxonomic units (OTUs) (Callahan et al., Reference Callahan, McMurdie and Holmes2017) or clusters of sequences above a similarity threshold, and thus an ASV is equivalent to a 100% similar OTU. Taxonomic information was assigned to ASVs using a naive Bayes taxonomy classifier against the SILVA database release 132 with a 99% identity sequence (Quast et al., Reference Quast, Pruesse, Yilmaz, Gerken, Schweer, Yarza and Glöckner2013).

Reads outside the kingdom Bacteria or assigned to mitochondria or chloroplasts were removed. Phylogenetic analyses aligned sequences with MAFFT (Katoh et al., Reference Katoh, Misawa, Kuma and Miyata2002) and generated a rooted phylogenetic tree with FastTree2 (Price et al., Reference Price, Dehal and Arkin2010) using default settings via QIIME 2. We generated an Alpha rarefaction curve with R package vegan to confirm that sample richness had been fully observed (Supplementary Figure S3).

Samples with extremely low number of reads (8000) were removed. This resulted in 6409 ASVs (later reduced to 1980 ASVs present in at least two individuals and with an abundance of at least 10 counts per individual) and 173 individuals: 138 Agta, 21 BaYaka and 14 Palanan farmers.

Identification of the Agta social microbiome

In our Agta sample, we first selected a set of strong social links (the top 25% of the weight distribution from each camp and multi-camp). For each ASV, we calculated the proportion of strong links ($f_{\rm A}^{\rm s}$) where a given ASV A was present. Next, we calculated the same proportion in the complementary set of 75% weak social links ($f_{\rm A}^{\rm w}$). We then computed for each ASV A the score

$$s_{\rm A}^{\rm s} = \displaystyle{{\,f_{\rm A}^{\rm s} -f_{\rm A}^{\rm w} } \over {\,f_{\rm A}^{\rm s} + f_{\rm A}^{\rm w} }}$$

or normalised difference between the two proportions. This score can be paired with the Z-score

$$\displaystyle{{\,f_A^s -f_A^w } \over {\sqrt {\,f_A( {1-f_A} ) ( \displaystyle{1 \over {n^s}} + } \displaystyle{1 \over {n^w}}) }}$$

which quantifies the deviation from the null hypothesis that the two proportions are equal, with n ^s and n ^w as respectively the numbers of strong and weak links and f _A as the proportion of total links that share ASV A. We then selected as affected by social interaction 137 ASVs with s > 0.5 and p < 0.05. The p-values were adjusted by false discovery rate.

We performed a sensitivity analysis to investigate the consequences of varying the threshold defining strong vs. weak social links in the Agta social network. Instead of 137 ASVs resulting from selecting dyads at top 25% of the weight distribution, we obtained:

top 45% – 166 ASVs;

top 35% – 156 ASVs;

top 25% – 137 ASVs;

top 15% – 138 ASVs;

top 5% – 99 ASVs.

The list shows that whether the strong set consists of a very reduced number of dyads with very strong weights (top 5%), or instead includes nearly half of the dyads (top 45%), the number of ASVs significantly responsive to social contact varies from 99 (5%) to 166 (8%), representing a small fraction of the total of 1980 ASVs found in the Agta. Therefore, setting the threshold at the top 25% did not affect our results and conclusions.

Portable radio sensor devices

Our wireless sensing devices store all between-device communications within a specified distance and have been described in detail elsewhere (Migliano et al., Reference Migliano, Page, Gómez-Gardeñes, Salali, Viguier, Dyble and Vinicius2017, Reference Migliano, Battiston, Viguier, Page, Dyble, Schlaepfer and Vinicius2020; Page et al., Reference Page, Chaudhary, Viguier, Dyble, Thompson, Smith and Migliano2017, Reference Page, Myers, Dyble and Migliano2019). We used the UCMote Mini (with a TinyOS operating system) sealed into wristbands or belts, labelled with a unique number and identified with coloured string to avoid accidental swaps. The devices require no grounded infrastructure and collect interactions even when individuals are away from camps. Individuals arriving at a camp after the start of data collection were given a device and the entry time was recorded, while those leaving a camp before the end of data collection had their exit time recorded. To prevent swaps individuals were checked twice daily, and device numbers were checked upon return. Any swaps were later corrected by reassigning data to the correct individuals.

Data were later downloaded via a PC side application in Java. Data were limited to 5 a.m. to 8 p.m. We ran raw data through a stringent data-processing system in Python to prevent data corruption. Data were matched to ID numbers and start-stop times of each sensor. The result was a matrix with the number of recorded beacons for all possible dyads and their weights.

For the camp-level experiment, all individuals from four camps wore sensors from 5 to 7 days. Each device sent a message every 2 minutes that contained its unique ID, a time stamp and the signal strength. Messages are stored by any other mote within a 3 m radius, a frequently used threshold (Brent et al., Reference Brent, Lehmann and Ramos-Fernández2011; Isella et al., Reference Isella, Stehlé, Barrat, Cattuto, Pinton and Van den Broeck2011). Detection of such close-range interactions is our proxy for social interactions. In previous studies, we have provided evidence that our sensor data capture interactions resulting from joint activities such as foraging, gathering, child care and socialising among others, through validation of sensor data by parallel observations obtained during camp scans (Migliano et al., Reference Migliano, Page, Gómez-Gardeñes, Salali, Viguier, Dyble and Vinicius2017; Page et al., Reference Page, Chaudhary, Viguier, Dyble, Thompson, Smith and Migliano2017). For the multi-camp experiment, adult individuals from two areas (consisting of seven and three camps respectively) wore sensors for 1 month.

Effect of dyad category on bacterial sharing

The bacterial sharing network was constructed by defining link weights as the number of social bacteria shared by two individuals. Dyads in the network were classified into: (i) levels of kinship (mother–offspring, father–offspring, siblings, r = 0.5; other kin, r = 0.25 or r = 0.125; non-kin, r = 0.0625; or lower, spouses, friends, defined as non-kin at the top 25% distribution of social dyadic weights, and other non-kin); and (ii) residence (same or different camp, same or different household). The mean weights were calculated for each dyadic category (Supplementary Table S2). Next we produced 1000 network randomisations based on a single-step ID swap between nodes. For example, if dyad 1 consisted of two spouses in the real network, randomisation preserved dyad 1 and its weight, but randomly replaced the two nodes (potentially changing the dyadic classification to siblings, friends, and so on). We calculated the mean weights for each dyadic category in the 1000 randomised networks, and then calculated one-sample t-tests with the mean weight in the real network as the test value. We repeated the analysis for non-social bacteria (Supplementary Table S3).

Reinforcement analysis

In a multilayer network, reinforcement analysis measures the overlap in links between different layers to quantify the probability of finding a link on a layer conditioned on the weight of the same link on another (Battiston et al., Reference Battiston, Nicosia and Latora2014). Reinforcement between two network layers α and α′, or P(α ^′∨α), is defined as

(1)$$P( {\alpha^{\prime}\vee \alpha } ) = \displaystyle{{\mathop \sum \nolimits_{ij} a_{ij}^{[ {\alpha^{\prime}} ] } w_{ij}^{[ \alpha ] } } \over {\mathop \sum \nolimits_{ij} w_{ij}^{[ \alpha ] } }}$$

where $a_{ij}^{[ {\alpha^{\prime}} ] }$ is the adjacency matrix of conditioned layer α′, and $w_{ij}^{[ \alpha ] }$ is the adjacency matrix of the conditioning layer α. We split the weight of social links into three tertiles, and computed equation (1) for each. We obtained increasing values of reinforcement from the lower to the higher tertile, providing evidence of an effect of social contact on bacterial sharing.

ASV diversity metrics

To distinguish between the effects of lifestyle and shared ecology on the microbiome, we compared the diversity of the oral microbiome of Agta hunter–gatherers with neighbouring Palanan farmers in the Philippines and BaYaka hunter–gatherers in Congo. Using the 6409 ASVs dataset we calculated the number of observed ASVs in each population with R package Phyloseq (version 1.30.0) and Faith's Phylogenetic Diversity index was calculated with R package picante using the generated rooted phylogenetic tree. To estimate the number of shared ASVs in the Agta, BaYaka and Palanan farmers, we sampled a random subset of 10 samples for each population without replacement, calculated the shared ASVs between the populations and repeated this procedure 100 times. Global differences between groups and pairwise comparisons were assessed by Kruskal–Wallis and Wilcoxon rank sum tests respectively and plotted by the R package ggpubr. Pairwise p-values were adjusted by false discovery rate.

Classification of oral bacteria as pathogens

ASVs were classified as oral pathogens if they have been reported as etiological agents of periodontitis (Pérez-Chaparro et al., Reference Pérez-Chaparro, Gonçalves, Figueiredo, Faveri, Lobão, Tamashiro and Feres2014; Socransky et al., Reference Socransky, Haffajee, Cugini, Smith and Kent1998) or dental caries (Simón-Soro et al., Reference Simón-Soro, Guillen-Navarro and Mira2014; Simón-Soro & Mira, Reference Simón-Soro and Mira2015). For gum disease pathogens, these included the classical ‘red’ and ‘orange’ complex of periodontal pathogens and the recent update (Pérez-Chaparro et al., Reference Pérez-Chaparro, Gonçalves, Figueiredo, Faveri, Lobão, Tamashiro and Feres2014) based on systematic review and metanalysis. For caries pathogens, the list of active microorganisms detected through metatranscriptomics of cavities was used. Common oral commensals potentially causing endocarditis or systemic infections in immunocompromised patients only were not considered pathogens. Bacteria reported as etiological agents of lower respiratory infections (e.g. pneumonia, whooping cough, bronchitis or sinusitis) and biofilm-mediated infections (e.g. lactational mastitis, medical implant biofilm infections, chronic lung infections, osteomyelitis or chronic wounds) were also considered pathogens, including organisms present in healthy carriers (Leung & Hon, Reference Leung and Hon2018; Natsis & Cohen, Reference Natsis and Cohen2018). Bacteria causing urinary tract infections or sexually transmitted diseases transiently found in the oral cavity were also considered pathogens (Jung et al., Reference Jung, Ehlers, Lombaard, Redelinghuys and Kock2017). Bacteria were classified as oral if detected in more than 10% of the population in oral samples according to the Human Oral Microbiome database. If a bacterial species or genus had been isolated from the oral cavity of another animal species, it was also classified as oral.

For assignment of bacteria as pathogenic or non-pathogenic, we used species-level ASVs, given the multiple cases where species from the same genus had a different assignment. If taxonomic classification of an ASV was only possible at the genus level, it was considered a pathogen if: (i) >90% of named species within the genus were pathogenic; or (ii) the genus included a major pathogenic species even if the remaining species were not classified as oral by the Human Oral Microbiome Database (Chen et al., Reference Chen, Yu, Izard, Baranova, Lakshmanan and Dewhirst2010). ASVs with a top hit to a sequence classified as ‘Oral taxa’ in databases but without a species assignment were discarded from the analysis. Cases where taxonomic classification of the ASV was only possible at the family level or higher were also discarded.