All participating patients with IPF were diagnosed between January 2011 and December 2013 at Asan Medical Center, Seoul, Republic of Korea and met the diagnostic criteria of the American Thoracic Society (ATS)/European Respiratory Society/Japanese Respiratory Society/Latin American Thoracic Association statement (1). Samples of lung tissues from patients with IPF (n = 20; 10 non-survivors [cause of death: AE = 1, disease progression = 2, unknown = 7]) were aseptically obtained at the time of surgical biopsy for diagnosis, and those from age- and gender-matched controls (lung cancer patients; n = 20) with no histological evidence of disease collected aseptically at the time of surgery were obtained from the Bio-Resource Center of Asan Medical Center. None of the patients with IPF or the controls had been treated with antibiotics, steroids, anti-fibrotic agents, or probiotics within 1 month before undergoing surgery. Lung tissues were procured under protocol #2016-1366. This study was conducted in accordance with the Declaration of Helsinki (2013) and was approved by the Institutional Review Board of Asan Medical Center (2018–1096). Written informed consent was obtained from all study participants.

Clinical and survival data of all patients were retrospectively collected from medical records, telephone interviews, and/or the National Health Insurance of Korea. Spirometry, total lung capacity (TLC) determined by plethysmography and diffusing capacity for carbon monoxide (DLco) measured according to published recommendations are expressed as ratios (%) of normal predicted values (23–25). The patients with IPF underwent 6-min walk tests (6MWT) according to the ATS guidelines (26). Baseline clinical data at the time of IPF diagnosis were collected within one month of sample acquisition.

Tissue samples were frozen in liquid nitrogen immediately after collection and stored at −80°C. Genomic DNA was extracted from lung tissues using Mo Bio PowerSoil® DNA Isolation Kits (Mo Bio Laboratories, Carlsbad, CA, USA) according to the instructions of the manufacturer. The variable V3 and V4 regions of the 16S rRNA genes were amplified using the following specific forward and reverse primers with overhang adapters: 5′-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGCCTACG GGNGGCWGCAG-3′ and 5′-GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGGACTACHVG GGTATCTAATCC-3′, respectively (27). The PCR proceeded using a 5 ng/μl DNA template, 2 × KAPA HiFi HotStart Ready Mix (KAPA Biosystems, Wilmington, MA, USA), and two amplicon PCR forward and reverse primers. The PCR protocol comprised initial incubation at 95°C for 3 min, followed by 25 cycles of 95°C for 30 s, 55°C for 30 s, and 72°C for 30 s, then 72°C for 5 min, and retention at 4°C. After PCR clean-up and index PCR, 300 bp paired-end sequences were pooled on the Illumina MiSeq platform (Illumina Inc., San Diego, CA, USA) as described by the manufacturer for all sample sequencing (28, 29). Distilled water provided in the PCR kit was used as a negative control, and no amplification was identified during the processes.

Fast Length Adjustment of Short (FLASH) reads, http://ccb.jhu.edu/software/FLASH/), were used for 16S rRNA gene by merging pairs of reads when the original DNA fragments were shortened than two times the length of the reads (30). Pre-processing and clustering were performed using the CD-HIT-operational taxonomic units (OTU; http://weizhongli-lab.org/cd-hit-otu/). Short reads (56,825) were filtered out and extra-long tails were trimmed. After filtering, the remaining reads were clustered at 100% identity. Chimeric reads (254,891) were filtered, and secondary clusters were recruited into primary clusters. After excluding reads with all other noise (5,292,165), the remaining reads (3,484,551) were clustered algorithm into OTUs at a cutoff of 97% (31, 32). Feature tables, such as abundance and representational sequence files, were created using UCLAST in Quantitative Insights Into Microbial Ecology (QIIME1; https://qiime.org) software (33). Taxonomy was assigned based on information about organisms with the closest similarity to the representative sequence of each OTU in the Basic Local Alignment Search Tool (BLAST), version 2.4.0, the NCBI 16S microbial reference database. Taxonomy was not assigned when the query coverage of the best match in the database was <85%, and the identity of the matched area was <85%.

Continuous data were analyzed using Mann–Whitney U tests, and categorical data were analyzed using Fisher exact tests. The decline rate of lung function and exercise capacity for one year was estimated by linear regression analysis. Correlations between the relative abundance of the microbiome and clinical parameters were assessed using Spearman's correlation coefficients (r). The risk of microbial relative abundance for a diagnosis of IPF was expressed as odds ratio (OR) with 95% CI using binary logistic regression. In addition, the risk of microbial relative abundance for IPF mortality was presented as hazard ratio (HR) with 95% CI using Cox proportional hazards regression analyses. Alpha diversity indices that estimate the number of unique OTU in each sample are represented using four indices; Observed estimated the actual number of different taxa evident in a sample, Chao 1 non-parametrically estimated the richness of the species (34), Shannon estimated richness and evenness of species present in a sample considering the distribution of strains belonging to each species (35), and Inverse Simpson measured the probability that two randomly selected objects in a sample belong to the same species (36). Principal coordinates analysis (PCoA), based on weighted UniFrac methods to obtain phylogenetic and quantitative indices for assessing abundance differences among groups (IPF vs. controls, survivors vs. non-survivors), was conducted for all samples using QIIME1 (37). The exploratory and differential microbial compositions were analyzed using QIIME1. All data were statistically analyzed using SPSS version 24.0 (IBM Corp., Armonk, NY, USA), and values with p < 0.05 (two-tailed) were considered statistically significant.

Among 20 patients with IPF, the mean age was 63.3 yr, and 65.0% were male (Table 1). Lung function (FVC, DLco, and TLC) was worse in the patients with IPF than in the controls, whereas demographics, lung function, and exercise capacity during the 6MWT did not significantly differ between IPF non-survivors and survivors.

Alpha diversity indexes, such as Observed, Chao 1, Shannon, and Inverse Simpson, did not differ between the IPF and control groups (Supplementary Figure A1). However, the PCoA plot revealed dissimilarity in the weighted UniFrac distance between the IPF and controls (Figures 1A–C), especially between five of the patients with IPF (non-survivors, n = 3; Figures 1A,C red circles) and controls, indicating more heterogeneity in the microbial distribution.

Among the 10 most frequent taxa, the genus Ralstonia was the most prevalent in the IPF and control groups, followed by Nocardia and Pelomonas (Supplementary Figure A2). On the contrary, Lactobacillus, Enterobacter, Tetragenococcus, and Neisseria were frequently identified in IPF, whereas Haemophilus, Caulobacter, Bradyrhizobium, and Thermomonas were prevalent in the controls. The relative abundance of Lactobacillus (0.91 [IPF] vs. 0.06% [control], p = 0.009), Paracoccus, (0.13 vs. 0.02%, p = 0.013), and Akkermansia (0.08 vs. 0.02%, p = 0.001) was higher, whereas that of Caulobacter (0.73 vs. 0.95%, p = 0.040), Azonexus (0.25 vs. 0.37%, p = 0.021), and Undibacterium (0.01 vs. 0.04%, p = 0.011) was lower in patients with IPF than in the controls (Figure 2). Logistic analysis adjusted by age and gender independently associated a diagnosis of IPF with lower relative abundance of the genera, Pelomonas (OR, 0.352; 95% CI, 0.139–0.891; p = 0.027), and Azonexus (OR, 0.013; 95% CI, 0.000–0.926; p = 0.046; Table 2).

Alpha diversity did not significantly differ between non-survivors and survivors of IPF (Supplementary Figure A3). However, the PCoA plot showed that the distribution of microbes differed between non-survivors and survivors (Figure 3). Among the 10 most frequent taxa, Ralstonia and Nocardia were the most common in both groups (Supplementary Figure A4). The genus Streptococcus was more abundant in non-survivors compared with survivors. In addition, the genera Neisseria, Haemophilus, Rothia, and Rubrobacter were frequently detected in non-survivors, while the genera Tetragenococcus, Enterobacter, Lactobacillus, and Caulobacter were prevalent in survivors. The relative abundance of genera Bifidobacterium (2.77 [non-survivors] vs. 0.68% [survivors], p = 0.003) and Olsenella (0.51 vs. 0.41%, p = 0.013) was significantly higher in non-survivors than in survivors (Figures 4A,B).

The median follow-up period for patients with IPF was 3.0 yr (interquartile range: 1.5–5.4 yr), and the median survival period was 3.1 yr. Unadjusted Cox analysis significantly associated the relative abundance of the Streptococcus, Sphingomonas, Veillonella, and Clostridium genera with IPF mortality. Neisseria and Granulicatella were also marginally associated with IPF mortality (Table 3). A multivariable model adjusted for age, and FVC selected a higher relative abundance of the Streptococcus (HR, 1.993; 95% CI, 1.019–3.901; p = 0.044), Sphingomonas (HR, 57.590; 95% CI, 1.714–1934.881; p = 0.024), and Clostridium (HR, 37.189; 95% CI, 1.228–1126.474; p = 0.038) genera as independent predictors of IPF mortality.

The relative abundance of Curvibacter and Thioprofundum was positively associated with FVC in patients with IPF, whereas Anoxybacillus, Enterococcus, Akkermansia, and Clostridium negatively correlated with FVC (Figure 5A and Supplementary Table A1). The relative abundance of Thermomonas and Peptoniphilus positively correlated with DLco, whereas that of Granulicatella and Rhodoferax was positively correlated with TLC.

The relative abundance of the Aquabacterium, Nakamurella, and Peptoniphilus genera was positively correlated, whereas that of the Fusobacterium, Anaerococcus, and Phycicoccus genera was negatively correlated with distance during the 6MWT (Figure 5B and Supplementary Table A2). The relative abundance of genus Rhodoferax and Lactococcus was positively correlated with resting and lowest oxygen saturation (SpO₂) during the 6MWT.

We estimated the decline rate in lung function and exercise capacity for one yr and compared them between survivors and non-survivors (Table 4). Non-survivors had a faster decline rate in DLco, and distance and the lowest SpO₂ during 6MWT compared with survivors. The relative abundance of Granulicatella and Paracoccus genera was positively correlated, while that of the Novosphingobium genus was negatively correlated with the decline rate in FVC (Figure 6A and Supplementary Table A3). The relative abundance of Bifidobacterium was positively associated, whereas Streptococcus was negatively associated with the decline rate in DLco. The relative abundance of Lactobacillus, Staphylococcus, Granulicatella, and Selenomonas genera was positively correlated with the decline rate of TLC.

The relative abundance of Staphylococcus and Variovorax was positively associated, while that of Legionella, Anoxybacillus, Acidocella, and Hyphomicrobium genera was negatively associated with the decline rate in distance during 6MWT (Figure 6B and Supplementary Table A4). The relative abundance of Beijerinckia, Mycobacterium, and Microbulbifer genera was positively correlated, whereas that of Enterobacter genus was negatively correlated with the decline rate in resting SpO₂. The relative abundance of genus Staphylococcus was positively correlated with the lowest SpO₂ during 6MWT.