This cross-sectional and monocentric study was approved by the Ethics Committee of the Umbria Region (CER Umbria, approval No. 3990/21) with a waiver of informed consent. Nasopharyngeal swabs from individual adults presenting at the Hospital of Perugia (Umbria, Italy) for SARS-CoV-2 molecular testing were collected. The inclusion criteria consisted of collecting only specimens from the population living in Umbria, with the exclusion of people younger than 18 years or with a history of previous SARS-CoV-2 infection. Samples were collected from January to March 2021, at a time when the predominant circulating variants of concern (VOCs) were B.1.1.7 and P.1. During this interval, the vaccine campaign was in its early phase, such that the vast majority of samples were from unvaccinated individuals.

Sterile synthetic-head, plastic-shaft swabs were used to collect specimens and placed in Universal Transport Medium (UTM, Copan, Brescia, Italy). The swabs were processed for SARS-CoV-2 testing according to standard procedures/Centers for Disease Control and Prevention (CDC) guidelines. The swabs were tested by TaqPath RT-PCR COVID-19 Kit (Thermo Fisher Scientific, Waltham, MA, USA), which targets S (spike), N (nucleocapsid) and ORF1ab (open reading frame) genes, and Allplex™ SARS-CoV-2 assay (Seegene, Seoul, Korea) which targets E (envelope), N (nucleocapsid), and RdRp (RNA-dependent-RNA-polymerase) genes.

A total of 540 nasopharyngeal swabs (NPS) were collected and stratified into 18 groups based on age (18–40, 41–60, and $>$60 years), gender (male and female), and viral load (negative, low, and high viral load), each group containing at least 30 samples. Once the specimens were tested for SARS-CoV-2 viral load, the residual material was used for DNA extraction and 16S rRNA sequencing. Briefly, DNA extraction was performed using QIAsymphony® DSP Virus/Pathogen Kit (QIAGEN, Milan, Italy) in accordance with the manufacturer’s instructions. For 16S targeted sequencing, the V4-V5 portions of the 16S rRNA gene were amplified using the primer set 515YF-926R and sequenced paired-end (Illumina MiSeq platform).

The primary study end-point was the change from baseline in both diversity indexes (alpha and beta) and abundances of taxa of the microbial signatures. To define sample size, we hypothesized a taxon proportion change from baseline ranges at least equal to 5%, using a non-parametric two-sample, two-sided T-test. In order to evaluate significant differences between means of groups of interest, a minimum sample size of each group was set as equal to 30. Mann-Whitney U statistics (test power = 95%) was used with a large effect size. It follows that each of the 18 groups of interest would include at least 30 individuals in order to achieve a power of 95%.

For comparison, groups of interest were determined by age (three classes), gender (two classes), and viral load (three classes). Study end-point objectives required a differential analysis to be applied between pairs of groups with constraints on the average percentage of positive individuals. Assuming a positivity frequency of 10% (at the time of study design) with a small effect size, the minimum sample size for each of the 6 age-gender combinations was equal to 663, as determined by applying a binomial test to each group (test power = 95%; significance level = 5%). Power analysis was conducted with R version 3.4.4 (package pwr https://cran.r-project.org/web/packages/pwr/index.html). As a result, the minimum number of swabs to be collected would be 663 $\times$ 6 = 3978 (approximated to 4000 for practical purposes) in order to achieve a test power of 95%.

Within positives, the study required differential comparisons between high and low viral loads for each age-gender combination. Assuming an equal probability of the high and low viral load with a medium effect size, a minimum sample size of positive viral load within each of the 6 age and gender combinations was equal to 49, as computed by applying a binomial test to each group (test power = 95%). Thus, the minimum number of positive individuals with high and low viral load was 49. This minimum sample size was already warranted by previously assumed event probabilities with a confidence level equal to 95% in that a binomial test with 49 positives out of 663 individuals provided an estimated probability equal to 7.4% with a 95% CI = (5.5%, 9.6%), and thus rejecting the null hypothesis of a 10% positives in the population. This means that sampling 663 individuals in a population with 10% positives will guarantee more than 49 positive persons with a confidence at least equal to 95%.

The 16S double-strand target gene sequencing was performed on Illumina Miseq V3 at the Bio-Fab Research company. The sequenced amplicon relates to the V4–V5 regions (515F-Y: GTGYCAGCMGCCGCGGTAA; 926R-jed: CCGYCAATTYMTTTRAGTTT) and obtained raw fastq data was properly filtered. In particular, reads with less than 100 nucleotides (nt) were removed, and primer sequences were detected and then clipped. Thus, sequenced data were imported into the next-generation microbiome bioinformatics platform Qiime2 [19] in a genomic cloud-computing environment based on Refs. [20, 21]. Data were first denoised, dereplicated, and filtered by both phiX reads and chimera (consensus) by using the q2-dada2 quality control method [22] for detecting and correcting (wherever possible) Illumina amplicon sequence data. Sequence variants (SVs) were obtained by considering the specific sequence error profiles truncating bases at the first instance of a quality score less than, or equal to, 2 in the reads. Strands with a number of incorrect bases larger than 2% were discarded, and only reads with a minimum overlap of 20 nt were retained and joined. SVs were assigned taxonomy using a Naive Bayes classifier model trained on the Silva138 99% database trimmed to the V4–V5 region of the 16S. A phylogenetic tree was constructed via sequence alignment with Multiple Alignment using Fast Fourier Transform (MAFFT) [23] by filtering the alignment and applying FastTree [24] to generate the tree. Rarefaction curves of the Shannon index emphasize a good sequencing quality as the richness index does not increase significantly with the sampling depth for each sample. Thus, the rarefaction strategy (to 10${}^4$ reads) was pursued by comparing microbial composition between groups [25]. Adonis and permutational multivariate ANOVA analyses were performed to evaluate differences in microbial profiles (beta-diversity) using non-parametric permutational (9999 permutations) multivariate analysis of variance with a p-value of 0.05.

The alpha-diversity of each sample was assessed using Shannon and Chao1 indexes. Corresponding statistical significances in sample group comparison were determined using a Kruskal-Wallis test [26]. Also, the between-samples beta-diversity was assessed based on the 16S rRNA, and estimates were calculated on the SVs using Jaccard and Bray-Curtis distances [27]. Permutational multivariate analysis of variance (MANOVA) [28] with 999 permutations was used to test significant differences between sample groups based on both Jaccard and Bray-Curtis distance matrices. 16S SVs were also collapsed into Phylum, Class, Order, Family, Genus, and Species taxonomic levels, and their association with the group of interest was evaluated by using LEfSe (Linear discriminant analysis Effect Size) [29]. LEfSe was applied with default alpha values for the Anova and Wilcoxon test (0.05), and the Linear Discriminant Analysis (LDA) effect size has been evaluated by setting the absolute value of the logarithmic LDA threshold equal to 3.5. Other LEfSe parameters have been set to the default.

SCFA was measured in nasopharyngeal swab solution, inactivated with guanidinium. A 300 µL volume was transferred into a 10 mL screw cap vial (Agilent Technologies, Santa Clara, CA, USA) containing 130 mg of NaCl and 5 µL of deuterated internal standard solution (an aqueous solution containing d3 acetic acid, d5 propanoic acid, d7 butanoic acid and d9 pentanoic acid); the pH of the solution was adjusted to 3.5 with 10% perchloric acid (10 µL). The solution was placed at 65 °C for 5 min, keeping the mixture under magnetic stirring (500 rpm); head space-Solid Phase MicroExtraction (HS-SPME), using a PDMS/CAR/DVB 2 cm fiber (Supelco, Sigma Aldrich Merck), was employed to extract SCFA (30 min fiber exposure at 65 °C, under magnetic stirrer agitation) and transfer the analytes in the gas chromatography (GC) injection port. Before each analysis, the fiber was conditioned in the GC injector at 260 °C for 15 min.

Samples were analyzed by gas chromatography-mass spectrometry (GC-MS) using a 6890 GC and a 5973 A MSD (Agilent Technologies, Santa Clara, CA, USA) with an electron ionization (EI) source. A Zebron ZB-waxplus GC column (Phenomenex) was used. Splitless injection (3 min valve off) was performed, and data were acquired in SIM mode. The quadrupole analyzer was calibrated before each analytical session.

For SCFA quantitation, a linear calibration curve was recorded in a blank nasopharyngeal swab solution inactivated with guanidinium. Using six 300 µL aliquots, scalar amounts of the chemical standard aqueous solution of SCFA and 5 µL of deuterated internal standard solution were added and then processed as for the actual samples. The concentration range was the following, expressed as ng/µL concentration in the nasopharyngeal swab solution: acetic acid 4–250, propionic acid 0.1–50, isobutanoic acid 0.08–15, butanoic, isopentanoic, and pentanoic acid 0.01–5, hexanoic acid 0.25–5. The peak area ratio (PAR) for the specific ions of each analyte and its corresponding deuterated IS was measured to construct the calibration curve for each short-chain acid and calculate the sample concentration. For isobutanoic acid, d7 butanoic acid was used as internal standard and d9 pentanoic acid for isopentanoic and hexanoic acids. To a 300 µL volume of blank nasopharyngeal swab solution inactivated with guanidinium, the 5 µL volume of internal standard solution was only added to check if it was contaminated with SCFA; the measure was done in triplicate. The presence of acetic acid was detected, and its concentration was calculated by the standard addition method; it was estimated as being equal to 10.8 $\pm$ 0.2 ng/µL, and this was subtracted from the measured concentrations of acetic acid in actual samples.

The residual nasopharyngeal swabs after SARS-CoV-2 testing were used for DNA extraction and 16S rRNA sequencing. The taxonomic composition of the nasopharyngeal microbiome revealed that Firmicutes, Bacteroidota, Proteobacteria, Actinobacteriota, and Fusobacteriota were the most abundant phyla (Fig. 1), while Prevotella, Streptococcus, and Veillonella dominated at the genus level (Supplementary Fig. 1), in line with previous studies [18, 30].

To identify the individual factors that influence microbiome composition in the context of SARS-CoV-2 infection, individuals were stratified according to age, gender, and SARS-CoV-2 diagnostic test results. In the case of positivity, samples were further stratified according to the viral load, as defined in the Materials and Methods section. The analysis of alpha diversity revealed distinct effects of age, gender, and viral load on both richness and evenness of genera, as measured by Observed Operational Taxonomic Units (OTUs) and Shannon indexes, respectively (Fig. 2).

More specifically, when the effect of age was considered, a decrease in the richness of genera was observed in older individuals across almost all groups, while a parallel decrease in evenness was apparently restricted to females (Fig. 2B). The effect of the SARS-CoV-2 test result (Fig. 2A) and gender (Fig. 2C) was less evident without any apparent pattern among the different groups. Age (Supplementary Fig. 2), gender (Supplementary Fig. 3), and viral load (Fig. 3) also influenced the compositional structure of the nasopharyngeal microbiome, as revealed by the Jaccard and Bray-Curtis indexes.

As for alpha diversity, the compositional structure at the genus levels was mostly affected by SARS-CoV-2 viral load across all groups (Fig. 3), followed by age (Supplementary Fig. 2). The effect of gender appeared to be restricted to the high viral load across all age groups (Supplementary Fig. 3).

With regards to SARS-CoV-2 viral load, it should be noted that, in each combination of age class and gender, with the exception of the elderly male, the Jaccard and Bray-Curtis indexes are not equal for at least one viral load comparison. It follows that each population is characterized by a microbial composition capable of identifying differences among viral loads, thus justifying an associative analysis of microbial taxa to each population. It is worth mentioning, however, that the absence of significant differences in beta indexes, as observed in elderly males, does not imply the absence of significant associations and, therefore, the presence of a specific microbial composition in the same way as the other age class-gender groups.

In order to evaluate which variables contribute most and significantly to the structure of the microbial community as an indicator of viral load, the Adonis test was applied to the Jaccard and Bray-Curtis diversities using a multifactorial model that also considers the interactions between gender, age, and viral load in a sequential manner. All the variables were found to interact significantly with each other for both metrics (Table 1) since the p-value of the F-statistic (Pr($>$F)) is 0.05 for all the considered interactions. Moreover, the R${}^2$ parameter emphasizes that age and viral load are the two terms most contributing to explain the distances variability, while gender gives a small although significant contribution.

All in all, the sequencing data allows us to identify a microbial compositional signature at the genus level that distinguishes the viral load in each combination of age class and gender.

The Principal Coordinate Analysis (PCoA) plots of the beta indexes shown in Fig. 4 help in visualizing the Adonis results. The figure shows the distribution of the PCoA coordinates of the BrayCurtis and Jaccard distances evaluated over the 18 groups and organized for each of the six populations given by the combination of gender and age classes. Coordinates of each sample (dots) of the population considered and grouped by degree of infectious load, together with their distance (dashed lines) from the corresponding centroid (diamond), allow us to visualize which of the main two components explains, on average, differences in infectious load within each population. The axis labels also report the significances among viral loads explained by the corresponding PCoA component (the corresponding pairwise comparisons are shown in Supplementary Table 1). This analysis highlights that, for each of the six populations, at least one metric has a component that significantly distinguishes viral load, thus confirming that age and gender modulate, on average, microbial composition for both genera types and proportions.

Collectively, these results indicate that age, gender, and SARS-CoV-2 test result variably affect nasopharyngeal microbial composition arguing for the identification of individual microbial signatures that could predict the risk of infection.

To identify individual microbial signatures associated with SARS-CoV-2 infection, we performed high-dimensional class comparisons using linear discriminant analysis of effect size (LEfSe) [29] at the genus level. Because significant differences as a function of SARS-CoV-2 test results were detected in females across all age groups, we focused our analysis on females only. On applying the LEfSe differential abundance test, we could identify a common signature in females across all age groups (Fig. 5 and Supplementary Fig. 4). Specifically, we observed that patients positive for SARS-CoV-2 infection—and irrespective of their viral load—had a relative predominance of Lachnospirales-Lachnospiraceae over negative individuals.

Collectively, these results suggest that the expansion of Lachnospirales-Lachnospiraceae may represent a microbial signature associated with SARS-CoV-2 infection in females and the consistent presence across all age groups—each characterized by a different microbiome composition—may indicate a causative role in the process.

Lachnospiraceae have been implicated in the production of metabolites able to affect human health and disease, such as short-chain fatty acids (SCFAs) [31], although with differences at the strain level [32]. We therefore performed targeted metabolomics for SCFAs in selected samples and compared their levels between positive (low and high viral loads) and negative females. While the most abundant SCFAs with shorter chains, such as acetic, propanoic, and butanoic acids, did not differ between the two groups, the branched SCFA isobutanoic acid, as well as the SCFAs with longer chains—i.e., valeric and hexanoic acids—were significantly reduced in the SARS-CoV-2-positive females (Fig. 6).

Altogether, the targeted metabolomics analysis of SCFAs showed that the expansion of Lachnospirales-Lachnospiraceae in SARS-CoV-2-positive females is associated with a reduction in selected SCFAs, thus providing a combined compositional and functional microbial profile that characterizes females at infection.