If not mentioned otherwise, all chemicals were purchased at the highest purity available from Sigma-Aldrich (Taufkirchen, Germany) or Carl Roth (Karlsruhe, Germany).

We utilized Melanoma 3D skin models (MLNM-FT-A375-AFAB) procured from MatTek (Ashland, MA, USA). These were maintained in six-well plates (Techno Plastic Products AG, Trasadingen, Switzerland) containing 2.5 mL of antibiotic-free maintenance medium (MLNM-FT-MM-AFAB) at 37 °C and 5% CO₂. Following 2 days of recovery with daily medium exchange upon arrival, the models were inoculated with bacteria collected from the volar forearm from a healthy volunteer using sterile omni swabs (QIAGEN GmbH, Hilden, Germany) as described previously (Sowada et al., 2014). The volar forearm exhibits the highest diversity of bacterial species (Grice et al., 2009) and was consequently selected for the collection. The volunteer gave written consent according to German federal government standards. The study has been approved by the Institutional Review Board of the German Federal Institute for Risk Assessment (BfR) with the SFP number 1323–107 and the internal procedural reference number BfR-182.

Briefly, an area of 1 cm² skin/melanoma model was sampled and the swabs were incubated in 1 mL phosphate-buffered saline (PBS; Ashland, MA, USA) for 10 min at 32 °C and 600 rpm. The swabs were removed using sterile tweezers and the liquid was centrifuged. The cells were resuspended in PBS and 15 μL bacterial suspension were applied on the models, while 15 μL sterile PBS were used for the uncolonized control models. To ensure a successful colonization, melanoma models were colonized twice at intervals of 2 days, i.e., on day 0 and day 2 (Figure 1A).

All melanoma models were cultivated aerobically for a period of 12 days with culture medium being exchanged every day. Culture medium was aliquoted and stored at −80 °C until analyses. The models were harvested either 4 or 12 days post inoculation. Bacteria were removed from the models for quantification and sequencing analysis. Epidermis and dermis were separated using sterile tweezers, and then snap frozen with liquid nitrogen and stored at −80 °C. Melanoma models utilized for tissue sections were transferred into embedding medium (Sakura Finetek, Torrance, CA, USA) and rapidly frozen in liquid nitrogen.

The quantity of bacteria was determined by taking surface imprints as outlined by Lemoine et al. (2020). In brief, 2 cm² of a sterilized velvet cloth per skin model were soaked in sterile PBS and carefully applied on the surface of the skin models. The velvet was incubated in 1 mL PBS for 30 min at room temperature at 600 rpm to detach the bacteria, wrung out using tweezers and serially diluted to count colony-forming units (CFUs) on Lysogeny Broth (LB) agar plates.

Metagenomic DNA was isolated from 350 μL of both skin swabs (inocula), surface imprints from non-colonized control, and colonized melanoma models, as well as DNA extraction reagents to capture the background (sterile PBS, water, and kit control), using the QIAamp UCP DNA Micro Kit (QIAGEN GmbH, Hilden, Germany) following the manufacturer’s instructions. Nucleic acid concentrations were measured using a UV–Vis Spectrophotometer (NanoDrop-ND 1000, Peqlab, Erlangen, Germany) and Qubit (Invitrogen, Carlsbad, CA, USA).

The V1-V3 regions of the 16S rRNA genes were targeted for amplicon library preparation using primer 27F-YM (5′-AGA GTT TGA TYM TGG CTC AG-3′) and 534R (5′-ATT ACC GCG GCT GCT GG-3′) (Frank et al., 2008; Fierer et al., 2005). Briefly, after a 2-step PCR, which adds adapter and barcodes for sequencing, AMPure XP beads (Beckman Coulter, Krefeld, Germany) were used for amplicon purification. All amplicons were pooled equimolarly and 10% (vol/vol) PhiX standard (Illumina) was added. The final library was sequenced on a MiSeq system using PE300 cartridges v3 with 2 × 275 cycles (Illumina Inc., San Diego, CA, USA) following the manufacturer’s instructions (Reitmeier et al., 2020).

Raw reads of the 16S rRNA gene amplicon data were trimmed 5 bp either side (Haider et al., 2024) and further processed using BiotaWiz, with built-in the Integrated Microbial Next-generation sequencing 2 (IMNGS2) pipeline, an updated version of IMNGS (Lagkouvardos et al., 2016). IMNGS2 is based on UPARSE (Edgar, 2013) and additionally includes the Taxonomy Informed Clustering algorithm for enhanced taxonomic resolution, generating species-like operational taxonomic units (SOTUs) (Kioukis et al., 2022). Briefly, denoised zero-radius OTUs (zOTUs) were clustered within each identified genus to produce SOTUs. Only SOTUs with a relative abundance ≥ 0.25% in at least one sample were kept for further analysis (Reitmeier et al., 2021). Moreover, potential contaminants identified in negative controls and non-inoculated melanoma models were thoroughly removed. SOTU sequences were identified by EzBioCloud’s 16S rRNA gene-based ID (Yoon et al., 2017). Assessment of α-diversity based on species richness and taxonomic binning were conducted using Rhea in R version 4.2.3 (Lagkouvardos et al., 2017). For visualization of the taxonomy, taxa with a minimum abundance of 0.5% were retained and plotted using ‘ggplot2’ version 3.5.2 (Wickham, 2016). Non-metric multidimensional-scaling (NMDS) plot of β-diversity was calculated by PERMANOVA and visualized using the ‘vegan’ package version 4.4.0 in R (Dixon, 2003).

The extraction of total RNA from the epidermis of the melanoma models was carried out using a TRIzol-based protocol (Chomczynski and Sacchi, 1987). The epidermis was disrupted twice at 25 Hz for 5 min in a TissueLyser II (Qiagen, Hilden, Germany) with TRIzol™ Reagent (Invitrogen, Carlsbad, CA, USA) in accordance with instructions from the manufacturers. The RNA concentration was then measured using a UV–Vis Spectrophotometer (NanoDrop-ND 1000, Peqlab, Erlangen, Germany). RNA-integrity (RIN) was determined on an Agilent 2100 Bioanalyzer System (Agilent Technologies, Waldbronn, Germany) using the Agilent RNA 6000 Nano Kit (Agilent Technologies, Waldbronn, Germany) following the manufacturer’s instructions.

Only RNA samples, which had a RIN > 7 as well as A₂₆₀/A₂₈₀ ratios > 1.5 and A₂₆₀/A₂₃₀ ratios > 1, were used. A minimum of 500 ng RNA per sample were sent for microarray gene analysis. The samples were labeled with a terminal deoxynucleotidyl transferase (TdT) using a proprietary DNA labeling reagent at ATLAS Biolabs (Berlin, Germany), who performed the microarray analysis on Human Clariom S Assays (Applied Biosystems, Foster City, CA, USA).

For statistical data analysis and visualization of microarray data, R version 4.4.0 was used (Team RC, 2020). The raw data (CEL files) were normalized and summarized using the robust multiarray average (RMA) algorithm from the R package ‘oligo’ version 1.68.2 (Carvalho and Irizarry, 2010). The genes were associated with vendor-provided annotations and control genes were removed after background correction. Genes with low expression were removed prior to further analysis, leaving 21,250 probes corresponding to 19,358 genes (intensity > 4 in at least two samples). One biological replicate from the non-colonized control melanoma models on day 4 was removed as a likely outlier following quality checks. The remaining data were subjected to differential gene expression analysis using the R package ‘limma’ version 3.60.2 (Ritchie et al., 2015). The linear models were fitted and moderated t-statistics computed using the eBayes function, followed by Benjamini-Hochberg correction for exclusion of false positives. Differentially expressed genes were considered significant if they showed an adjusted p-value (adj. p) < 0.05 and a log2 fold change > 0.5 as well as log2 fold change < −0.5. Probabilistic principal component analysis (PCA) was applied on Pareto-scaled and scaled data using the R package ‘pcaMethods’ version 1.96.0 (Stacklies et al., 2007). Volcano plot was displayed using the R package ‘ggplot2’ version 3.5.1.

Gene ontology (GO) term enrichment analysis for biological processes was performed using a hypergeometric test by the R package ‘clusterProfiler’ version 4.12.0 (Wu et al., 2021).

For functional interpretation of differential gene expression results, gene lists were subjected to Ingenuity Pathway Analysis (IPA, version 111,725,566, Qiagen Bioinformatics, Redwood City, CA, USA). Statistical significance was assessed using the integrated Fisher’s exact test, with adjustments made through Benjamini-Hochberg correction.

The reverse transcription of 1 μg RNA was carried out using the high-capacity cDNA reverse transcription kit (Applied Biosystems, Foster City, CA, USA). Reverse transcribed genes were quantified using the fast SYBR Green mix (Applied Biosystems, Foster City, CA, USA) on a QuantStudio 3 Real-Time PCR instrument (Applied Biosystems, Foster City, CA, USA). The specific primers were obtained from Eurofins (Ebersberg, Germany) with the following sequences: CDH1 forward: 5′-AAG AAG CTG GCT GAC ATG TAC GGA-3′, CDH1 reverse: 5’-CCA CCA GCA ACG TGA TTT CTG CAT-3′, CDH2 forward: 5’-CCT CCA GAG TTT ACT GCC ATG AC-3′, CDH2 reverse: 5′-GTA GGA TCT CCG CCA CTG ATT C-3′, HPRT forward: 5′-GTT CTG TGG CCA TCT GCT TAG-3′, HPRT reverse: 5’-GCC CAA AGG GAA CTG ATA GTC-3′. Hypoxanthine–guanine phosphoribosyl transferase (HPRT) was used as reference gene. The relative transcript levels were calculated using the 2^-ΔΔCT method (Pfaffl, 2001).

Secretion of lactate dehydrogenase (LDH) was assessed using the Cytotoxicity Detection Kit (LDH; Roche, Basel, Switzerland) following the instructions from the manufactures. The absorbance was measured at wavelengths of 490 nm and 690 nm utilizing an Agilent Biotek Synergy 2 plate reader (Thermo Scientific, Schwerte, Germany).

The measurement of cytokine secretion was carried out using a customized LEGENDplex™ multi-analyte bead-based multiplex assay plate according to the instructions of the manufacturer (Biolegend, San Diego, CA, USA) with following analytes: IL-13, IL-8, IL-1a, IL-6, IL-10, IFN-y, VCAM-1, ICAM-1, PIGF, GM-CSF, TIM-3, sRAGE, VEGF. A BD FACS ARIA III (Becton Dickinson, Franklin Lakes, NJ, USA) flow cytometer was used for detection and the LEGENDplex™ data analysis software Suite (Biolegend, San Diego, CA, USA) for data analysis.

The excreted melanoma marker MIA were measured using the ELISA kit for Melanoma Inhibitory Activity Protein 1 (MIA1) (SEH650Hu) from Merck (Darmstadt, Germany) in accordance with the protocols provided by the manufactures.

Cryopreserved tissues were sectioned, mounted to super frost slides (Thermo Scientific, Waltham, MA, USA) and fixed in ice-cold methanol. After evaporization, the tissues were rehydrated in DPBS (Pan-Biotech, Aidenbach, Germany). The cryosections were permeabilized in 0.5% Triton X-100 in DPBS for 9 min, washed twice with DPBS, and blocked with 0.2% Tween and 3% bovine serum albumin (BSA) in DPBS for 40 min. Slides were incubated with primary antibody for 2 h, washed, and incubated with secondary antibody for 50 min. The following antibodies were applied and diluted in blocking solution: Cadherin-1 1:30 (33–4,000; Thermo Scientific, Waltham, MA, USA), Cadherin-2 1:30 (18–203; Abcam, Cambridge, UK), Alexa fluor 488 Goat anti-Rabbit IgG (H + L) 1:200 (A-11008, Thermo Scientific, Waltham, MA, USA) and Alexa fluor 488 Goat anti-Mouse IgG (H + L) 1:400 (A-11001, Thermo Scientific, Waltham, MA, USA).

After antibody staining, slides were washed three times using 0.5% Triton X-100 in DPBS and then once in DPBS. Sections were counterstained with 1 μg/mL Hoechst in DPBS for 10 min and mounted in Fluor Save Reagent (Sigma-Aldrich, Taufkirchen, Germany). Stained slides were examined using a LSM700 confocal microscope (Carl Zeiss, Oberkochen, Germany) with the “Tile Scan” function in the ZEN 2012 black edition software (Carl Zeiss, Oberkochen, Germany), using a 2×2 tile configuration to create composite images. Relative fluorescence intensity was analyzed using the ZEN 3.1 lite blue edition software (Carl Zeiss, Oberkochen, Germany). The whole image was marked as region of interest, for which the fluorescence intensities were extracted.

To assess cell apoptosis, tissue sections were stained with the TUNEL In Situ Cell Death Detection Kit, Fluorescein kit (Roche, Basel, Switzerland), following the manufacturer’s instructions. The stained tissues were then analyzed using the fluorescent microscope BZ-X (KEYENCE, Neu-Isenburg, Germany).

Unless otherwise specified, data visualization was conducted using the R package ‘ggplot2’ version 3.5.2, and heatmaps with clustered dendrograms were generated with ‘pheatmap’ version 1.0.13 within the R environment version 4.4.3.

Statistical tests were conducted using R programming or GraphPad Prism version 10.1.2 (Graph Pad, La Jolla, CA, USA) with significance levels of *p < 0.05, **p < 0.01, ***p < 0.001, and ****p < 0.0001. The specific details of each statistical analysis can be found in the corresponding method descriptions and figure legends.

The experimental setup graphic and the summarizing figure were graphically illustrated using BioRender.

To ensure successful colonization with bacteria, melanoma models were inoculated with skin bacteria twice at intervals of 2 days, i.e., on day 0 and day 2 (Figure 1A). The melanoma models were cultivated for 12 days, with scheduled measurements of cytotoxicity and bacterial growth on days 4 and 12 after bacterial inoculation. During the entire cultivation period, neither the control nor the colonized melanoma models exhibited significant reduction in cell viability, as evidenced by the absence of DNA fragmentation in the TUNEL assay (Figure 1B) and consistent LDH release (Figure 1C).

The colonized melanoma models exhibited viable colonization, with mean counts of 1 × 10⁶ CFUs on day 4 and 1 × 10⁷ CFUs on day 12 per skin model (Figure 2A). Concomitantly, microbial composition was monitored by 16S rRNA gene amplicon sequencing. With respect to the inocula, the α-diversity was clearly reduced on day 4 and showed a further decline by day 12 (Figure 2B). The taxonomic composition of the skin swabs was initially similar, but changed during cultivation by day 4 and day 12 (Figure 2C). Both skin swabs exhibited a diverse microbiota, encompassing at least 26 distinct genera (Figure 2D). Despite the high α-diversity, the day 0 inoculum comprised ∼50% Streptococcus, whereas the day 2 inoculum contained 25% Streptococcus. The abundance and composition changed during co-cultivation (Figure 2D). On day 4, the bacterial composition varied among replicates. The predominant genera in the high CFUs melanoma models were either Streptococcus (Replicate 1) or Micrococcus (Replicate 2 and 3), while the two melanoma models with low CFUs (Replicate 4 and 5) were comparably diverse with a significant presence of Bacillus. On day 12, bacterial composition was consistent across replicates, with Streptococcus as predominant genus, likely attributable to culture-effects.

To characterize the influence of the skin bacteria on the physiology of the melanoma models, the transcriptomic responses of the epidermis and cytokine release were analyzed.

For a deeper understanding of the microbial contribution to melanoma biology, expression of selected melanoma markers was assessed and compared to the gene signature for metastatic potential of melanoma proposed by Gerami et al. (2015). According to this study, the metastatic signature mainly relies on the downregulation of 24 relevant genes. Of these 24 genes, 11 transcripts were downregulated in our melanoma models upon colonization (Figure 5A). Of those, ID2, CXCL14, DSC1, GJA1, ROBO1, TRIM29, and CLCA2 were significantly downregulated on day 12, suggesting that the regulation toward metastasis development was more pronounced by day 12.

Transcriptome data were further analyzed using Ingenuity Pathway Analysis (IPA), which identifies biological processes from a given data set, including biological pathways, gene networks, and disease mechanisms. IPA revealed an induction of melanoma-associated pathways in colonized epidermis on day 4 and 12 (Figure 5B). Interestingly, genes from the MAPK and PI3K signaling pathways relevant for melanoma progression, were significantly higher expressed in the colonized melanoma models. Furthermore, genes related to tumor microenvironment and S100 family signaling, which both impact the invasion and survival of tumor cells, were also significantly increased. In contrast, PTEN signaling was downregulated, which is consistent with PI3K pathway induction, as PTEN negatively regulates PI3K (Tamguney and Stokoe, 2007).

In line with this, the melanoma-derived growth regulatory protein (MIA) was detected at increasing levels over the course of cultivation, reaching significance on day 12 (Figure 5C).

As the GO term enrichment analysis showed that cell–cell tight junction assembly was affected in the colonized melanoma models (Figure 3B), we quantified the expression of Cadherin-1 and Cadherin-2, which both play a role in cell adhesion. Therefore, we examined mRNA levels of the corresponding genes and visualized both proteins by immunofluorescent staining (Figure 6).

The relative expression of CDH1 was clearly extenuated upon colonization (Figure 6A). In concordance with this, intensities of protein levels of Cadherin-1 were significantly downregulated (Figures 6B,C). Concomitantly, the relative expression of CDH2 was upregulated on day 12 after bacterial inoculation (Figure 6D) with protein staining revealing a significant increase for Cadherin-2 (Figure 6E), which extended into the dermal skin layers in the colonized melanoma models (Figure 6F).

Overall, the data presented here indicate a strong impact of bacterial colonization on melanoma biology. Transcriptomic changes translate, among others, into enhanced expression and excretion of metastasis markers and a switch from Cadherin-1 to Cadherin-2.