Predicting how biological communities respond to disturbance requires understanding the forces that govern their assembly. We propose using human skin piercings as a model system for studying community assembly after rapid environmental change. Local skin sterilization provides a ‘clean slate’ within the novel ecological niche created by the piercing. Stochastic assembly processes can dominate skin microbiomes due to the influence of environmental exposure on local dispersal, but deterministic processes might play a greater role within occluded skin piercings if piercing habitats impose strong selection pressures on colonizing species. Here we explore the human ear-piercing microbiome and demonstrate that community assembly is predominantly stochastic but becomes significantly more deterministic with time, producing increasingly diverse and ecologically complex communities. We also observed changes in two dominant and medically relevant antagonists (Cutibacterium acnes and Staphylococcus epidermidis), consistent with competitive exclusion induced by a transition from sebaceous to moist environments. By exploiting this common yet uniquely human practice, we show that skin piercings are not just culturally significant but also represent ecosystem engineering on the human body. The novel habitats and communities that skin piercing produce may provide general insights into biological responses to environmental disturbance with implications for both ecosystem and human health.

Human research ethics approval: Protocols for study participant recruitment, data security, sample collection, and associated procedures were approved by the McGill University Research Ethics Board Office (REB-1 #70-0617).

Sample collection: From October 2019 to March 2020, we recruited 28 individuals who were receiving earlobe piercings at Tattoo Lounge in Montreal, Quebec, Canada and received their written, informed consent to participate in the study. Following standard ear-piercing protocols, we sterilized the ear lobe skin area to be pierced with a benzalkonium chloride antiseptic towelette (Jedmon Products) immediately before piercing. We pierced earlobes using a sterilized beveled hollow needle (Ruthless/Precision) and then inserted a 5/16” surgical steel grade (316L) piercing labret stud composed of chromium, nickel, and molybdenum. Both needle and stud were dipped in a water-based lubricant jelly (Personelle, Jean Coutu) to minimize friction and cleaned off after using a cotton-tipped swab. We collected skin swab samples using the DNA/RNA Shield Collection Tube w/Swab – DX (Zymo Research), which was used to preserve nucleic acids within samples at ambient temperatures. The piercer collected samples from the earlobe to be pierced and an adjacent unsterilized part of the ear farther up the ear but still part of the earlobe skin to serve as a temporal control. Samples were collected both before and after the piercing event (defined as a three-part process that includes A) skin sterilization followed by B) skin piercing and then C) insertion of the metal stud). Study participants were then instructed to self-sample both the piercing and the adjacent skin control while wearing gloves over the following 2 weeks at specified timepoints: 12 hours, 1 day, 3 days, 1 week, and 2 weeks. Additionally, environmental controls were collected by the piercer before the piercing and by the participant at the 1- and 2-week timepoints by waving a swab in the air for 30 seconds. In total, we collected 17 samples from each participant.

DNA extraction and amplicon sequencing: We extracted DNA from swabs using the DNeasy PowerSoil kit (QIAGEN) and then purified using the OneStep PCR Inhibitor Removal kit (Zymo Research). Skin swab samples and environmental controls were processed with a DNA extraction negative control included within each batch of 24 extractions. This work was carried out in a laboratory facility designed to handle low-copy and highly degraded environmental DNA samples through mitigation of contamination factors (e.g., no exposure to PCR products, regular deep cleaning, and strict usage protocols to limited trained personnel). The V1-V3 region of the 16S rRNA gene was PCR amplified using the primers 27F (5'-AGAGTTTGATCCTGGCTCAG-3') and 518R (5'-ATTACCGCGGCTGCTGG-3'). Library preparation, quality control, and high throughput sequencing with Illumina MiSeq v2/v3 were conducted at Génome Québec and the McGill Genome Centre (Montreal, Quebec, Canada).

Data processing: Raw sequences were processed using the QIIME2 bioinformatics pipeline. Primer sequences were trimmed using cutadapt before ASVs were generated using DADA2. Contaminant ASVs were identified using environmental and DNA extraction negative controls for each sequencing batch with the prevalence-based method at a classification threshold of P* = 0.5 within decontam. The unpierced control of each individual is only experimentally valid if it exhibits no significant differences from the microbiome of the skin to be pierced prior to piercing. Thus, statistical outlier individuals were defined as having an absolute difference in ASV richness between sample and control prior to piercing that was greater than 1.5 times the interquartile range across all individuals. A total of 1,047 contaminant ASVs and two statistical outlier individuals were removed resulting in 10,915 ASVs across 392 samples with a mean sequencing depth of 27,817 reads per sample. ASVs were aligned using MAFFT and phylogenetic trees were built using FastTree 2 based on Jukes-Cantor distances. For taxonomic assignment, the 27F/518R 16S rRNA primers were used to in silico extract the target V1-V3 amplicon from the SILVA 132 database. A naïve bayes classifier was trained using scikit-learn on the extracted database and then used to taxonomically assign ASVs from domain down to species. Assignments were accepted if classification confidence was at least 0.7.

Statistical analyses: ASV counts were normalized via Total Sum Scaling (TSS), and biodiversity indices, PCoA, and PERMANOVA (999 permutations) were calculated using the R ‘phyloseq’ and ‘vegan’ packages implemented within MicrobiomeAnalyst 2.0. Data was not rarefied to maximize the amount of data analyzed and the number of participants included in the study. Alpha and beta diversities were measured using Chao1 and Bray-Curtis dissimilarity, respectively. Betadisper was calculated separately using the R ‘vegan’ package version 2.6-2 and ‘ggstatsplot’ version 0.10.0 was used for plotting within RStudio Desktop version 2022.12.0+353 and R version 4.2.2. ASV co-occurrence networks were built using Random Matrix Theory (RMT)-based Spearman’s rank correlation through the Molecular Ecological Network Analysis Pipeline (MENA) implemented within the Integrated Network Analysis Pipeline (iNAP). Data was first filtered by retaining only ASVs present in >15% of samples and then log transformed before calculation of similarity matrices allowing a single timepoint lag for time-dependent interactions. Co-occurrence networks were visualized using Cytoscape version 3.9.1 keeping only nodes with valid genus-level taxonomic assignments and edges with a P-value < 0.05. The ‘iCAMP’ R package version 1.5.12 was used to calculate pNST and infer community assembly mechanisms by phylogenetic bin-based null model analysis. Bootstrapping tests with a resampling size of 1000 were used to assess significant pairwise differences between time points. Core microbiome community taxa were classified based on a minimum of 5% relative abundance across at least 50% of all samples.