A total of 18 patients with lung adenocarcinoma were enrolled from China-Japan Friendship Hospital, and the clinical information is detailed in Table 1 . Half of them were pathologically diagnosed with AIS or MIA, and other patients were diagnosed with IAC. Formalin-fixed paraffin-embedded (FFPE) samples of tumor tissues and paired adjacent normal tissues from all 18 patients were collected for whole exome sequencing (WES). Tumor tissues and BALF were obtained from 17 patients for metagenomic next-generation sequencing (mNGS). Among these, two tumor tissues and corresponding BALF samples were sequenced for 3 patients with multiple primary nodules. One patient had only BALF sample collected for mNGS sequencing.

Lung cancer patients with different degrees of pathological infiltration underwent surgical resection. During the surgery, tumor tissue, adjacent non-cancerous tissue, and BALF samples were collected. All specimens were stored at −80°C.

We transferred 2 ml BALF into a centrifuge tube, followed by the addition of lysis reagent and nuclease. The sample was incubated in a constant-temperature metal bath at 37°C for 30 minutes, then at 65°C for 10 minutes. Each fresh tumor tissue sample, approximately the size of a soybean, was placed into a centrifuge tube, and DTT and protease were added. The tube was incubated at 56°C to digest the tissue cells, followed by the addition of lysis reagent and nuclease, with a final incubation at 65°C. Both the alveolar lavage and tissue samples were transferred to lysis tubes. After adding zirconia beads, the tubes were homogenized for 30 minutes utilizing a homogenizer. DNA was eluted with 35 μL of elution buffer, and the quality of extracted DNA was assessed on Qubit^® 4.0 fluorometer (Thermo Scientific, Waltham, MA, USA).

The construction and sequencing of the DNA library were carried out following the manufacturer’s instructions (Illumina). Extracted DNA was fragmented utilizing a Covaris ultrasonic disruptor to produce fragments of approximately 200 bp. Fragment ends were repaired and A-tailed, followed by ligation with adapters containing barcodes. Subsequently, PCR amplification was performed. The library quality was assessed on Agilent Bioanalyzer 2100 system and quantified on Qubit 4.0 fluorometer (Thermo Scientific, Waltham, MA, USA). The library was sequenced using Nextseq550DX platform in Repugene Technology Co., Ltd (Hangzhou, China).

Genomic DNA extracted from formalin-fixed paraffin-embedded samples was fragmented into 150-300 bp, with paired tumor and normal tissues. After end-repair of the fragmented DNA and addition of an A-tail, adapters are ligated to both ends of the DNA fragments to construct the DNA library. Libraries with specific indices were pooled and subjected to liquid-phase hybridization using biotin-labeled probes. Exome capture was performed using Agilent SureSelect Human All ExonV6 Kit (Agilent Technologies, Santa Clara, CA, USA). The captured libraries were then linearly amplified by PCR, followed by quality control. DNA sequencing was performed on Illumina NovaSeq 6000 platform (Illumina Inc., San Diego, CA, USA) in Repugene Technology Co., Ltd. (Hangzhou, China), generating 150-bp paired-end reads with mean coverage of 200× for tumor tissue and 100× for adjacent non-cancerous tissue.

For the raw data of WES, clean reads were obtained following data filtering and alignment to the human reference genome (GRCh38). Singlenucleotide variants (SNVs) and insertions/deletions were identified using GATK (version 4.4). Additionally, MuTect2 (version 4.1) was employed to detect somatic mutations. To identify structural variants, gene fusion detection was performed using LUMPY (version 0.2.13), and copy number variation (CNV) was detected using CNVkit (version 0.9.9) (15).

For the metagenomic data, Kraken2 (version2.0.7) was employed for unique non-human sequence alignment to annotate microbial species, and count values for all species were normalized as relative abundance. MEGAHIT (version1.2.9) was used for de novo assembly of host-filtered sequences, and assembly results were statistically summarized, including final contig sequence lengths and other assembly statistics. MetaGeneMark (version3.38) was utilized for gene prediction on contig sequences from each sample using the MetaGeneMark_v1.mod model. MMseq2 (version2-13.45111) was utilized to cluster homologous genes.

Gene abundances were quantified using Salmon (version1.10.2). For diversity analysis, the vegan package (version 2.6.4) in R was used to compute α-diversity indices, such as Shannon and Simpson indices which reflect the richness and evenness of microbial communities within each sample, and a principal coordinates analysis (PCoA) to provide insights into microbial community differences between groups. Differential abundances of microbiota between groups were analyzed using rank-sum tests. To explore correlations between differential microbiota, Pearson correlations were calculated using the psych package (version 2.3.6). The metacoder package (version 0.3.6) was utilized to visualize the evolutionary relationships of microbial species. Based on the known relationships, taxonomy sets enrichment analysis (TSEA) was employed to explore associations between the microbiota taxa we identified and specific disease. Custom R functions were developed for dimensionality reduction and differential analysis of KEGG Orthology (KO) gene matrices (16). The KEGG KO library was curated, and differential KO genes were analyzed for KEGG enrichment using custom R functions. Spearman correlations between microbiota abundances and KO genes were calculated using the psych package (version2.3.6). Spearman correlations between microbiota abundance and mutation frequencies of genes were assessed using the psych package (version2.3.6).

Mutation-related analysis was conducted using cBioPortal (17) (https://www.cbioportal.org/), a comprehensive open web platform that integrates multiple datasets and offers a range of functions including data mining, data integration, and visualization. For key mutations associated with microorganisms, databases such as CPTAC, OncoSG, and TCGA were jointly applied to explore their mutational characteristics and clinical relevance. The specific analysis steps are displayed in the workflow.

Wilcoxon rank sum test was used for microbiota statistical comparison between two groups. The fisher test was performed for mutated genes statistical comparison between two groups. The hypergeometric test was used to assess statistical significance of KO genes enrichment. All statistical analyses were carried out by R (version 4.4.1). Statistical significance was defined as an adjusted P-value < 0.05.

Compared to AIS/MIA group, IAC tumors exhibited higher oncogenic mutation frequencies, indicating a heavier tumor mutation burden ( Figure 1A ). For each group, we calculated the proportion of samples with gene mutations out of the total number of samples. Top 15 genes with the highest mutation frequencies in the IAC group was visualized, and these mutations were not detected in AIS/MIA group ( Figure 1B ). The heatmap generated from metagenomic data illustrated different distribution patterns between the groups. Several species, such as Bosea sp. AS-1, Bosea sp. PAMC 26642, Bosea sp. RAC05, Bosea vaviloviae, Microbacterium paludicola, Rhodopseudomonas palustris, and Stenotrophomonas maltophilia, displayed decreased abundance in IAC group compared with AIS/MIA group ( Figure 1C ). Additionally, genes with differential mutation frequencies between IAC and AIS/MIA groups were selected for downstream analysis. Significant associations between α-diversity indices and mutation frequencies of specific genes suggested a potential link between the microbial community and genetic alternations ( Figure 1D ).

We identified six predominant microbial phylum in the tumor microenvironment, including Actinobacteria, Bacteroidetes, Firmicutes, Proteobacteria, Verrucomicrobia, and Fusobacteria. Increased proportions of Bacteroidetes and Firmicutes were observed in IAC group rather than AIS/MIA group ( Figure 2A ). At the genus level, 130 and 46 unique genera were identified in IAC and AIS/MIA groups, respectively, with 528 genera shared between the two groups. At the species level, 517 and 192 unique species were identified in IAC and AIS/MIA groups, respectively, and 1327 species were shared by the two groups ( Figure 2B ). Subsequently, α-diversity analysis revealed that IAC group was characterized by higher Shannon diversity and Simpson diversity scores than AIS/MIA group, despite with no significant difference ( Figure 2C ). Additionally, PCoA plot showed significant differences in tumor microbiota between IAC and AIS/MIA groups (P = 0.004) ( Figure 2D ). Several microbial genera detected in this study have previously been reported to be involved in lung cancer development, as revealed TSEA analysis, such as Blastomonas, Gemmatimonas, Mesorhizobium, Microbacterium, Mycobacterium, Mycoplasma, Parvimonas, Porphyromonas, Sphingomonas, Staphylococcus, and Veillonella ( Figure 2E ). Additionally, KO genes with increased abundance in the AIS/MIA group were significantly enriched in signaling transduction pathways, such as ABC transporters and bacterial chemotaxis, and metabolic pathways ( Figure 2F ). Subsequently, we explored interactions between microbes with significantly different abundance and genes with significantly different mutation frequencies. As a result, the interaction network uncovered the extensively positive correlations between specific microbes and KO genes, including Bosea sp. PAMC 26642 with K01432 and K02035, and Bosea sp. RAC05 with K01438 ( Figure 2G ).

Five predominant microbial phylum in the BALF microenvironment were detected, including Actinobacteria, Bacteroidetes, Firmicutes, Proteobacteria, and Fusobacteria. Compared to AIS/MIA group, IAC group showed increased Actinobacteria but declined Proteobacteria ( Figure 3A ). At the genus level, 47 and 122 unique genera were identified in IAC and AIS/MIA groups, respectively, with 320 genera shared between the two groups. At the species level, 217 and 337 unique species were identified in IAC and AIS/MIA groups, respectively, and 800 species were shared by the two groups ( Figure 3B ). Although there was no significant difference, we obtained higher Shannon diversity and Simpson diversity scores in AIS/MIA group rather than IAC group ( Figure 3C ). Slight difference between IAC and AIS/MIA groups was revealed by PCoA plot, and statistically significance did not reached ( Figure 3D ).

Based on the known relationships between microbes and specific diseases, we performed TSEA analysis, which suggested the involvement of some genus in lung cancer development, such as Achromobacter, Acinetobacter, Porphyromonas, Propionibacterium, and Veillonella ( Figure 3E ). KO genes with increased abundance in the IAC group were significantly enriched in more pathways associated with metabolism and signaling transduction, such as carbon metabolism, two-component system, and biosynthesis of cofactors ( Figure 3F ). Subsequently, we explored the interactions between microbes with significantly different abundance and genes with significantly different mutation frequencies. The interaction network uncovered the extensively positive correlations between specific microbes and KO genes, such as Bosea sp. PAMC 26642 with K03418 and K00449, and Bosea sp. RAC05 with K02031 ( Figure 3G ).

For the microbes with significantly different abundances between IAC and AIS/MIA groups, we visualized their profiles in two types of samples: tumor and BALF ( Figure 4A ). A total of 23 microbes exhibited consistent upregulation or downregulation in IAC group across both sample types ( Supplementary Table S1 ). Specifically, five microbes increased in IAC group, including Mycobacteroides abscessus, Mycolicibacterium aurum, Mycolicibacterium rhodesiae, Finegoldia magna, and Acinetobacter wuhouensis, while the remaining 18 microbes exhibited decreased abundance in IAC group ( Figure 4B ). The correlations among these 23 microbes were subsequently investigated. As expected, negative correlations were observed between upregulated and downregulated microbes, indicating potential antagonistic relationships. Of note, Mycobacteroides abscessus, which was upregulated in IAC group, displayed negative associations with all 18 declined microbes, suggesting its powerful influence ( Figure 4C ). The phylogenetic tree generated from the taxonomy analysis indicated that these 23 microbes primarily belonged to Proteobacteria, Alphaproteobacteria, and Rhizobiales taxa ( Figure 4D ). Among the five upregulated microbes in IAC group, two of them belonged to Mycolicibacterium genus, including Mycolicibacterium rhodesiae and Mycolicibacterium aurum. To further characterize their roles, we performed correlation analysis among microbes, KO genes, and KO pathways, and revealed the axes among Mycolicibacterium rhodesiae/Mycolicibacterium aurum KO0624 gene-KO04146 pathway, suggesting the potentially regulatory roles ( Figure 4E ).

To investigate the dysregulated characteristics within tumor microenvironment, we performed a correlation analysis between mutated genes and 23 significant microbes. Significant correlation values were widely observed between key microbes and mutated genes ( Figure 5A ). Of note, the strong correlations of TYW1, PTPRZ1, and GCNA genes with more than 20 different microorganisms suggested that these genes may play an crucial regulatory role in the microbial community of the lung. Among them, previous studies have found that PTPRZ1 is closely related to the occurrence and development of lung cancer. Through the joint analysis of multiple databases, the mutation situation of PTPRZ1 in lung adenocarcinoma is shown in Figure 5B . Mutations in the PTPRZ1 protein are closely related to TMB, Aneuploidy score, and Buffa hypoxia score ( Figures 5C-E ).