Participants
The sample consisted of 67 triads: heterosexual women (mean age = 24.2 years, SD = 3.3), their fathers (mean age = 51.5 years, SD = 6.3) and romantic partners (mean age = 26.7 years, SD = 4.4). The required sample size was estimated based on a priori power analysis in which, based on previous studies, we assumed an association between partner’s and opposite-sex parent’s characteristics to be equivalent to a positive correlation of value 0.25 (for details, see the study’s preregistration at https://osf.io/adrtc). Only women who met the following criteria were enrolled in the study: i) aged 18-35 years, ii) currently in a committed long-term (for a minimum of 6 months) heterosexual romantic relationship (actual mean length = 59 months, SD = 45.3, range = 8-233); iii) shared a common household with their biological father until at least 8 years of age (actual mean = 20.9, SD = 3.6, range = 8-30 years); and iv) not suffering from conditions affecting their sense of smell which may in turn influence their olfactory preferences31. The enrolment criteria for fathers and romantic partners were: i) aged under 65 years (fathers) or 18-40 years (partners) ii) non-smokers; iii) in good health and having no underlying condition which affects body odour (e.g., diabetes)30; iv) not currently using medication; and v) having unshaven axillae32.
Participants were recruited via an e-mail contact list compiled from our previous studies, social media (Facebook, Instagram), and flyers advertised at Charles University, Prague. To check whether participants met the enrolment criteria (e.g., age, health conditions), they completed an online contact form via the Qualtrics platform. Subsequently, we arranged a personal meeting with one individual from each triad and handed over the package with material for data collection. During this meeting, participants were informed about the aims of the study and the sampling procedure was fully explained by going over the written instructions which accompanied the package.
Procedure
Within the space of one week, each woman completed questionnaires to assess the quality of the relationship with her father during childhood (s-EMBU and overall assessment of childhood relationship quality) and her current attachment style (ECR-RS). Fathers and partners were asked to provide body odour samples for perceptual ratings, chemical analysis, and axillary microbiome analysis. Participants also provided facial photographs and voice recordings which are not the subject of this paper and will be reported elsewhere.
Data collection was carried all over the Czech Republic (32 villages, towns or cities, Fig. 1a), increasing the variability of socio-demographic characteristics of our sample (Extended data Table 1). Participants confirmed their participation by signing a consent form. Each triad received compensation for their participation in the sum of 1200 CZK (app. 44 EUR). The study was preregistered (https://osf.io/adrtc) and was approved by the Institutional Review Board of Charles University, Faculty of Science (no. 2017/20).
Body odour collection
Each father and partner obtained a package including printed instructions and materials for sample collection (non-perfumed liquid soap, tissues, six sterile cotton swabs, three glass vials, three Eppendorf tubes, a vial with the sterile physiological solution, three cotton pads separately enclosed in the labelled zip-lock plastic bags, surgical tape, white cotton T-shirt and a polystyrene transport box (15x20x10 cm). Participants were asked to refrain from consuming aromatic food (a list was provided), drinking alcohol, smoking or using drugs, using cosmetics (e.g., perfumes and deodorants) and from strenuous physical activities including sex during the 48 hours preceding sample collection. All these activities may influence body odour quality33. After collection of the samples, each participant completed a survey to document rule compliance (Supplementary Table 1). Although there were some minor violations two days before sampling, there were none reported in the 24h period leading up to the sampling.
In the evening on the sampling day, participants collected axillary samples for chemical analysis and then for microbiome analysis. They first washed their hands using non-perfumed liquid soap and dried them with the tissues. To collect samples for chemical analysis, participants wiped each axilla 20 times with a clean cotton swab and placed each swab in a 2 ml empty glass vial. The same procedure was used for collecting axillary microbiome samples, except that the swabs were first soaked in sterile physiological solution. In both types of sampling, they used a sterile cotton swab, unique for each axilla. Axillary microbiome swabs were placed inside an Eppendorf vial containing absolute (99.8%) ethanol. Finally, two control samples (one each for chemical and axillary microbiome analysis) were collected by waving a clean cotton swab 20 times in the air); these samples were used to assess possible effect of environment (e.g., specific odour of the household). All these samples were put into the pre-labelled zip-lock plastic bags and odourless polystyrene box and placed in a freezer (-20°C) immediately after sampling (Fig. 1b).
Subsequently, the participants took a shower and washed with the provided non-perfumed liquid soap. They then attached a cotton pad to each axilla using surgical tape. To minimise odour contamination from clothing and other extrinsic ambient odours, participants wore a new white 100% cotton T-shirt as the first layer of their clothing. They wore the pads for the following 12 hours (overnight)34,35. Participants collected a further control sample by placing a clean cotton pad on their bedside table for the same period that they wore the axillary pads. Next morning, they placed the cotton pads in separate, pre-labelled zip-lock plastic bags and added these to the polystyrene box together with the samples for chemical and axillary microbiome analysis, which was then put back into the freezer36. They returned all samples to the experimenters within a week (in most cases, the experimenters collected these from the participant’s homes). For an overview of the entire procedure, see Supplementary Fig. 1.
Questionnaires
The quality of relationship with the father during childhood was assessed using the s-EMBU scale (an acronym for ‘Egna Minnen Beträffande Uppofostran’, Swedish for ‘My memories of upbringing’)37. This scale is a widely used tool for assessing adults’ perception of their parents’ rearing behaviour. It consists of 23 items, loading onto three sub-scales: Rejection (7 items; possible range 7–28), Emotional Warmth (7 items; possible range 7–28), and (Over)Protection (9 items; possible range 9–36). Responses were indicated on a 4-point Likert-type scale (ranging from 1 – ‘No, never’ to 4 – ‘Yes, most of the time’). We used a Czech language version used in previous studies38. We further used a single item to explore overall quality of the relationship during childhood (‘How would you evaluate your relationship with your father during childhood?’), indicated on a 7-point Likert-type scale (ranging from 1 – ‘Very negative’ to 7 – ‘Very positive’).
Adult attachment style between the woman and her father was assessed using the Experiences in Close Relationships – Relationship Structures (ECR-RS) questionnaire39. This consists of 9 items about attachment orientation of the daughter to father and comprises two subscales: Anxiety (3 items; possible range 3–21) and Avoidance (6 items; possible range 3–42). Daughters completed a Czech version of the scale using a 7-point Likert-type scale (ranging from 1 – ‘Strongly disagree’ to 7 – ‘Strongly agree’). Mean scores for each questionnaire are reported in Extended data Table 1.
Perceptual Ratings
Raters
In total, 306 female raters (mean age = 24.0, SD = 4.0) took part in the rating sessions. They were recruited via an e-mail contact list compiled from our previous studies, social media (Facebook, Instagram) and flyers or by oral invitation around Charles University and the Czech Technical University. Only women in good respiratory health and without any olfactory malfunction took part in the rating sessions. Each rater could participate only in a single rating session. At the end of the session, raters were informed about the study goals and obtained 150 CZK (approximately 5.5 EUR) as compensation for their time.
Match-to-sample tests
To assess perceived similarity between father and partner body odour we used a match-to-sample paradigm. In total, 67 body odour sets were created. Each set consisted of a father’s body odour sample (the ‘target’), the partner’s body odour sample (the ‘match’) and body odour samples from three other individuals (‘distractors’). The distractor body odours were collected from independent male donors (N= 89, mean age = 25.8; SD = 5.4) who met the same conditions as partners (e.g., the same age, being non-smokers, having unshaven axillae) and who followed the same body odour sampling procedure as described above. Each set consisted only of samples collected from one axillary side (i.e., left/right). In total, 28 sets consisted of left axillary samples and 39 sets used right axillary samples. We sometimes also used partners’ body odour samples as distractors, but these were always the other axillary sample (in other words, if a sample from the left axilla was used as a match, the right axilla sample was used as a distractor).
In each session (13 sessions in total), 5 sets of stimuli were used (except the last two sessions consisting of 6 sets). The initial position of the match and the distractors was randomized for each rater (Fig. 1c). The match sample position was then rotated clockwise for each subsequent rater. The distractors used for each rater were randomly chosen from the pool of 15 (or 18 in the last two sessions) distractors available for the given rating session. Similarly, the order of the individual sets was randomized for each rater. Samples were randomized by a blinded experimenter who followed a randomization matrix generated before the rating session.
Rating session
Rating sessions took place in a quiet, ventilated room with a mean temperature of 22.5°C. The body odour samples were thawed one hour before the rating session36. Subsequently, the samples were enclosed in 250ml opaque labelled jars. Raters provided informed consent and were informed about the rating task. The origin of the body odour samples was not disclosed.
First, raters were asked to evaluate similarity of the body odours within each set (‘Please sniff the target and then the four other samples. Order the four samples from the most similar to the least similar as compared to the target’). Raters marked their scores into tablet using a purpose-built application MatchApp. The sniffing time was not restricted, most raters finished the odour rating task in 30 minutes. The body odour ratings were followed by assessments of facial photographs and voice recordings, results of which are not relevant to this study.
Statistical analysis of perceptual ratings
The perceived father-partner similarity in body odour was first evaluated using Chi-square tests - the perceived similarity rankings for each father-partner dyad were compared to the null hypothesis of uniform random distribution of the ranks (i.e., that the match (partner) should appear as the most similar of the four samples in 25 % of cases). Next, to analyse the relative frequency of the match appearing in each of the four ranks, while accounting for non-independence of individual rankings and target’s identity we employed an ordinal probit bayesian mixed model (fitted using R package brms)40 where parameter estimates and 95% posterior credible intervals were reported. The identity of the raters, the identity of the targets, and axillary side (left/right) were all included as random effects. For a subset of the partner samples, there was no difference between left and right axilla in similarity ratings (estimate [± 95% posterior credibility intervals] = 0.04 ± [-0.43~ 0.50]).
The potential effect of relationship quality with the father on perceived father-partner odour similarity was analysed by non-parametric paired Spearman’s correlation (Extended data Table 2). Correlations between individual measures of the relationship quality were generally moderate (Spearman's rho < 0.4) except for the strong negative association between the variables Overall relationship quality and EMBU Rejection (Spearman's rho = -0.59) and the positive association of Overall relationship quality and EMBU Emotional Warmth (Spearman's rho = 0.57) (Supplementary Table 3).
Axillary microbiome analysis
Metagenomic DNA was extracted using PowerSoil (Qiagen) kits and sequencing libraries were prepared using a two-step PCR approach, following Glenn et al.41. Gene-specific primers covering the V3-V4 variable region of bacterial 16S rRNA (i.e., S-D-Bact-0341-b-S-17 [CCTACGGGNGGCWGCAG] and S-D-Bact-0785-a-A-21 [GACTACHVGGGTATCTAATCC]) were used during the first PCR step42. The primers were extended at the 5’ end with “tails” serving as a priming site for the second PCR (i.e., the first 33 bp and 34 bp of read1 and read2 Nextera adapters for F and R primers, respectively). The first PCR was performed in 10 μl and consisted of 5 μl 1x KAPA HiFi Hot Start Ready Mix (Roche), each primer at 0.2 μM and 4.6 μl DNA. PCR conditions were as follows: initial denaturation at 95°C for 3 min followed by 35 cycles each of 95°C (30 s), 55°C (30 s), and 72°C (30 s), and a final extension at 72°C (5 min). Dual-indexed Nextera sequencing adapters were reconstructed during the second PCR, which followed the first PCR conditions except that: it was performed in 20 μl volume; the concentration of each primer was 1 uM; 6μl of the first PCR product were used as a template; and the number of PCR cycles was 12. Technical duplicates were prepared to account for noise due to PCR and sequencing stochasticity.
Products of the second PCR were run on 1.5% agarose gel and their concentrations were assessed based on gel band intensities using GenoSoft software (VWR International, Belgium). They were pooled equimolarly and size-selected with Pippin Prep (Sage Science, USA) at 520 - 750 bp. Resulting libraries were sequenced using MiSeq (Illumina, USA) and v3 chemistry (i.e., 2 × 300 bp paired-end reads). Sequencing data are available at the European Nucleotide Archive under the accession number of the whole project (will be completed upon the acceptance of the manuscript).
Samples were demultiplexed and primers were trimmed by skewer software43. Using dada244, we eliminated low-quality sequences (expected number of errors > 2), denoised quality-filtered fastq files, and constructed abundance matrix representing read counts for individual amplicon sequencing variants (hereafter ASV) in each sample. Next, we identified chimeric ASVs using uchime45 and gold.fna reference database and eliminated them from the abundance matrix. Taxonomic assignation of ASVs was conducted by RDP classifier (80% confidence threshold46) and Silva database (v. 132)47.
Technical replicates of axillary microbiota samples exhibited high consistency in both composition (Procrustean correlation, r = 0.99, p = 0.0001) and Shannon diversity (Pearson's r = 0.93, p < 0.0001). Consequently, we merged sequences corresponding to individual samples for all subsequent analysis. To avoid microbial diversity inflation due to PCR and sequencing artifacts, we eliminated all OTUs that were present just once within the set of PCR replicates for a given sample48. Prior to statistical calculations, we also excluded one sample that exhibited low consistency among its technical replicates and one sample of low sequencing coverage (< 1500 high quality sequences). The final dataset included 246 samples of left and/or right axillary microbiota coming from 132 individuals (66 fathers and 66 partners). These samples were represented by 3,129,530 high-quality reads that were distributed among 4,891 ASVs. The median sequencing coverage was 13,152 reads per sample (range = 1,750 – 34,343). Functional predictions of bacterial metagenome were conducted by PICRUSt2 pipeline49 using the default setup. Predicted metagenomes were then categorized to functional pathway50 and their predicted abundances in each sample were used in later statistical analyses. Weighted Nearest Sequenced Taxon Index (weighted NSTI) was 0.02755007 +/- 0.001274799 (+/- S.E) suggesting very high precision of metagenome predictions.
In addition to the axillary microbiota samples, we also collected 74 environmental microbiota samples. Due to lower DNA concentrations and PCR stochasticity, duplicates from the environmental samples had lower alpha diversity (r = 0.74, p < 0.0001) and compositional consistency (r = 0.68, p < 0.0001) compared to duplicates from the axillary microbiota, and often a lower number of reads. Therefore, we focused our analyses on 50 ambient environmental microbiota samples that had >1000 reads after all filtering steps.
Microbial alpha diversity did not significantly vary between left and right armpit (Generalized Linear Mixed Model (GLMM): Δ D.f. = 1, χ2= 1.31, p = 0.252 for ASVs richness and Δ D.f. = 1, χ2= 1.31, p = 0.252 for Shannon diversity) and there was also no significant difference in the number of detected ASVs between fathers and partners (GLMM: Δ D.f. = 1 χ2= 2.06, p = 0.151). However, partners exhibited slightly lower Shannon diversity compared to fathers (GLMM: estimate = -0.24 +/- 0.12, Δ D.f. = 1, χ2= 4.27, p = 0.039). NMDS (Extended Data Fig. 1a) and PERMANOVA analyses did not detect any systematic difference in the composition between left and right axilla, yet they revealed a slight compositional shift between fathers and partners explaining only < 2% of total variation in composition (Extended Data Table 3a and Extended Data Fig. 1b). Similarly, alpha diversity was closely correlated in the left and right axilla (Extended Data Fig 1c). According to GLMM analyses with distances to group-specific centroids as a response, profiles of left and right axilla were equally dispersed, but fathers exhibited higher dispersion in composition (i.e., inter-individual variation) than partners (Extended Data Table 3b and Extended Data Fig. 1d).
Samples from both left and right axillae were successfully profiled in the case of 114 individuals (58 fathers and 56 partners). Based on this subset, there was a tight correlation in alpha diversity between left and right axillae (Pearson’s correlation: r = 0.796, p < 0.0001 for ASVs richness and 0.797 for Shannon diversity). The same holds for consistency in the composition, where microbiota divergences calculated using ASVs presence/absence was more decisive (Procrustean correlation for Jaccard distances: r = 0.94, p = 0.0001) than divergence metrics utilizing relative ASVs abundances (Procrustean correlation for Bray-Curtis distances: r = 0.82, p = 0.00001, Fig. 2b and Extended Data Fig. 1c). Furthermore, consistency between left and right axillae was significantly higher for fathers than for partners (GLMM with Procrustean residuals as a response: Δ D.f. = 1, χ2= 5.44, p = 0.02 for Bray-Curtis and Δ D.f. = 1, χ2= 9.09, p = 0.003 for Jaccard dissimilarities, Extended Data Fig. 1d). For the subsequent analyses, profiles of left and right axillae of the same individual were aggregated.
Statistical analysis of microbiome data
Shannon indices and number of ASVs in each sample (hereafter ASVs richness) were calculated after dataset normalization by rarefaction (threshold = 1,750 reads per sample) and used as alpha diversity measures. We analysed differences in alpha diversity between fathers and partners and left vs. right axillae using GLMM with Gaussian error distribution. Individual identity nested within pair identity were specified as random effects. Pearson correlation was used to check for a correlation in alpha diversity between left and right axillae.
Differences in composition of microbiota profiles was analysed based on Bray Curtis dissimilarities accounting for ASV relative abundances and binary Jaccard dissimilarities accounting just for ASVs presence/absence. Both these metrices were calculated after the rarefaction of the ASV abundance matrix. First, we checked if there were any differences in composition between left and right axillae and between fathers and partners using Nonmetric Multidimensional Scaling (NMDS) and PERMANOVA with permutation constraint (i.e., ‘strata’) setup to within-pair level. To test if interindividual variation differed between left and right axillae, and between fathers and partners, we exacted distance to centroids for father’s and partner’s samples and used them after log transformation as a response in GLMM, assuming Gaussian errors and the same random structure as described above. ASVs overrepresented in fathers or partners were identified by differential abundance analyses, where read counts for each ASV were included as a response into GLMM with negative binomial errors and log-scaled total sequencing depth for a given sample as an offset. The model random structure remained unchanged. To achieve model convergence, these models were fitted using glmmTMB package51 only for ASVs detected in at least 10 samples. The Qvalue method52 was then used to control for false discoveries due to multiple testing, with the threshold being setup to q < 0.05.
The composition and alpha diversity of left and right axillary microbiota were tightly correlated at the individual level. We therefore aggregated the left and right axillary profiles for each individual and used the aggregated profiles to test if there was higher similarity between father-partner pairs than expected by chance. To do this, we compared average microbiota dissimilarity between actual father-partner dyads with a distribution of average dissimilarities drawn from randomly paired fathers and partners (n = 999 permutations). The outcome of these analyses were also expressed as Standardised Effect Sizes (SES)53 using the formula SES = (∑DIF(obs) – ∑DIF(sim)/nperm)/sd(∑DIF(sim)), where ∑DIF(obs) is the sum of observed within dyad dissimilarities, ∑DIF(sim) is the sum of simulated within-dyad dissimilarities, nperm is number of permutations, and sd is standard deviation. A simulation study was conducted to test the robustness of this permutation routine to false positives. First, using the R package AHMbook54, we simulated 1,000 random communities of 1,000 species mimicking father-partner dyads and thus consisting of two sample groups of 60 samples each. No a priori driver of correlation in community composition was specified within the 60 dyads. The frequency of false positive results estimated, based on these simulated datasets and the permutation routine described above, was low at 4.6% (Extended Data Fig. 1e). Analyses on predicted metagenome pathways utilized most of the above-described approaches, with two distinctions. First, instead of abundance matrix rarefaction, functional pathways abundances were normalized by their transformation to proportions in each sample. Second, we used only relative-based (i.e., Bray-Curtis) dissimilarities in all analyses on predicted metagenomes. All statistical calculations and data visualizations were conducted using R 4.0.255.
Chemical analysis
Sample extraction
The sampling for analyses of axillary chemical profiles was performed by the participants themselves as described above. The samples were transported and stored at -20°C until extraction, which took place as soon as possible after the reception of samples. Left and right axillary samples were successfully collected from all participants (67 father-partner dyads). To allow multiple analyses and storage of the irreplaceable samples, we selected the solvent extraction of the samples, in spite of the drawbacks inherent to this method when compared with direct thermal desorption, especially the lower representation of highly volatile and low molecular mass compounds. The analytes were extracted from the cotton swabs with 400 µl of distilled n-hexane for 24 hours at 8°C plus an additional 10 minutes of sonication. The extracts were transferred into clean vials, two internal standards (10 ng of 1-bromononane and 100 ng of 1-bromoeicosane) were added, and the extracts were either stored at -80°C or directly analysed.
Instrument setup
The samples were analysed using a GC×GC-TOFMS instrument Pegasus 4D (LECO Corp., St Joseph, MI, USA). The system is based on a quad-jet modulator with two cold liquid-nitrogen jets and two hot-air jets, which trap and refocus the analytes eluting from the first dimension column. The column set consisted of a nonpolar Rxi-5Sil MS (30 m, ID 0.25 mm, df 0.25 µm, Restek, Bellefonte, PA, USA) in the first dimension and a mid-polar Rxi-17Sil MS (1.5 ms, id 0.1 mm, df 0.1 µm, Restek, Bellefonte, PA, USA) in the second dimension. The columns were connected using SilTite® µ-Union (Trajan Scientific, Australia). The modulation period was 5 s and the hot pulse duration 600 ms. The temperature program was 50°C (1 min), 8°C/min ramp to 320°C (20 min). The secondary column was set 10°C higher. The modulator offset was 20°C. Helium was used as a carrier gas at a constant flow rate of 1 ml/min. Total run time was 54.75 min. 1µl of the extract was injected into a split/splitless injector heated to 200°C in splitless mode. Acquisition solvent delay was set at 400 s, while MS transfer line temperature was 280°C. The ion source temperature was 220°C and electron impact energy was -70 eV. Collected mass range was 29-600 amu and acquisition rate was 100 scans/s. The detector voltage was 1,500 V.
Measurement sequence
The samples were analysed in randomized order and injected using a 7683 Series Injector (Agilent Technologies). Before each sample, a blank was run using the same temperature program as for the sample runs. To guarantee the reliability of the chemical analyses and to avoid biases in chromatograms alignment, we only considered the samples which were analysed during a single measurement sequence, uninterrupted by any change in the analytical setup, such as column exchange, etc. The only intervention in the setup during the sequence measurement was the regular exchange of the split/splitless liner (Topaz liner, Splitless, Single Taper Gooseneck w/Wool, 4mm × 6.5 ×78.5, Restek), which was replaced after every 10 samples. The longest successful uninterrupted sequence included both left and right axillary samples of 42 father-partner dyads (168 samples), which were then used for data analysis.
Data processing
For peak detection, alignment and quantification, we processed the data generated by the operating software LECO ChromaTOF software v. 4.72 (LECO, St. Joseph, MI) with the development version of ChromaTOF-Tile (v. 0.27). This software applies the recently developed tile-based algorithm for GC×GC-TOFMS multiple chromatogram processing, which provides an effective alternative to the widespread pixel-based and peak table-based methods56,57. Tile-based algorithm offers two main advantages over the traditional approaches: it is computationally fast and it significantly reduces the problem of 2D retention time misalignments. ChromaTOF-Tile is primarily designed for supervised comparison of two or multiple pre-defined classes (groups) of samples. Therefore, for the purposes of unsupervised exploration of chemical profiles across all body odour samples with no à priori classification, we defined all measured samples as one group of 168 samples. To generate the second, reference (control) group, we used the chromatogram of a hexane solution of n-alkanes (C7‒C40, 1ng/µl each, Sigma Aldrich), obtained using identical chromatographic conditions as the samples. Upon importing into ChromaTOF-Tile, we multiplicated this chromatogram into 16 files with variable dilution factors (0.96‒1.03) to give rise to the reference group with identical and low coefficients of variation for all analytes.
To reduce the number of multiple thousand putative analytes initially detected by ChromaTOF-Tile, operated over the range of masses m/z 29‒600, to the most informative peaks characterizing body odour samples, we took advantage of F-ratio comparison between the two predefined classes calculated by ChromaTOF-Tile for each analyte’s peak volumes. At the same time, we optimized the signal-to-noise ratio (S/N ratio) for peak inclusion to eliminate artifact peaks. Both parameters (i.e., F-ratio and S/N ratio) were fine-tuned during the pilot study with six male participants sampled six times each (Supplementary information). Desirable separation of chemical profiles of the six participants and optimum clustering of the six samples per participant were obtained for F-ratio = 2 and S/N ratio = 50 (Extended data Fig. 5). These inclusion parameters were thus used for the main analysis. The resulting dataset was then manually curated to exclude the remaining column bleeding signals and artifacts, peak tailing-related hits, and analytes identified in the extracts of the holding sticks of the cotton swabs used for sampling.
Analytes identification
For tentative identification of individual analytes included in the main analysis, both MS spectra and 2D retention characteristics were used. Electron ionization MS spectrum (EIMS) of each compound was compared with mass spectral libraries NIST 2017, Wiley Registry™ of Mass Spectral Data (8th edition) and an in-house wax ester library58. EIMS spectra of all analytes were manually inspected, eventual erroneous library hits corrected, and new identifications proposed based on diagnostic fragments. In parallel, retention indices (RI) for both GC dimensions were retrieved from ChromaTOF-Tile output data and manually compared with NIST and other sources59,60. Tentative identification is provided only for compounds matching by both the EIMS and RI data (Supplementary Table Dataset 2).
Statistical analysis
The peak volumes for all included analytes across the 168 samples were square root transformed to reduce heteroscedasticity and related to the sum of peak volumes of each sample to correct for the differences in sample intensities. For subsequent calculations, three datasets with 84 samples each were prepared and independently processed. The first two included only right or only left axillary samples, while the third consisted of averaged values from left and right axillary sample for each analyte. Each dataset was subjected to Principal Component Analysis with autoscaling and centering, performed in Canoco for Windows v. 4.52. The factor scores for the first three principal components were plotted and used to calculate the matrix of Euclidean distances among all possible pairs of participants. These distances were then used as indicators of chemical dissimilarity among different samples. The differences between within-dyad father-partner distances and averaged distances of the corresponding father to all extra-pair partners were compared with zero expected under the null hypothesis of no influence of the pairing on chemical similarities, using Wilcoxon Signed Rank tests. Independently, we tested whether the observed within-dyad father-partner distances differ from distances obtained under random pairing of fathers with partners. To do so, we generated datasets of randomly paired father and partners (999 permutations) and compared the observed average father-partner distance with the distribution of average father-partner distances in the permuted datasets.
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
31. Schäfer, L., Schriever, V. A. & Croy, I. Human olfactory dysfunction: causes and consequences. Cell Tissue Res. 383, 569–579 (2021).
32. Kohoutová, D., Rubešová, A. & Havlíček, J. Shaving of axillary hair has only a transient effect on perceived body odor pleasantness. Behav. Ecol. Sociobiol. 66, 569–581 (2012).
33. Havlicek, J. & Lenochova, P. Environmental Effects on Human Body Odour. in Chemical Signals in Vertebrates XI (eds. Hurst, J. L., Beynon, R. J., Roberts, S. C. & D., W. T.) vol. 11 199–212 (Springer, 2008).
34. Fialová, J., Roberts, S. C. & Havlíček, J. Consumption of garlic positively affects hedonic perception of axillary body odour. Appetite 97, 8–15 (2016).
35. Havlíček, J., Lenochová, P., Oberzaucher E., Grammer, K. & Roberts, S. C. Does length of sampling affects quality of body odor samples? Chemosens. Percept. (2011).
36. Lenochova, P., Roberts, S. C. & Havlicek, J. Methods of Human Body Odor Sampling: The Effect of Freezing. Chem. Scences 34, 127–138 (2009).
37. Arrindell, W. A. et al. The development of a short form of the EMBU 1: Its appraisal with students in Greece, Guatemala, Hungary and Italy. Pers. Individ. Dif. 27, 613–628 (1999).
38. Šaffa, G., Štěrbová, Z. & Prokop, P. Parental Investment Is Biased toward Children Named for Their Fathers. Hum. Nat. 32, 387–405 (2021).
39. Fraley, R. C., Heffernan, M. E., Vicary, A. M. & Brumbaugh, C. C. The Experiences in Close Relationships-Relationship Structures Questionnaire: A Method for Assessing Attachment Orientations Across Relationships. Psychol. Assess. 23, 615–625 (2011).
40. Bürkner, P. C. brms: An R package for Bayesian multilevel models using Stan. J. Stat. Softw. 80, 1–28 (2017).
41. Glenn, T. C. et al. Adapterama I: universal stubs and primers for 384 unique dual-indexed or 147,456 combinatorially-indexed Illumina libraries (iTru & iNext). PeerJ 7, e7755 (2019).
42. Klindworth, A. et al. Evaluation of general 16S ribosomal RNA gene PCR primers for classical and next-generation sequencing-based diversity studies. Nucleic Acids Res. 41, 1–11 (2013).
43. Jiang, H., Lei, R., Ding, S. W. & Zhu, S. Skewer: A fast and accurate adapter trimmer for next-generation sequencing paired-end reads. BMC Bioinformatics 15, 1–12 (2014).
44. Callahan, B. J. et al. DADA2: High-resolution sample inference from Illumina amplicon data. Nat. Methods 13, 581–583 (2016).
45. Edgar, R. C., Haas, B. J., Clemente, J. C., Quince, C. & Knight, R. UCHIME improves sensitivity and speed of chimera detection. Bioinformatics 27, 2194–2200 (2011).
46. Wang, Q., Garrity, G. M., Tiedje, J. M. & Cole, J. R. Naïve Bayesian classifier for rapid assignment of rRNA sequences into the new bacterial taxonomy. Appl. Environ. Microbiol. 73, 5261–5267 (2007).
47. Quast, C. et al. The SILVA ribosomal RNA gene database project: Improved data processing and web-based tools. Nucleic Acids Res. 41, 590–596 (2013).
48. Pafčo, B. et al. Metabarcoding analysis of strongylid nematode diversity in two sympatric primate species. Sci. Rep. 8, 1–11 (2018).
49. Douglas, Gavin, M. et al. PICRUSt2: An improved and extensible approach for metagenome inference. bioRxiv 672295 (2020).
50. Ye, Y. & Doak, T. G. A parsimony approach to biological pathway reconstruction/inference for genomes and metagenomes. PLoS Comput. Biol. 5, 1–8 (2009).
51. Brooks, M. E. et al. glmmTMB balances speed and flexibility among packages for zero-inflated generalized linear mixed modeling. R J. 9, 378–400 (2017).
52. Storey, J. D. & Tibshirani, R. Statistical significance for genomewide studies. Proc. Natl. Acad. Sci. U. S. A. 100, 9440–9445 (2003).
53. Gotelli, N. J. & McCabe, D. J. Species co-occurrence: A meta-analysis of J. M. Diamond’s assembly rules model. Ecology 83, 2091–2096 (2002).
54. Kery, M. & Royle, J. A. Applied Hierarchical Modeling in Ecology: Analysis of distribution, abundance and species richness in R and BUGS: Volume 2: Dynamic and Advanced Models. (Academic Press, 2020).
55. R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing. (2020).
56. Marney, L. C. et al. Tile-based Fisher-ratio software for improved feature selection analysis of comprehensive two-dimensional gas chromatography-time-of-flight mass spectrometry data. Talanta 115, 887–895 (2013).
57. Parsons, B. A. et al. Tile-Based Fisher Ratio Analysis of Comprehensive Two-Dimensional Gas Chromatography Time-of-Flight Mass Spectrometry (GC × GC-TOFMS) Data Using a Null Distribution Approach. Anal. Chem. 87, 3812–3819 (2015).
58. Urbanová, K., Vrkoslav, V., Valterová, I., Háková, M. & Cvačka, J. Structural characterization of wax esters by electron ionization mass spectrometry. J. Lipid Res. 53, 204–213 (2012).
59. Adams, R. P. Identification of Essential Oil Components by Gas Chromatography/mass Spectrometry. (Allured Publishing Corporation, 2007).
60. Stránský, K., Zarevúcka, M., Valterová, I. & Wimmer, Z. Gas chromatographic retention data of wax esters. J. Chromatogr. A 1128, 208–219 (2006).