We retrospectively analyzed 158 patients who were hospitalized in the First Affiliated Hospital of Guangzhou Medical University from March 2021 to March 2023 and showed imaging manifestations of diffuse parenchymal lung lesions. Of them, 21 were finally diagnosed with lung cancer, 93 with lung infection, and 44 with other conditions (except malignancy and infections). The clinical characteristics of the patients and metagenomic sequencing results of BALF were summarized and analyzed.

This study was registered with the Chinese Clinical Trial Register (www.chictr.org.cn) under registration number ChiCTR-CCC-12002950. All enrolled patients signed informed consent. Figure 1 presents the study design and the full inclusion and exclusion criteria.

A routine clinical status assessment was performed before bronchoalveolar lavage with strict control of indications and contraindications. Briefly, surface anesthesia was performed using nebulized inhaled lidocaine, followed by the injection of 1–2 ml 2% lidocaine into the lavage lung segment through the biopsy hole for local anesthesia. Then, 60 ml room-temperature sterilized saline was rapidly injected through the operated orifices in divided injections. Immediately thereafter, the BALF was collected by suction under negative pressure (<100 mmHg), and the total recovery rate was ≥30%. The collected BALF was divided into two sterile sealed containers for pathogenetic and cytological analysis, respectively (Society, C. R, 2017). Each sample was saved in a 2-6°C cooler, and immediately transported to a -86°C ultra-low-temperature refrigerator for storage until it was thawed before sequencing.

Nucleic acids were extracted from the BALF, and the DNA library was prepared using the PathoLib™ Genomic Library Construction Kit from Willingmed Technology (Beijing Co., Ltd). The concentrations of nucleic acid, library, and fragment length were measured simultaneously. The library was diluted to the appropriate concentration for sequencing, denatured by adding NaOH to a final concentration of 0.1 nM, mixed with Illumina HT1 solution, and then sequenced on an Illumina NextSeq 550 System (Eisenstein, 2015). Negative controls, positive controls, and internal controls were set up during the aforementioned steps. All reagents were provided by Willingmed Technology (Beijing Co., Ltd).

Relative quantitative PCR assay using the SYBR method. Amplification/oligonucleotide primer pairs were designed for Lactobacillus acidophilus based on the characteristics in the nucleotide sequences of the 16S-23S rRNA spacer regions (Forward primer 5’-TCTAAGGAAGCGAAGGAT-3’, Reverse primer 5’-CTCTTCTCGGTCGCTCTA-3’) (Tilsala-Timisjärvi and Alatossava, 1997). 18sRNA was used as an internal reference gene (Forward primer 5’-GCAATTATTCCCCATGAACG-3’, Reverse primer 5’-GGCCTCACTAAACCATCC-3’). PCR reaction conditions: 95°C for 1 minute; 95°C for 20 seconds, 56°C for 1 minute, for a total of 40 cycles. DNA levels were calculated using the 2-^ΔΔCT method.

FastQC v0.12.1 (de Sena Brandine and Smith, 2021) was used for quality assessment, and fastp v0.23.4 for quality control of sequencing data. Trimmomatic (Bolger et al., 2014) was used to remove primers, connectors, and low-quality sequences. Host sequences were removed using Kneaddata v0.12.0 (http://huttenhower.sph.harvard.edu/kneaddata) and Bowtie2 v2.5.1 (Langmead and Salzberg, 2012), and the resultant non-host sequences underwent downstream analysis. MetaPhlAn2 (Truong et al., 2015) was used to analyze the composition of the microbiome community and HUMAnN 2 (http://www.huttenhower.org/humann2) to obtain the species abundance and metabolic pathway function information of the microbiome. The above analysis was performed using the EasyMetagenome pipeline (https://github.com/YongxinLiu/EasyMetagenome). The cooccurrence network of the microbiota was generated using Cytoscape v3.7.0 and visualized in a circular layout (Sam Ma et al., 2015).

Data were shown as the mean (standard deviation) or the median (interquartile range) for continuous variables and number (%) for categorical variables. Continuous variables were compared between groups by the one-way ANOVA or Kruskal–Wallis test, and categorical variables were analyzed using the Pearson Chi-square test or Fisher exact test. A two-sided P value< 0.05 was considered statistically significant. Statistical analyses were performed using SPSS Version 22 (IBM, Corp., Armonk, NY, USA) statistical software. Box plots were drawn using GraphPad Prism v9.5 (GraphPad Software, Inc., La Jolla, CA, USA). Principal co-ordinates analysis (PCoA) Plot-Adonis were drawn for the three groups using the Bray-Curtis distance. MedCalc Statistical Software v19.0.4(MedCalc Software bvba, Ostend, Belgium; https://www.medcalc.org; 2019) was used to analyze and compare the receiver operating characteristic (ROC) curves. A microbiome cooccurrence network was established using SparCC (Sparse Correlations for Compositional data) (von Mering et al., 2012) and visualized using Cytoscape software (v3.9.0, National Institute of General Medical Sciences, Bethesda, MD, USA) (Shannon et al., 2003).

The study cohort comprised 158 patients (average age 62.76 years; 98 males and 60 females). The final diagnoses included 21 cases of lung cancer, 93 cases of lung infection, and 44 cases where the diagnosis was neither malignancy nor infection (others). Of the 21 patients with lung cancer, 16 had adenocarcinomas, 4 had squamous cell carcinomas, and 1 had small cell carcinoma. The clinical staging of these patients was as follows: 5 patients at stage I, 8 patients at stage II, 7 patients at stage III, and 1 patient at stage IV. Body mass index (BMI) was significantly lower in patients with lung cancer compared with patients with lung infection and other diagnoses, and peripheral blood leukocyte counts, and neutrophil percentages were significantly higher in patients with lung infection compared with the other two groups. The characteristics of the study participants and all P values are presented in Tables 1 and 2 .

Metagenomic sequencing of BALF from the 158 patients yielded 310,149,065 sequences, of which 287,595,289 (92.73%) were classified. The average number of sequences per sample was 1,962,968, including 82.13% human, 3.9% microbial, 4.7% low quality, 2.0% duplicates, and 7.27% unclassified.

Taxonomic analyses of relative abundance were conducted at the kingdom, phylum, genus, and species levels. At the kingdom level, Bacteria were dominant in all study subjects. Eukaryota were present in a smaller proportion in the lung infection and the others groups, but almost absent in the lung cancer group ( Figure 2A ). At the phylum level, Firmicutes, Actinobacteria, Fusobacteria, Bacteroidetes, and Proteobacteria were present in high abundance in all subjects ( Figure 2B ). At the genus level, the lung cancer group was dominated by Lactobacillus and Streptococcus (both belonging to Firmicutes), while the abundance of all other genera was lower. In the lung infection group, the genera with higher abundance were Lactobacillus, Streptococcus, Corynebacterium, Pseudomonas, Staphylococcus, Veillonella, Neisseria, Acinetobacter, and Klebsiella. In the others group, Actinomyces, Prevotella, Fusobacterium, Leptotrichia, Corynebacterium, and Rothia were in high abundance, and the proportion of each genus was relatively balanced ( Figure 2C ). At the species level, the lung cancer group was dominated by Lactobacillus acidophilus and Streptococcus mitis. In the lung infection group, common pathogenic bacteria included Pseudomonas aeruginosa, S. oralis, Acinetobacter baumannii, Staphylococcus aureus, Corynebacterium accolens, Enterococcus faecium, Candida albicans, and S. salivarius. In contrast, the composition of microbial species in the others group maintained its genus-level characteristics and was more balanced, with a relatively high proportion of F. nucleatum ( Figure 2D ). Thus, the component ratio of LRT microbiota taxa in patients with lung cancer was imbalanced and tended to be homogeneous. Table 3 shows the relative abundance of microbial taxa. Heat maps ( Supplementary Figure 1 ) were presented to demonstrate the microbiome composition of each sample.

We identified 341 types of bacteria at the species level (lung cancer group, 162; lung infection group, 318; and others group, 173). Of them, 168 (49.27%) distinct species were identified in the lung infection group, 3 (0.88%) in the others group, and none in the lung cancer group ( Figure 3A ).

Next, we compared the diversity across different groups of BALF microbiomes. Alpha diversity was lower in the lung cancer group and lung infection group at the genus level (Shannon, Simpson, Choa1, and richness index: all P< 0.001; Figures 3B–E ), detailed data are presented in Table 4 . The results of Principal co-ordinates analysis (PCoA) showed significant differences in the composition of the bacterial flora among the three groups ( Figure 3F ). These findings indicated a decreased number of species and diminished microbial diversity in the LRT microbiome of patients with lung cancer.

We searched for LRT microbial markers in patients with lung cancer by analyzing differences among all phyla, genera, and species in terms of relative abundance using the linear discriminate analysis effect size (LEfSe). The results showed that specific taxa were enriched in each group ( Figures 4A–C ). Firmicutes, Lactobacillus, and L. acidophilus were significantly enriched in the lung cancer group. In the lung infection group, P. aeruginosa, C. striatum, S. oralis, S.epidermidis, C. albicans, E. faecium, and other common pathogenic bacteria were significantly enriched. In the others group, Actinomyces, Prevotella, Fusobacterium, Leptotrichia, Rothia, and so on were significantly enriched.

The results of the relative quantitative PCR assay showed that the relative DNA levels of L.acidophilus in the BALF of patients with lung cancer were significantly higher than those in the other two groups, as shown in Figure 5A . The DNA levels of each group and the p-values for two-by-two comparisons are shown in Table 5 . This indicates the enrichment of microbial biomarkers in the LRT of patients with lung cancer. Table 3 shows the exact values of the abundance of microbial biomarkers.

Our dataset was used to establish a diagnostic model for predicting lung cancer risk. The abundance of L. acidophilus in BALF and BMI were used as candidate indices for lung cancer risk prediction and were included in multivariate logistic regression analyses to establish prediction equations (P< 0.05, respectively; Figure 5B ). Table 6 presents detailed information on the equations.

The predictive performance of the model was evaluated according to the area under the curve (AUC). The ROC value was 0.845 (P< 0.001) for BMI and 0.914 (P< 0.001) for L. acidophilus %. The combination of the two indices showed a higher ROC value than each one alone (AUC = 0.965, accuracy = 93.04%, sensitivity = 95.24%, and specificity = 92.70%; P< 0.001; Figure 5C ). We also found a statistically significant difference in AUC between L. acidophilus % or BMI alone compared with their combination (P< 0.001, respectively). Table 7 shows the ROC values.

We used data on the BMI and L. acidophilus % in BALF of three typical cases showing Diffuse lung parenchymal lesions on CT imaging who had final diagnoses of lung cancer ( Figure 6A ), lung infection ( Figure 6B ), and congenital lung malformation ( Figure 6C ) and input these into the above prediction equation. The results showed that the probability of lung cancer in these three patients was 98.61%, 0.10%, and 0.09%, respectively.

To explore potential microbiome coexistence and coexclusion relationships, we performed a co-occurrence network analysis. We selected the top 15 most abundant genera in each group and used them to build and estimate a network based on the relative abundances of microbiome genera using SparCC (Sparse Correlations for Compositional data) with an r threshold = 0.7 and a p threshold = 0.01. The degree of connectivity between genera was weaker in the lung cancer group. Lactobacillus, the most abundant genus, had few connections with other genera, while Streptococcus, the second most abundant genus, had connections with the other eight genera ( Figure 7A ). The degree of connectivity between genera was also weaker in the lung infection group, with Neisseria showing more connections to other genera. ( Figure 7B ). The genera showed a strong degree of connectivity in the others group, with Prevotella being mutually exclusive with the other six genera ( Figure 7C ). This indicated the weakened interaction between microbes in the LRT of patients with lung cancer or lung infections. The numbers presented in Tables 8 – 10 represent the correlation coefficient using SparCC for the top 15 genera in BALF across the three groups.