2.1 Volunteers recruitment
The study aimed to enroll Maskne patients and healthy individuals(HCs) aged 16-25 residing in Hunan Province, China, during the summer months from June 2023 to September 2023. The average low and high temperatures during the study period range approximately between 23-26°C and 33-38°C, respectively, with a relative humidity of around 70% to 85%. Participants who had been living in Hunan for several years, not undergone antibiotic treatment in the three months prior to sampling, and willing to avoid any other medications during the testing period were recruited.Participants were requested to choose a testing time point that did not overlap with their menstrual time. Written informed consent was obtained from each participant before registration. All procedures involving human participants in this study adhered to the principles of the Helsinki Declaration. The study had been approved by the Institutional Review Board of the Third Xiangya Hospital, Central South University (Protocol Number: [ 21141 ]).
Inclusion Criteria:
1Residing in Hunan Province, China.
2Patients:
Individuals diagnosed with Maskne who had completed a three-month washout period under medical supervision, during which they did not wear any facial masks and did not use any skincare or medications that influence the skin. Patients were classified as moderate or severe according to the Global Acne Grading System (GAGS) during the period of illness, regardless of their current skin condition.
HCs:
Individuals without acne or any other skin diseases, such as atopic dermatitis, allergic dermatitis, eczema, etc.
3Age between 16 and 25 years.
Exclusion Criteria:
Participation in other clinical studies in the past three months;pregnancy or lactating women;presence of any skin disease (e.g., atopic dermatitis, psoriasis, bruises, eczema, urticaria, rosacea, etc.);deformities or tattoos in the area of study or open wounds,use of medications within the past three months, including antibiotics, hormonal medications, retinoids, adapalene, alpha hydroxy acids, beta hydroxy acids, salicylic acid, tretinoin, or other skincare products with known efficacy in treating acne or skin were excluded.
2.2 Samples collection:
Within one week, participants should avoid swimming in chlorinated pools, as well as using hot water, saunas, or ultraviolet.Considering the time required for microbiota to regain balance after facial cleaning,patients refrained from any facial cleansing for 24 hours before the study beginning.For 12 hours before the study beginning, participants should avoid touching their faces with hands, tissues, towels, or other items.Volunteers were not allowed to wear any type of facial mask during the washout period three months before the study.
Using a sterile skin swab moistened with physiological saline, random samples were collected from the cheek area within the mask-covered region. The swab sample covered an area of 2.5 x 2.5 cm².During collection, the skin should be swabbed alternately in a cross pattern about 25 times. Samples should be collected before and after wearing the mask for 4 hours .Swabs were placed in phosphate-buffered saline and immediately stored at -20°C before DNA extraction.
In order to express conveniently, we abbreviate normal people before and after wearing masks to NB and NA, and abbreviate Maskne patients before and after wearing masks to B and A in the following paper.
2.3 DNA extraction, PCR amplification, and sequencing
F medial primer:5'-TTCCCTACACGACGCTCTTCCGATCT-barcodeF1Specific primers-3'
F Lateral primers:5'-AATGATACGGCGACCACCGAGATCTACAC-
TCTTTCCCTACACGACGCTC -3'
R medial primer:5'-GAGTTCCTTGGCACCCGAGAATTCCA- barcodeR1Specific primers-3'
R Lateral primers:5'-CAAGCAGAAGACGGCATACGAGAT- barcodeR2 -
GTGACTGGAGTTCCTTGGCACCCGAGA-3'
16S V3-V4 specific segment :
357F 5'-ACTCCTACGGRAGGCAGCAG-3’
806R 5'-GGACTACHVGGGTWTCTAAT-3'
The PCR amplification products were add to 2% agarose gel for electrophoresis and subsequently excised and collected from the gel. AxyPrepDNA Gel Collection Kit from AXYGEN was used for the recovery process. Following recovery, Illumina high-throughput sequencing and bioinformatics analysis were carried out.
2.4 Data processing
The sequencing work was conducted by Tinygene Biotech Co., Ltd. (Shanghai, China).The obtained paired-end reads (PEreads) were initially sorted based on barcodes for each sample. Subsequently, sequences quality control and filtering were performed. The sequences were then merged based on their overlap relationships,then the quality control and filtering were performed again. The specific steps are outlined below:
Sequence quality control using Trimmomatic (version: 0.38): If the average quality score within the window dropped below 20, the bases starting from the window were trimmed off the end, and reads below 50 bp after filtering were discarded.A sliding window approach was employed, using a window size of 50 base pairs (bp).
Adapter and primer trimming using cutadapt (version: 1.16): cutadapt software was used to process sequencing adapters and primers.
Sequence Merging using FLASH (version: 1.2.11): PE reads were merged into a single sequence using FLASH, taking into account the overlap relationship between PEreads. The minimum overlap length was set to 10 bp, and the maximum allowed mismatch rate in the overlap region was 0.2. Additional parameters were set as follows: maxambig=0, maxhomop=8, minlength=200, maxlength=485.Merged sequences that didn't meet the criteria were filtered out.
The resulting optimized sequences were subjected to Operational Taxonomic Units (OTUs) clustering analysis and taxonomic classification.Using USEARCH (version: 8.1.1861) software to cluster the assembled sequences into OTU,the main steps were as follows:
Utilizing UPARSE, clustering was carried out at a 97% sequence similarity level, resulting in representative sequences for each OTU.UCHIME and the existing chimera database, golddatabase (v20110519), were used for sequence alignment to remove chimeric sequences generated during PCR amplification from the representative sequences of OTUs.The usearch_global method was employed to align all sequences back to the representative sequences of OTUs, generating abundance tables for every OTU in each sample.Following the acquisition of representative OTU sequences, mothur (classify.seqs, version: 1.39.5) software was used to annotate these sequences with species information by aligning them against a reference database. The confidence threshold of 0.6 was set for the annotations. The databases used for alignment were Silva128 (bacteria) and Silva119 (archaea).OTUs without annotation results or annotated to species not relevant to the analysis project were filtered out. For example, if the analysis project focused on bacterial 16S samples, OTUs annotated as archaea were removed. Software used for this step included USEARCH (version: 8.1.1861) and mothur (version: 1.39.5).
2.5 Statistical analysis
All statistical analyses were performed using the R V3.6.0 environment. With no special instructions, the statistical results were visualized using the ‘ggplot2’ package . Alpha diversity was measured using the function ‘diversity’ in the package “Vegan” based on a flat taxonomy table. Gini-Simpson diversity index was obtained by subtracting the value of the classical Simpson index from 1. Beta diversity was compared using principal coordinate analysis (PCoA) based on Bray-Curtis distances. Redundancy analysis (RDA) was also conducted using Vegan . Beta diversity across sample groups was compared by PERMANOVA with permutations of 999. The ANOSIM test was selected to assess significance between groups following the criteria of Wang et al. , where R > 0 and p < 0.05were considered statistically significant. Relative abundance data at the genus level (excluding unclassified) were divided by group to construct separate network graphs. The ‘cor.test’ function from R version 3.6.3 was employed to compute Spearman correlation coefficients, with p-values corrected using the Benjamini-Hochberg (BH) method.Species pairs with correlation coefficients above 0.8 and p-values below 0.05 were selected to establish correlations. The ‘ggraph version 2.0.3’ package was used for network graph visualization, and the ‘group_components’ function determined Modules within the network. Nodes and edges within the same Module were assigned the same color. The ‘centrality_degree’ calculated degrees, the ‘triangles’ calculated local_triangles, and the ‘group_edge_betweenness’ computed Cluster attributes. Visualization involved the use of ‘ggplot2 version 3.4.2’ and ‘ggsignif version 0.6.0’ to create box plots, with the wilcox.test method assessing inter-group differences.The ‘centrality_hub’ computed hub_score attributes to select the Module with the highest hub_score value as the hub_network to visualization.
In addition, biomarkers of sample groups were discovered by Linear Discriminant Analysis (LDA) Effect Size (LEfSe)3. The strategy for multi-class analysis was set one-against-all, and the threshold on the logarithmic LDA score for discriminative features was set to 3.0.
It needs to be emphasized that A and B, NA and NB, are paired sample groups