This prospective study was conducted from July 2018 to December 2019 in the adult ICUs of the First Affiliated Hospital, College of Medicine of Zhejiang University, China. The patient inclusion criteria were as follows: (1) patient was mechanically ventilated and hospitalized in the ICU; (2) Acute Physiology and Chronic Health Evaluation (APACHE) II score was greater than 12; and (3) collection of endotracheal deep aspirate (ETA) specimens was possible. Exclusion criteria were as follows: (1) age < 18 years; (2) comorbidity of chronic lung disease or lung cancer; (3) hospital stay < 24 h; and (4) co-infection with other bacteria. VAP was defined by the criteria of the Centers for Disease Control and Prevention (CDC) of the United States based on clinical, laboratory, radiological, and microbiological data (Horan et al., 2008). Patients meeting VAP criteria with positive culture for CRAB ETA were assigned to the CRAB VAP group (CRAB-I). Respiratory tract CRAB colonization (CRAB-C) was defined as CRAB-positive culture from ETA without VAP. Control patients (CRAB-N) were patients with neither CRAB VAP nor CRAB colonization. Approval was obtained from the ethical board of the hospital (reference number: 2016-458-1).

The aspirate samples were transported to the microbiology laboratory within 2 h. ETAs were inoculated with a calibrated loop (0.001 ml). A gram-stained smear was prepared for all specimens and examined microscopically. If the number of squamous epithelial cells is less than 10/ each low power field, it is judged as qualified. MacConkey agar (Oxoid, United Kingdom) and incubated aerobically at 37°C overnight. Acinetobacter baumannii was identified by using biochemical methods and matrix-assisted laser desorption-time of flight mass spectrometry (MALDI-TOF MS) (Bruker, Bremen, Germany). The confirmed A. baumannii from clinical and sputum samples were stored at −80°C for further study. The resistance of A. baumannii to antibiotics, including Piperacillin/tazobactam (TZP), Cefoperazone/sulbactam (CSL), Ceftazidime (CAZ), Ceftriaxone (CRO), Cefepime (FEP), Ciprofloxacin (CIP), Levofloxacin (LVX), Imipenem (IPM), Meropenem (MEM), Trimethoprim/sulfamethoxazole (TMX/SXT), Amikacin (AMK), and Gentamicin (GEN), was analyzed according to the recommendation of Clinical and Laboratory Standards (CLSI, 2018). Resistance to Tigecycline (TGC) and Polymyxin (POL) was examined by the broth microdilution method and evaluated according to the European Committee on Antimicrobial Susceptibility testing (2016). Escherichia coli ATCC25922 was used as a quality control strain. CRAB was defined as AB strains that were non-susceptible to imipenem or meropenem. For CRAB-C group, we selected samples and strains cultured for the first time for metagenomics and whole-genome sequencing. If the patients were in CRAB-I group, we selected the samples of the latest infection state and the cultured strains for sequencing.

Clinical data including demographic variables, length of ICU stay, length of hospital stay, comorbidities, previous invasive procedures (central line insertion, intubation, continuous renal replacement therapy, and surgery under general an aesthesia), as well as lengths and types of antibiotic treatments and severity of illness [Acute Physiology and Chronic Health Evaluation (APACHE II)], were recorded (LeGall et al., 1986). Data on the levels of serum markers for liver and renal function [e.g., bilirubin, alanine transaminase (ALT), aspartate transaminase (AST), gamma glutamyl transpeptidase (γ-GT), urea, and creatinine] as well as those of blood biomarkers of infection other than core temperature [e.g., white blood cell (WBC) and CRP] were also collected.

Total genomic DNA was extracted from samples using the cetyl trimethylammonium bromide (CTAB) method. We performed PCR amplification and pooled purified amplicons and then carried out the paired-end sequencing on an Illumina NovaSeq PE250 platform. After demultiplexing and trimming of the barcode and primer sequence using FLASH (V1.2.7) (Magoč and Salzberg, 2011), the paired-end raw read data of each sample was acquired. Subsequently, quality control was carried out using Qiime (Caporaso et al., 2010; Bokulich et al., 2013) and effective tags were obtained for analysis. OTUs, clustered with a 97% similarity cut-off using Usearch (Version 7.0) (Edgar et al., 2011), were taxonomically annotated using Mothur and SILVA database (Quast et al., 2013). MUSCLE (Version 3.8.31) software was used for rapid multiple sequence alignment and the phylogenetic analysis for all OTU representative sequences (Edgar, 2004). We then constructed the abundance matrices at the levels of phylum, class, order, family, and genus were constructed for each sample. Alpha diversity (Shannon and Simpson index within a sample) and Beta diversity (Bray–Curtis dissimilarity matrix across samples) were calculated using the phyloseq R package (version 1.32.0). PCoA and NMDS analysis based on Bray–Curtis dissimilarity was performed using the Vegan R package (version 2.5.6) and visualized using the ggplot2 R package (version 3.3.2). Analysis of similarity (ANOSIM) was used to test whether there is significant. We also calculated the co-occurrence relationships between bacteria using Python-based SparCC tool with SparCC correlation method and visualized using Cytoscape (Version 3.8.0) (Csardi and Nepusz, 2006). Sample names and the corresponding NCBI accession numbers are listed in Supplementary Table 1.

Following fragmentation of microbial DNA, metagenomic sequencing was performed on an Illumina NovaSeq 6000. After paired-end Illumina sequencing, we employed a previously reported bioinformatics pipeline to detect and profile the airway microbiome (Dickson et al., 2014; Kelly et al., 2016; Faner et al., 2017; Emonet et al., 2019). The low-quality sequences were filtered out or trimmed using PRINSEQ-lite (Version 0.19.3). And then de novo assemblies were generated using SPAdes genome assembler (Version 3.11.1) (Nurk et al., 2017) and coding sequences were predicted using MetaGeneMark (Version 3.38). Taxonomy assignments of both the clean reads and coding sequences were performed by Kaiju classifier (Version 1.7.2) (Menzel et al., 2016) with the National Center for Biotechnology Information Refseq database. We also utilized tool StrainPhlAn (Truong et al., 2017) to detected the multiple CRAB strains in ETA metagenomics samples. The functions of coding sequences were obtained using DIAMOND software (Version 0.9.30) with Kyoto Encyclopedia of Gene and Genomes (KEGG) database (Buchfink et al., 2015; Kanehisa et al., 2017). Furthermore, virulence genes in metagenomics sequences were identified by comparing the coding sequences against the Virulence Factor Database using DIAMOND software (Chen et al., 2016). The correlation coefficient between bacteria and virulence genes was generated using Python-based SparCC tool with SparCC correlation method. We performed multivariate linear regressions with feature selection, using the lasso penalized maximum likelihood technique in the “glmnet” R package (Version 4.0.2). Sample names and the corresponding NCBI accession numbers are listed in Supplementary Table 1.

The genomic DNA of 46 isolates was extracted using a Qiagen DNA purification kit (Qiagen, Hilden, Germany) and sequenced on an Illumina HiSeq 4000-PE150. For each isolate, de novo assembly of reads was performed using SPAdes (version 3.11) genome assembler. CheckM was used to check genome completeness and possible contamination into the genome (Parks et al., 2015). We annotated the assemblies using Prokka (version 1.14.6) (Tatusova et al., 2016). GFF format files produced by Prokka were subjected to Roary (version 3.11.2) (Page et al., 2015) to construct core genome. Meanwhile, we analyzed whether there were multiple CRAB ST clones in one ETA sample using the raw reads through MetaMLST (Zolfo et al., 2017). A maximum-likelihood phylogenetic tree was constructed using RaxML software (version 8.2.12) with 1,000 bootstraps replicates (Price et al., 2010). Average nucleotide identity (ANI) was calculated by using a pyani (version 0.2.10). ANI values above 95% between genomes of these isolates denote the same source (Jain et al., 2018). We then conducted MLST analysis according to the Institute Pasteur scheme (MLST-IP) and Oxford Database (MLST-OD) (Gaiarsa et al., 2019). Clonal complexes (CCs) were assigned by eBURST and were defined as single locus (Hall, 2013). CCs were named according to the number of the predicted founder ST. Furthermore, we detected Antimicrobial resistance genes (ARG) and AB virulence genes by comparing genome assemblies against the ResFinder antibiotic resistance gene database using the Abricate software (version 0.8) and against the Virulence Factor Database using the DIAMOND software (version 0.9.30), respectively (Seemann, 2017). We also determined the capsular polysaccharide (KL) and lipooligosaccharide outer core (OCL) synthesis of the A. baumannii using Kaptive software (version 0.5.1) (Wyres et al., 2020). Strain names and the corresponding NCBI accession numbers are listed in Supplementary Table 1.

The Linear Discriminant Analysis Effect Size (LEfSe) program was used to identify bacterial taxa and virulence genes that were differently abundant between sample types (Segata et al., 2011). Normally distributed continuous variables were represented as means ± standard deviation (SD) and compared using Student’s t-test, whereas non-normally distributed continuous variables were represented as median and interquartile range (IQR) and compared using Mann–Whitney U-tests. Categorical variables were compared by the χ² test or two-tailed Fisher’s exact test, as appropriate.

All sequencing data during the current study are available in the Sequence Read Archive (SRA). Metagenomics data and the 16S rRNA gene data are under BioProject PRJNA 681291, and the WGS data under BioProject PRJNA 679997.

A total of 101 patients were screened. According to the inclusion criteria, 64 patients were enrolled, and 52 patients completed the entire study protocol (Figure 1): 24 with CRAB-I, 22 with CRAB-C, and six with CRAB-N. The CRAB-I and CRAB-C patients did not differ in terms of age, sex, or severity indices (APACHE II scores), but C-reactive protein (CRP) and 30 day mortality were higher in CRAB-I (Supplementary Table 2). There was no significant difference in the prior antibiotic therapy between CRAB-C and CRAB-I. All CRAB isolates were highly resistant to all antibiotics except amikacin, polymyxin, and tetracycline (Supplementary Figure 1).

The Linear Discriminant Analysis Effect Size analysis using a logarithmic LDA score cut-off of 2.5 identified 46 and 55 different pathways in CRAB-C and CRAB-I, respectively, relative to CRAB-N, of which 40 and 45, respectively, increased. These pathways included fatty acid degradation, oxidative phosphorylation, nicotinate and nicotinamide metabolism, transport, porphyrin and chlorophyll metabolism, benzoate degradation, and biofilm formation in Vibrio cholerae (Figures 5A,B). Moreover, signaling proteins (KEGG pathway ko99995, Figure 5C) were more active in CRAB-C than in CRAB-I.

In the three groups of patients, we detected a total of 1,628 virulence genes, divided into 69 functional groups. The top 10 most abundant functional virulence groups in each group accounted for 26.5% of all factors, and the composition of the “toxicity” functional group differed between the three groups of patients (Figure 6A). These factors tended to be involved in iron uptake, siderophore biosynthesis, immune evasion, and biofilm formation. Ninety-five virulence genes (relative abundance > 0.1% for each) differed in abundance between CRAB-N and CRAB-C, and 105 differed between CRAB-N and CRAB-I (Supplementary Table 4). Eleven genes, including AB57_ 0984, AB57_ 0990, AB57_ 0992, and mymA, were more abundant in CRAB-I than in CRAB-C (Figure 6B and Supplementary Table 4). Virulence gene networks constructed from CRAB-I had a higher complexity index than those from CRAB-C and CRAB-N patients (56,206.99 vs. 44,722.75 vs. 11,052.12) (Figure 6C). More virulence genes associated with Acinetobacter were detected in CRAB-I than in CRAB-C and CRAB-N, with functions including immune evasion, iron uptake, and VFDB-unclassified (Figure 6D). Consistent with this association, the levels of VFs and the relative abundances of genera (Acinetobacter and Methylorubrum) or species exhibited a strong and significant positive correlation (R2 = 0.529, P = 1.1e-06; R2 = 0.755, P = 7.6e-12, Pearson’s correlation; Figures 7A,B), indicating that differences in the abundance of VFs were driven by differences in the species present in each group of patients.

We calculated quality of microbial genomes and found the completeness was 99.85% ± 0.27% and the contamination was 0.96 ± 0.43. The average nucleotide identity (ANI) of the 46 CRAB strains was >95%, indicating that they belonged to the same species (Supplementary Table 1 and Supplementary Figure 6A). MLST of 46 CRAB was dominated by ST208 (30.43%), followed by ST938 (15.22%). eBURST analysis revealed that seven ST types (87.5%) clustered in the same clonal complexes (CCs) (CC92) (Supplementary Figure 6B). The KL types were mainly KL9, KL2, KL93, and KL7, and all strains belonged to the OCL1 type (Figure 8). We observed no statistically significant difference in ST or KL type between the CRAB-I and CRAB-C groups (P = 0.478 and 0.444; respectively).

All CRAB isolates harbored bla_ADC–25, bla_OXA–23, and bla_OXA–66. All isolates harbored more than one oxacillinase gene, and 28 (60.8%) harbored the class A β-lactamase gene bla_TEM. The number of resistance genes did not differ significantly between the two groups (Figure 8). Annotation and analysis of virulence genes showed that strains from CRAB-I had more virulence genes than those from CRAB-C (Figure 8), and chi-square tests (Table 1) revealed that AB57_0990, Lpxl, and ABZJ_00085 were more abundant in CRAB-I, consistent with the metagenomics analysis (Figure 6B).