The study cohort consisted of a subset of hospitalized subjects enrolled in our Clinical Humoral Biological Sample Library. We collected 128 cases that, according to their CT scanning characteristics, were suspected lung cancer (LC) cases, and excluded the inappropriate patients in light of our clinical research design (#2021LC2115). Details of inclusion and exclusion criteria and workflow are displayed in Table 1 ; Figure S1 . A final diagnosis of LC depended on pathological characteristics of tissue samples from electronic bronchoscopy-mediated needle aspiration biopsy after BALF collection. At enrollment, we included patients with lung cancer who had not been treated with pharmacological interventions for the previous 3 months, such as anti-tumor regimes, antibiotics, probiotics intake, and other potential preparations that might affect local and extensive microbial compositions. Exclusion criteria included patients with concomitant infectious or inflammatory respiratory diseases, tumor-associated obstructive pneumonia, and patients using glucocorticoid drugs in the preceding 6 months. All patients fully understood the objectives and were volunteers for potential inspection risks. Each subject signed an informed consent approved by the Ethics Committee of the First Affiliated Hospital of the Air Force Medical University; the Academic Integrity Supervision Committee of Air Force Military Medical University carried out supervision of the whole course within the study.

Samples processed for microbiota analysis were collected from patients consulting for medical assistance in our center who needed electronic bronchoscopy-mediated needle aspiration biopsy to reach a definite diagnosis. Before that, bronchial and alveolar lavage fluid was obtained from normal lung segment (NLS) and tumor-burden lung segments (TBLS) successively within the same lung lobe. Each lavage was treated with preheated sterile physiological saline for 50-60ml, maintaining a stable recovery rate of >60%. All samples intended for microbial analysis were under centrifugation at 4°C 12000rpm for 40 min. Centrifugal sedimentation and supernatant were segregated and restored at -80°C for microbial and targeted metabolomics analysis concurrently until processing. All processes strictly abided by sterile operating standards.

BALF precipitation samples (1-3mg) were weighed in 2 ml microcentrifuge tubes and placed on ice. Total DNA from the lower respiratory tract microbiotas was extracted using the QIAamp Fast DNA Stool Mini Kit (QIAGEN, Germany) per the manufacturer’s instructions (see the QIAamp Fast DNA Stool Mini Kit Handbook, www.qiagen.com/handbooks). The degradation degree and potential contamination of the DNA were analyzed using 1% agarose gels. The DNA purity was determined using the NanoPhotometer^® spectrophotometer (IMPLEN, CA, USA). DNA samples were further diluted with sterile water to an OD value between 1.8 and 2.0, measuring with the Qubit^® dsDNA Assay Kit in Qubit^® 2.0 Fluorometer (Life Technologies, CA, USA). One microgram of qualified DNA was used to construct the library via NEBNext^® Ultra DNA Library Prep Kit for Illumina (NEB, USA). DNA samples were fragmented to 350 bp by sonication, and then the DNA fragments were end-polished, A-tailed, and ligated with the full-length adaptor for Illumina sequencing with further PCR amplification. Libraries were analyzed for size distribution using the Agilent2100 Bioanalyzer (Agilent, USA) and quantified via real-time PCR to keep size distribution of DNA fragments >3nM. The libraries were then sequenced on an Illumina PE150 HiSeq platform.

Raw data obtained from the Illumina PE150 sequencing platform were preprocessed by Readfq (V8, https://github.com/cjfields/readfq) to obtain clean data for subsequent analysis. The clean data were utilized for assembly analysis with MEGAHIT software (v1.0.4-beta in a –presets meta-large (–end-to-end, –sensitive, -I 200, -X 400) parameter settings, and the Scaftigs were obtained by breaking the resulted scaffolds from the N junction. All the sample details on the quality of their assemblies are present in Table S1 .

The Scaftigs (≥ 500 bp) were submitted to predict the open reading frame (ORF) using MetaGeneMark (V2.10; http://topaz.gatech.edu/GeneMark/) to filter out the excessive information with a length less than 100nt, and CD-HIT software (V4.5.8; http://www.bioinformatics.org/cd-hit/) to eliminate redundancy. Clean data of each sample was aligned to the initial gene catalog by using Bowtie2 (V2.2.4; https://bowtie-bio.sourceforge.net/bowtie2/) to calculate the number of reads of the genes on each sample alignment, with parameter settings: –end-to-end, –sensitive, -I 200, -x 400. Genes with reads ≤2 in each sample were filtered out to finally determine the gene catalog (Unigenes) for subsequent analysis ( Tables S2, S3 ). Based on the number of reads aligned and the length of the gene, the abundance of each gene in each sample was calculated by the following formula:

in which r is the number of gene reads on alignment, and L is the length of the gene (Qin et al., 2010). Based on the abundance of each gene in the gene catalog in each sample, basic information statistics, core-pan gene analysis, correlation analysis between samples, and Venn diagram analysis of gene number were performed.

The obtained unigenes were used to blast the sequences for the bacteria, fungi, archaea, and viruses, which were extracted from the NR database (V20180102; https://www.ncbi.nlm.nih.gov/) of NCBI using DIAMOND software (V0.9.9.110; https://github.com/bbuchfink/diamond/). We used the lowest common ancestor (LCA) algorithm to obtain the number of genes and abundance information for each sample in each taxonomic hierarchy (kingdom, phylum, class, order, family, genus, and species). DIAMOND software was also used to blast unigenes to functional databases, including the KEGG (V20180101; http://www.kegg.jp/kegg/) databases, for the blast results, and the best blast hit was used for subsequent analysis.

According to the alignment results, the relative abundance at different functional levels was calculated (the relative abundance at each functional level was equal to the sum of the relative abundance of genes annotated at that functional level). The gene number table of each sample at each taxonomy level was derived from the result of functional annotation and gene abundance table. The number of genes with a certain function in a sample was equal to the number of genes whose abundance was non-zero among the genes annotated with this function. Based on the abundance table at each taxonomy level, annotated genes statistics, relative abundance overview, and abundance clustering heat map were carried out, combined with PCA and NMDS analysis of dimension reduction, ANOSIM analysis of inter-/intra-group differences based on functional abundance, metabolic pathway comparative analysis, as well as Metastat and LEfSe analysis on the inter-group functional difference.

SCFA contents in BALF supernatant were detected by Metware Biotechnology Co., Ltd. (Wuhan, China) with gas chromatography-tandem mass spectrometry analysis. Briefly, BALF samples were thawed and vortexed for 1 min prior to analysis. A total of 50μL of samples were mixed with 100μL of phosphoric acid (0.5% v/v) solution, vertexing for 3 min and ultrasonicating for 5 min. After that, the mixture was centrifuged at 12000 rpm for 10 min at a temperature of 4°C. The supernatant was collected and used for GC-MS/MS analysis. Agilent 7890B gas chromatograph coupled to a 7000D mass spectrometer with a DB-5MS column (30m length × 0.25mm inner diameter × 0.25μm film thickness; J&W Scientific, Folsom, CA) was used. Helium was used as the carrier gas, at a flow rate of 1.2mL/min. Injections were made in the spitless mode, and the injection volume was 2μL. The oven temperature was held at 90°C for 1 min, raised to 100°C at a rate of 25°C/min, raised to 150°C at a rate of 20°C/min, and held at 150°C for 0.6 min. Then, the temperature was further raised to 200°C at a rate of 25°C/min and held at 200°C for 0.5 min. After running for 3 min, all samples were analyzed in multiple reaction monitoring mode. The temperature of the injector inlet and transfer line were held at 200°C and 230°C, respectively.

The random forest algorithm was applied to elucidate the influence of candidates on lung cancer prediction by repeated cross-validation. Further analyses were carried out in R software (v3.5.2). The LASSO logistic regression model was performed to select the most useful prognostic risk factors for SCFA candidates in BALFs collected from lower respiratory tracts. All samples were identified using dummy variables. We used R software version 3.6.1 and the “glmnet” package (R Foundation for Statistical Computing, Vienna, Austria) to perform the LASSO logistic regression analysis.

The significance of the differences between groups was analyzed using the Wilcoxon rank-sum test and ANOSIM with P value< 0.05 (5% level of probability) with VEGAN of R package being considered to be significant and denoted as follows: ^* P<0.05, ^** P<0.01, and^*** P<0.001. The statistical significance was adjusted for multiple testing using FDR correction with the cutoff adjusted p-value< 0.05 unless otherwise stated. The receiver operating characteristic curve (ROC) analysis was performed using the R project, and the discriminative power of the predictor was assessed by calculating the area under the receiver operating characteristic curves (AUC). A variable with an AUC above 0.7 was considered useful. Significant differences between corresponding subgroups were determined via an unpaired t-test and a false discovery rate approach using the two-stage linear step-up procedure with a false discovery rate (Q) of 1%. Testing conditions were analyzed individually, without assuming a consistent SD. Statistical analysis was performed with GraphPad Prism (V9.0.0 for Windows; www.graphpad.com).

From May 2022 to December 2022, we collected 128 patients with highly suspected lung cancer based on computed tomography scanning (CT) with typical malignant imaging features, including solitary or multifocal mass nodular shadow, unsmooth edges with a burr, and microvascular insertion, in light of independent judgment from our Pulmonary Nodule Diagnosis and Treatment Center. Typical CT scanning and corresponding 3D view of the targeted lesion within a representative patient among this cohort was displayed as follows ( Figures 1A, B ). All subjects were evaluated to undergo lung malignant lesion biopsy after bronchoalveolar lavage in adjacent normal segments of the ipsilateral lobe and tumor-burden lung segment via electronic bronchoscope ( Figure 1C ). After the exclusion of benign lesions and other interference factors of sample acquisition, only those patients with pathological diagnoses of malignancy were successfully enrolled, with follow-up sequencing and analysis being carried out ( Figure 1D ). The demographics of the participants are shown in Table 1 and specific inclusion criteria and other exclusion criteria are displayed as a flowchart in Figure S1 . Since the samples were also taken as the self-control of the same patient, we did not set up a blank control group in this study.

To determine compositional diversity between tumor-burden lung segment (TBLS) and ipsilateral normal lung segment (NLS), we drafted those precipitations of BALF samples collected from corresponding pulmonary segments or subsegments profiling with shotgun metagenomics sequencing, generating 1.3Gbp of sequencing data on average, and further analyzed their alpha diversity indices for the subset of final enrolled samples. Consequently, multidimensional scaling (MDS), an ordination plot based on Bray-Curtis dissimilarities, revealed distinct lower respiratory tract microbial compositions among both groups at the species level (Stress=0.1311; ADONIS P ^**=0.001; ANOSIM P ^**<0.001), with the majority of TBLS samples overlapping with the NLS subjects ( Figure 2A ). Additionally, the alpha-diversity comparison of indicated groups also demonstrated low taxonomic abundance in the TBLS-BALF subgroup by Simpson index (^** P<0.001, Wilcoxon rank sum test), which had no significance in the Shannon index ( Figure 2B ). Across the board, however, the lower respiratory tract microbiome at both the phylum and genus levels rarely fluctuated no matter which samples we sequenced ( Figures S2A, B , 2E, F ). Other beta diversity analyses seemed to reach the same conclusion as mentioned above ( Figures S2C, D ). These results suggested a perspective that despite restricted loaded biomass, minor alterations in the lower respiratory tract microbiota, especially several key species, facilitated a microbiota prone to oncogenesis and tumor development.

To further explore the differences among species that presented spatially in NLS and TBLS at the time of microscopic examination, we identified 18 differentially abundant microbial species in the comparison between both groups (FDR P<0.05, Wilcoxon rank-sum test) ( Figure 2C ). Meanwhile, linear discriminant analysis effect size (LefSe) was performed to uncover the potential the tumor-related species biomarkers. We compared the microbiota compositions of the above candidates by the LDA score of the species (log10) to enlighten the distribution diagram of species differences ( Figure 2D ), finding that the relatively abundant microbial species were differentiated in TBLS and NLS ( Figure S3A ). CIRCOS plot of taxonomic abundance within each sample also verified the outcomes mentioned above ( Figure S2G ). A species co-abundance network among this differential genus between both lung segments further suggested that the high abundance of C. jejuni in TBLS might promote the dominance of Firmicutes and impede Bacillota by their intra-phylum positive associations along with the negative associations with Bacillota species ( Figure S3B ). Particularly, a Cladogram based on differential candidates also revealed that specific taxa related to lung cancer differed from those in normal lung segments, characterized by genus enrichment of Campylobacter, Enterobacter, Debaryomyces, and Fusobacterium in tumor-burden lung segments, which were replaced by Bacillus, Klebsiella, and Acinetobacter in normal lung segments ( Figure 2E and Table S4 ), indicating the consistency of pathogenic microbial genus from biological evolutional perspectives. Collectively, these results further illustrated that compositional variations existed in cancer-loaded segments, some of which were quite distinct from those in healthy lower respiratory tract. Given the transient and significantly variable nature of normal lung microbiota in a relatively open environment (Dickson et al., 2015), the presence of a specific community could signal an ongoing pathological process providing bacteria with nutrients, a process that also deserves additional attention.

Except for the direct cytotoxic effects of the majority of viruses and quite limited bacteria species, metabolites accounted for the interaction between microorganisms and hosts (Bhatt et al., 2017; Sepich-Poore et al., 2021). Short-chain fatty acids derived from the intestine are important protective lipid metabolites released by anaerobic or facultative anaerobic microbiomes to regulate distant primary tumors (Kim et al., 2016). Despite the extensive literature on the inhibitory function of gut microbiome-derived SCFAs, several lower respiratory tract microbiota at the distal end of the tumor lesion could utilize SCFAs to regulate the local ecological environment (Jin et al., 2018; Yue et al., 2020). Correspondingly, to examine the dominant SCFAs in lung cancer blockade, except for the influences from the gut microbiome, we further detected SCFAs in BALF samples mentioned above to screen out the predictive components of SCFAs in lung cancer initiation or those associated with clinical diagnosis. To our surprise, SCFAs were generally expressed at a low level in the lower respiratory tract and were slightly increased in the TBLS group but with no significance ( Figure 3A ). This outcome seemed difficult to confront in light of the probiotic effects of SCFAs in preventing tumor process, and inevitable bias or other unknown correlated noise could have contributed to the outcome. If anything, the release of SCFAs-oriented from the lower respiratory tract within different lung segments of the same lung cancer patient was prone to be identical, regardless of the differentiated microbial composition, which was in line with previous studies.

From a practical perspective, however, exploring the predictive value of a single metabolite under sophisticated circumstances in lower airways seemed unacceptable, due to the potent interactional multiplicities between the microbiome and the host. Thus, we reconstructed a machine learning-based multivariate prediction model to clarify the predictive function of SCFAs. LASSO regression coefficient profiles of the seven SCFA candidates showed that priorities for prediction were given to combined metabolites of three SCFAs, namely, CA, VA, and IBA ( Figures 3B , S4A ), which was also confirmed by the Random forest prediction model and ROC curve based on repeated cross-validation from SCFA candidates ( Figures 3C–E , S4B ). Despite restricted accuracy of under 50%, the predictive value of this combined model should be highlighted, probably because it presented a new lung cancer diagnostic approach based on metabolic exhalation detection, deserving further validation in clinical settings.

The production of SCFAs bears a tight correlation with anaerobic or facultative anaerobic microbiome in guts, supported by sufficient findings that the fluctuation of microbial metabolites may be attributed to microbiome compositional diversity (Asnicar et al., 2021). Next, we implemented an integrated analysis of the candidate microbial species and SCFAs, in order to screen out dominant SCFA-associated microbes in tumor-burden lower respiratory tract. As a consequence, CCA profiling showed that the potential correlation between SCFAs and differential microbes mattered in tumor-burden segments with merely low efficiency ( Figure 4A ), partially due to restricted abundance and sample capacity. Heatmap of microbial species and SCFAs might present explicit correlations of differential microbes and SCFAs ( Figure 4B ), indicating a positive SCFA correlation with Brachyspira hydysenteriae and a negative connection with Pseudomonas at the genus level. These results further illustrated that the microbial-metabolic prediction model facilitated cancer screening and diagnosis by bronchoscopy-dependent BALF examination, which still deserves detailed evaluation in a large-scale population. Furthermore, as to significant correlations with clinical characteristics mentioned in other studies (Ubachs et al., 2021), we found that SCFAs and differential microbes were bound up with various clinical factors ( Figures 4C, D ), including sex, smoking status, TNM stages, and tumor gradings, although these correlations might be triggered indirectly by other unverified factors. Owing to the lack of experimental verifications of indicated candidates correlated with these characteristics, additional preferences should be given in clinical studies to further demonstrate the underlying role of the lower respiratory tract microbiome.

The Bray-Curtis dissimilarities based on KEGG pathway abundances illustrated the marginally separate clusters of NLS and TBLS (ANOSIM, ^** P<0.01) ( Figure 5A ). The KEGG pathway enrichment analysis of the metagenomic data showed that activated pathways in TBLS overlapped with those in NLS, whereas minor differences were detected only in environmental information processing and metabolism-related cascades, including cellular community-prokaryotes, signaling transduction, membrane transport, metabolism of cofactors and vitamins, and carbohydrate metabolism ( Figures 5B , S5A, B ). It is reasonable to speculate that microbe-mediated host interactions were achieved by microbial metabolites, which might induce oncogenesis or other tumor processing in a complicated microenvironment in lower respiratory tracts, further validated by a restricted proportion of functional cascades based on the KEGG pathways ( Figures 5C , S5C ). These results were also in accordance with various previous works that showed that the utilization of complex metabolites to induce local chronic inflammatory stimulation may be one of the dominant factors in microbial-mediated tumor development and progression (Hosseinkhani et al., 2021).