A total of 78 participants who came to the Department of Thoracic Surgery of the Third Affiliated Hospital of Harbin Medical University were recruited, including 43 patients with early-stage non-small cell lung cancer and 35 healthy relatives of these patients ( Table 1 ). Sixty-three serum (35 NSCLC and 28 healthy) and seventy-eight stool (43 NSCLC and 35 healthy) samples were collected. Fecal and serum samples were collected according to protocols approved by the local ethics committee, and written informed consent was obtained from all participants.

Collate clinical parameter information (including age, gender, body mass index (BMI), tumor stage, novel adenocarcinoma grade, smoking, family history, etc.), and exclude any unhealthy conditions by electrocardiogram and chest X-ray results. The main exclusion criteria were as follows: (1) ≤ 18 years old or > 80 years old, (2) individuals who had received antibiotics or probiotics in the past 3 months, (3) underlying diseases such as diabetes and hypertension.

Fecal and serum samples were collected in the morning after an overnight fast (≥8 h). The stool samples were divided into 3 equal parts (200 mg each), placed in sterile cryovials, and immediately transported to the laboratory for storage at -80°C. Blood samples were collected in coagulation tubes. After the blood was collected, the blood was gently mixed up and down for about ten times, and then centrifuged at 1800g for 10 minutes. The supernatant (serum) was collected in a 1.5 ml centrifuge tube, centrifuged at 13,000 g for 2 min, and the supernatant was transferred to a cryovial and stored at -80°C for further analysis.

DNA from different samples was extracted using the CTAB according to manufacturer ‘s instructions. Analyze the integrity and fragment size of the extracted DNA using 1% agarose gel electrophoresis. And NanoDrop 2000 (Thermo Scientific, USA) was used to measure the extracted DNA quality.

PCR amplification was performed using the following primers: 341F (5’-CCTACGGGNGGCWGCAG-3’) and 805R (5’- GACTACHVGGGTATCTAATCC -3’). The 5’ ends of the primers were tagged with specific barcodes per sample and sequencing universal primers. And then the PCR products were purified by AMPure XT beads (Beckman Coulter Genomics, Danvers, MA, USA) and quantified by Qubit (Invitrogen, USA). The amplicon pools were prepared for sequencing and the size and quantity of the amplicon library were assessed on Agilent 2100 Bioanalyzer (Agilent, USA) and with the Library Quantification Kit for Illumina (Kapa Biosciences, Woburn, MA, USA), respectively. The libraries were sequenced on NovaSeq PE250 platform.

Paired-end reads were assigned to samples based on their unique barcodes and truncated by cutting off the barcode and primer sequence. Merge paired-end reads using FLASH (Reyon et al., 2012). Raw reads were quality filtered according to fqtrim (v0.94) under specific filter conditions to obtain high quality clean labels. Chimeric sequences were filtered using Vsearch software (v2.3.4) (Caporaso et al., 2010). After dereplication using DADA2, the representative sequence with single-base accuracy is obtained, that is, the ASV (Amplicon Sequence Variants) feature table and feature sequence. Alpha diversity and beta diversity were calculated by QIIME2 after random normalization to the same sequences, and the graphs were drawn by R package. Blast was used for sequence alignment, and the characteristic sequences of each representative sequence were annotated with the SILVA database. LEfSe (Segata et al., 2011) analysis and Wilcoxon rank sum tests were used to identify genera that were differentially abundant between groups of subjects. Other diagrams were implemented using the R package (v3.5.2) and GraphPad Prism software.

Metabolites in serum samples were extracted using 80% methanol buffer. 400 µL of pre-chilled 80% methanol was added to 100 µL of the sample, vortex for 1 min, incubated for 5 min at room temperature, then overnight at -20°C. After centrifugation at 4000 g for 20 min, the supernatant was transferred to a new 96-well plate. QC samples were prepared by pooling together 10 μL of each extract. Metabolites were stored at -80°C prior to liquid chromatography-mass spectrometry (LC-MS) analysis (Want et al., 2006; Barri and Dragsted, 2013).

UHPLC-MS/MS analyses were performed using a Vanquish UHPLC system (Thermo Fisher, Germany) coupled with an Orbitrap Q ExactiveTMHF-X mass spectrometer (Thermo Fisher, Germany). The sample was injected onto a Hypesil Gold column (100 × 2.1 mm, 1.9 μm) using a 12-minute linear gradient at a flow rate of 0.2 mL/min. The eluents for positive polarity mode were 0.1% formic acid–water (A) and methanol (B). The eluents used for negative polarity mode were 5 mM pH 9.0 ammonium acetate (A) and methanol (B). The Q ExactiveTM HF-X mass spectrometer was operated in positive/negative mode with a spray voltage of 3.5 kV, a capillary temperature of 320°C, a sheath gas flow of 35 psi, an auxiliary gas flow of 10 L/min, S-lens RF class 60, Auxiliary gas heater temperature 350°C.

Statistical analysis was performed using statistical software R (R version R-3.4.3), Python (Python 2.7.6 version) and CentOS (CentOS 6.6 version). When the data were not normally distributed, area normalization was used for positive state transformation method.

These metabolites were annotated using the following databases: the KEGG database (https://www.genome.jp/kegg/pathway.html), the HMDB database (https://hmdb.ca/metabolites) and the LIPIDMaps database (http://www.lipidmaps.org/). Partial least squares discriminant analysis (PLS-DA) was performed in metaX (Wen et al., 2017). Metabolites with VIP > 1 and P value < 0.05 and fold change (FC) ≥ 1.2 or FC ≤ 0.833 were considered differential metabolites. Volcano plots were used to filter metabolites of interest based on log2 (fold change) ≥ 0.263 or log2 (fold change) ≤ -0.263, and -log10 (P-Value) metabolites from ggplot2 in R language. For cluster heatmaps, data were normalized using z-scores of regions of differential metabolite intensity and plotted in R by the heatmap package. The functions of these metabolites and metabolic pathways were investigated using the KEGG database. Metabolic pathway enrichment of differential metabolites was carried out.

Patient characteristics were expressed as mean ± standard deviation (SD), differences between groups were compared using the χ² test or independent samples t-test. Wilcoxon rank-sum test (for two groups) and Kruskal-Wallis test (for more than two groups) were used to compare differences among microbial groups. Student’s t-test and fold change analysis were used to compare metabolites between groups. The relationship between microorganisms and metabolites was assessed using Spearman rank correlation analysis. Values of P < 0.05 were considered as statistically significant.

To determine whether gut microbial changes were associated with early-stage NSCLC, we examined different groups of fecal microbiome samples, including 43 NSCLC patients and 35 healthy individuals, by 16S rRNA gene sequencing. All patients with non-small cell lung cancer are in the early stage and have not developed distant metastasis, including stage 0 (adenocarcinoma in situ, AIS) (2.33%), stage I (79.07%), and stage II (18.6%). The detailed clinical characteristics of all participants are shown in Table 1 . There were no significant differences in age, gender, smoking status body mass index (BMI) and family history between the two groups (P > 0.05).

Using amplicon sequence variants (ASVs) to track the dynamics of bacterial abundance in feces from different groups, Venn plots visualized the number of ASVs shared and unique between the healthy control (HC) group and the early-stage NSCLC group ( Figure 1A ). We found the two groups shared 1821 ASVs, and the early-stage NSCLC group had more unique ASVs than HC group ( Supplementary Table 1 ). And then we analyzed the community structure of gut microbes ( Supplementary Table 2 ). At the phylum level, Firmicutes, Bacteroidota, Proteobacteria, and Actinobacteriota were the main components in both the HC group and early-stage NSCLC group, with the abundance of Firmicutes and Proteobacteria being higher in the NSCLC group ( Figure 1B ). At the family level, compared to the HC group, the abundance of Lachnospiraceae, Bacteroidaceae, and Enterobacteriaceae were higher in early-stage NSCLC group, while the abundance of Bifidobacteriaceae, Prevotellaceae and Veillonellaceae were lower ( Figure 1C ). At the genus level, apart from the similar abundance of Faecalibacterium in both HC group and NSCLC group, Bacteroides and Escherichia-Shigella were slightly more abundant in early-stage NSCLC group, while Bifidobacterium, Megamonas, Prevotella_9 and Dialister were relatively lower ( Figure 1D ).

Next, statistical analysis of microbial abundance was performed. Both early-stage NSCLC group and HC group showed comparable numbers of observed OTUs (operational taxonomic units). The Shannon and Simpson indexes both showed that community diversity was similar among the two groups. The Chao1 index showed no significant differences in community richness between early-stage NSCLC group and HC group. These data suggest that global community alpha diversity is similar between early-stage NSCLC group and HC group ( Figure 2A ). When comparing microbial community structure, beta diversity showed differences between the two groups ( Figure 2B ).

We next compared gut microbes with significant differences in expression abundance between HC and early-stage NSCLC groups at the phylum and genus levels. In total, 1 phylum and 12 genera were significantly decreased in the abundance of early-stage NSCLC patients ( Figure 3B ), while 14 genera were conversely enriched ( Figure 3A ). Desulfobacterota, the only phylum with significant differences between these two groups, was more abundant in the HC group. At genus level, Agathobacter, Blautia, Clostridium, an uncharacterized genus of family Muribaculacea, Cetobacterium, an uncharacterized genus of family Pasteurellaceae and eight other genera were significant abundant in early-stage NSCLC than in HC group, whereas Lanchnoclostridium, Prevotella, Lachnospia, Catenibacterium, Oscillospira, UGG-003, Lachnospiaceae_UGG-010, an uncharacterized genus of family vandinBe97, Acidaminococcus, Prevotellaceae_NK3B31_group, Oxalobacter were significantly enriched in the HC group.

Linear discriminant analysis (LDA) effect size (LEfSe) (Segata et al., 2011) was then used to generate cladograms to reveal differences in taxa abundance between early-stage NSCLC and HC ( Figure 4A ). There were 25 and 8 bacterial taxonomic clades that were significantly different in HC and early-stage NSCLC groups, respectively [log10 (LDA score) > 3] ( Figure 4B ). We found that Clostridia class was significantly higher in the early-stage NSCLC group. Agathobacter and Blautia were the prominent gene level biomarkers for early-stage NSCLC group. For healthy controls, the Desulfovibrionia and Negativicutes were the abundant class, and Prevotella_9, Prevotella, Lachnospira and Catenibacterium were the most prominent genus level biomarkers. Overall, these findings demonstrated that the early-stage NSCLC group had relatively lower microbial abundance than the HC group and was sufficient enough to distinguish healthy individuals from early-stage NSCLC patients.

Adenocarcinoma was the main type of pathology in NSCLC (38/43, 88.37%). In this study, after ruled out adenocarcinoma in situ (1 case) and minimally invasive adenocarcinoma (8 cases), we classified invasive pulmonary adenocarcinoma (IPA) into three groups (grade 1, n =7; grade 2, n =12; grade 3, n =10) according to the new grading system proposed by the IASLC (Moreira et al., 2020; Deng et al., 2021). Differences in microbial composition between healthy patients and IPA patients with different grades under the new grading system were analyzed at the phylum and species levels ( Supplementary Figure 2 ). Relative abundance analysis showed clear distinctions between high and low abundance taxa, and used color gradients to reflect similarities and differences in the composition of multiple samples at each taxonomic level. As shown in Supplementary Figure 2 , according to the change of the color gradient, the differences between the four groups of samples can be seen intuitively. The data showed that the dominant flora of IPA patients in each group was different from that of healthy people, suggesting a correlation between the florae features and the histopathological process of invasive lung adenocarcinoma.

Next, we analyzed biomarkers between IPA patients with different grades and healthy controls by multi-level LEfSe ( Figure 4C ). There were significant differences in 7, 7, 3, and 4 bacterial taxonomic clades in the healthy group and in the invasive pulmonary adenocarcinoma grades type 1, type 2, and type 3 group, respectively [log10 (LDA score) > 3] ( Figure 4D ). The key species were Erysipelatoclostridium in HC; Blautia and Marinobacter in IPA grade type 1; Dorea in IPA grade type 2; and Agathobacter in IPA grade type 3. These results showed that the fecal gut microbiota was specific for a novel graded type of invasive lung adenocarcinoma.

Previous studies have revealed that gut microbiota has a significant impact on blood metabolite profiles (Wikoff et al., 2009; Wilmanski et al., 2019). To further explore changes in gut microbe-host interactions, we performed LC-MS/MS-based non-targeted metabolomic analysis of serum from healthy individuals and patients with early-stage NSCLC. A total of 866 metabolites were identified and quantified, including 553 positive ions and 313 negative ions ( Supplementary Table 3 ).

Supervised multivariate statistical analysis using partial least squares discriminant analysis (PLS-DA) to maximize screening for differential metabolites across groups. The PLS-DA score plot showed a clear separation between the HC group and the early-stage NSCLC group ( Figure 5A ). Permutation tests indicated that the data were not overfit, the R2Y and Q2 values were 0.87 and -0.4, respectively, validating the OPLS-DA model ( Figure 5B ). The differentially expressed metabolic ions are screened by p value of the t-test and variable difference contribution (VIP), where VIP≥1.0, P< 0.05 as the filter condition. A total of 123 different metabolites were identified in serum between the HC and early-stage NSCLC groups ( Supplementary Table 4 ), most of which were upregulated. Figures 6 revealed the changes in these metabolites. In the HC group, the abundant metabolites mainly included carbohydrates and carbohydrate conjugates (Deoxyribose 5-phosphate), steroids and steroid derivatives (Testosterone sulfate, Glycocholic acid, D-erythrose 4-phosphate) and amino acids, peptides and analogs (Ergothioneine). In contrast, metabolites with higher levels in the early-stage NSCLC group mainly included sphingolipids (D-erythro-sphingosine 1-phosphate, palmitoyl sphingomyelin), fatty acyl (Avocadyne 1-acetate, 12(S)-HETE, 20-Carboxy-Leukotriene B4, Thromboxane B3, 6-Keto-prostaglandin f1alpha, Sebacic acid, Tetradecanedioic acid) and glycerophospholipids (LPC 20:2, LPC 18:0, LPC 18:4, LPE 20:2, LPC 20:1, LPC 16:1, LPC 20:0, LPA 18:2, LPC 17:1, LPC 17:2, LPC 19:0). In order to show the relationship between samples and the expression differences of metabolites in the two groups more intuitively, we performed hierarchical clustering analysis, and the results for the top 50 metabolites of p values with significant differential expression were shown in Supplementary Figure 3 .

The KEGG pathway enrichment analysis was then performed on the differentially abundant metabolites ( Supplementary Table 5 ). The results showed that the differential metabolites of early-stage NSCLC and HC were mainly involved 20 pathways ( Figure 7 ), including sphingolipid metabolism, sphingolipid signaling pathway, primary bile acid biosynthesis, the pentose phosphate pathway, carbon metabolism, arginine biosynthesis, phenylalanine, tyrosine and tryptophan biosynthesis, etc.

To further investigate microbiota-metabolite interactions associated with early-stage NSCLC, we assessed correlations between 27 genera and 32 metabolites ( Figure 8A ; Supplementary Table 6 ). The results showed that the abundance of several microbial genera in the early-stage NSCLC group were positively correlated with serum metabolite levels (Sperman’s correlation analysis, P<0.05, Figure 8A ). Then, based on the above microbiome data, a co-occurrence network was constructed to elucidate the major interactions between the early-stage NSCLC associated microbiome and metabolites ( Figure 8B ). The results showed correlations between Muribaculacead, Clostridium, Blautia, Agathobacter and the related metabolites. From the graph, Muribaculacea and Clostridium seemed to be the core genera given that they were positively correlated with metabolites enriched in early-stage NSCLC and negatively correlated with certain metabolites enriched in HC (eg, Deoxyribose 5-Phosphate and Testosterone sulfate).