Antibiotic Susceptibility, Virulome, and Clinical Outcomes in European Infants with Bloodstream Infections Caused by Enterobacterales

Mortality in neonates with Gram-negative bloodstream infections has remained unacceptably high. Very few data are available on the impact of resistance profiles, virulence factors, appropriateness of empirical treatment and clinical characteristics on patients’ mortality. A survival analysis to investigate 28-day mortality probability and predictors was performed including (I) infants <90 days (II) with an available Enterobacterales blood isolate with (III) clinical, treatment and 28-day outcome data. Eighty-seven patients were included. Overall, 299 virulence genes were identified among all the pathogens. Escherichia coli had significantly more virulence genes identified compared with other species. A strong positive correlation between the number of resistance and virulence genes carried by each isolate was found. The cumulative probability of death obtained by the Kaplan-Meier survival analysis was 19.5%. In the descriptive analysis, early age at onset, gestational age at onset, culture positive for E. coli and number of classes of virulence genes carried by each isolate were significantly associated with mortality. By Cox multivariate regression, none of the investigated variables was significant. This pilot study has demonstrated the feasibility of investigating the association between neonatal sepsis mortality and the causative Enterobacterales isolates virulome. This relationship needs further exploration in larger studies, ideally including host immunopathological response, in order to develop a tailor-made therapeutic strategy.

Very few data are available on the impact of different treatment regimens on clinical outcome in neonates with GN-BSIs. Previous studies conducted in both adults and children showed conflicting results on the impact of resistance profiles, appropriateness of empirical treatment and clinical characteristics on patients' mortality [7][8][9][10][11][12][13][14][15]. In the last decade, the implementation of modern bioinformatics to assist next-generation sequencing data analysis greatly improved the knowledge on the genetic characterization of pathogenic strains that may serve as target for new therapies [16]. There has been a growing body of evidence about the role of virulence factors (VFs) in the pathogenesis of invasive infections. Enterobacterales employ many strategies to enhance invasiveness, overcome host defenses and cause infections. Different strains can use alternative VFs with similar functions during the infection process, with this plasticity leading to unique combinations of such factors [17]. Some VFs are disease-specific whereas others seem to play different roles in different types of infection [18]. This plasticity enables pathogenic strains to colonize and infect different tissues and hosts. Some major classes of VFs, such as capsule, siderophores and fimbriae, have been characterized well [18]. However, several other factors were recently identified and have yet to be defined to fully understand their mechanisms of action and clinical significance (summarized in Table 1). To achieve this goal, a genomic approach can be used to identify genes encoding specific virulence determinants. In adults, several genes present in the great majority of bacteremic strains and involved in virulence have been identified [19][20][21]. These are presumably essential for the infection process. However, the specific VFs that are relevant in causing neonatal GN-BSIs are not well defined yet, partly due the variability among the few studies available so far [2]. These have mostly been conducted on neonates and children with Escherichia coli bacteremia, with virtually no data available on other Enterobacterales [22][23][24][25]. Also, in the recent years, several data have been published investigating a potential role of the bacterial virulome, defined as the set of genes contributing to the bacterial virulence, in determining the outcome of patients with both Gram-positive and GN-BSIs. Again, these studies demonstrate a mixed picture reporting a significant correlation between virulence factors and mortality in both children and adults [2,19,[26][27][28][29].  Clarifying the role of the main determinants leading to adverse outcomes could help to define targeted interventions to decrease mortality. With this study, we aimed to investigate potential associations between patient characteristics, pathogen characteristics and antibiotic treatment regimen on the clinical outcome of neonates/infants affected by culture-proven GN-BSIs.

Microbiological Data
The species IDs identified on matrix-assisted laser desorption/ionization time-of-flight (MALDI-TOF) mass spectrometry were all confirmed by sequencing. The isolate-specific accession numbers are indicated in Supplementary Table S1. The most frequently isolated pathogen was E. coli followed by Enterobacter spp. and Klebsiella spp. ( Table 2). The percentage of multidrug-resistance (MDR) isolates was 30% (26/87). Based on interpretation of the in vitro susceptibility profile, 16/87 (18%) were suspected of producing an extendedspectrum beta lactamase (ESBL) enzyme, and only one isolate (Klebsiella pneumoniae) was resistant to carbapenems. Susceptibilities of single species to the investigated antibiotics are presented in Table 3. A total of 50 different sequence types (STs) were found, with ST131 and ST90 as the most frequent in E. coli and Enterobacter cloacae, respectively. The median number of classes of resistance genes carried per isolate was 5 (IQR 4-5) whereas the median number of classes of virulence genes was 7 (IQR 7-9). Twenty-five isolates harbored blaTEM-type genes, two non-Klebsiella spp. strains the blaSHV-type determinant, and two E. coli strains the blaCTX-M-type genes. One K. pneumoniae ST17 carried blaVIM-12 gene, and two Enterobacter asburiae (ST484) the mcr-9 determinant.
Overall, the genome sequencing identified 299 different virulence genes among all the pathogens. There was a strong positive correlation between the number of resistance and virulence genes carried by each isolate (Rho = 0.79; p = 0.001) (Figure 1).
E. coli strains showed the highest mean number of virulence genes (105 vs. <65 in the other species overall), mainly those involved in fimbriae production (p < 0.0001) ( Table 4). The genes that were more frequently carried by the isolates are summarized in Table 5. Among the most represented strains (E. coli, E. cloacae, K. pneumoniae, K. oxytoca) the following virulence genes were carried by all the isolates: gal, gnd, rcs (capsule); ompA (cell invasion); ent, fep (iron metabolism); acr (pumps). Some genes were shown to be strainspecific. Among them, the adhesion mrk gene cluster was sequenced in all Klebsiella spp. isolates as well as the secretion system's exe and impA.tss genes. Conversely, genes encoding for motility and chemotaxis proteins (che, flg and fli) were only carried by E. coli and Enterobacter spp. The clb, esp and hly genes encoding for toxin proteins were sequenced only in E. coli strains. One hypervirulent K. pneumoniae with a hypermucoviscous phenotype harboring the rmpA and wca genes was found [33].
Overall, the genome sequencing identified 299 different virulence genes amon the pathogens. There was a strong positive correlation between the number of resis and virulence genes carried by each isolate (Rho = 0.79; p = 0.001) (Figure 1). E. coli strains showed the highest mean number of virulence genes (105 vs. <65 i other species overall), mainly those involved in fimbriae production (p < 0.0001) (Tab The genes that were more frequently carried by the isolates are summarized in Ta Among the most represented strains (E. coli, E. cloacae, K. pneumoniae, K. oxytoca) th lowing virulence genes were carried by all the isolates: gal, gnd, rcs (capsule); ompA invasion); ent, fep (iron metabolism); acr (pumps). Some genes were shown to be s specific. Among them, the adhesion mrk gene cluster was sequenced in all Klebsiella isolates as well as the secretion system's exe and impA.tss genes. Conversely, genes en ing for motility and chemotaxis proteins (che, flg and fli) were only carried by E. col Enterobacter spp. The clb, esp and hly genes encoding for toxin proteins were seque only in E. coli strains. One hypervirulent K. pneumoniae with a hypermucoviscous ph type harboring the rmpA and wca genes was found [33].

Determinants of 28-Day Case-Fatality
The cumulative probability of death obtained by the Kaplan-Meier survival analysis was 19.5% with the greater percentage of deaths happening in the first week. In the descriptive analysis, early age at onset (p = 0.002), culture positive for E. coli (p = 0.029), number of classes of virulence genes carried per isolate (p = 0.022) and GA (weeks) at the onset (p = 0.003) were significantly associated with mortality (Table 6). By Cox multivariate regression, none of the investigated variables was significant (Table 7).

Discussion
This study included 87 European neonates and infants younger than 90 days with GN-BSIs due to Enterobacterales. Overall, 299 virulence genes were identified in these pathogens. Among the different organisms, E. coli had significantly more virulence genes identified compared with other species. Gal, gnd, rcs, ompA, ent, fep, and acr virulence genes were identified from all the pathogens, likely being essential for the infection process. Conversely, some other genes were shown to be strain-specific. A strong positive correlation between the number of resistance and virulence genes carried by each isolate was found. By survival analysis, the 28-day probability of death was 19.5%. In the descriptive analysis, early age at onset, GA at the onset, culture positive for E. coli and number of classes of virulence genes carried by each isolate were significantly associated with mortality whereas discordant therapy was not related to mortality. By Cox multivariate regression, none of the investigated variables was significant.
Many studies have been conducted in both adults and neonates trying to define the main determinants of mortality in patients with GN-BSIs. AMR has been broadly investigated in the adult population, with the majority of studies reporting a significant association between multiple resistance to antibiotics and patients mortality [7,[34][35][36]. Some large cohorts have not found a clear correlation between AMR and adverse outcome [8,12,13]. Data from this small study did not confirm a significant impact of resistance profile on neonatal mortality [5].
In recent years, an increasing number of studies are being conducted trying to investigate the impact of virulence genes on the outcome of patients with GN-BSIs. Among them, E. coli was the most frequently investigated pathogen followed by K. pneumoniae. Independent risk factors associated with 30-day mortality among adult patients with ESBL-producing E. coli bacteremia included siderophores iroN and iss positivity [21], the siderophore fyuA gene, and the presence of the afimbrial adhesin afa gene [19,37]. In a large prospective study investigating the main determinants for adverse outcome in patients with K. pneumoniae BSIs, the cytotoxicity pks gene cluster carriage by causative strains was an independent risk factor for 30-day mortality when accompanied by MDR [38]. Lastly, the siderophore-related iutA gene was found to be an independent predictor of the 30-day mortality in K. pneumoniae bacteremia [39]. However, almost all of these studies were conducted with pre-selected virulence genes searched by polymerase chain reaction (PCR) rather than sequencing the entire bacterial virulome. This led to a wide heterogeneity, with each group analyzing different genes and hampering the comparison of the results. Very little data have been published on the relationship between bacterial virulence factors and BSI mortality in children, and almost all in patients with E. coli bacteremia. In a prospective cohort of 43 septic neonates, the adhesin hek/hra gene was found to be significantly more frequent in isolates from newborns who died than in isolates from survivors [2]. On the other hand, in a cohort of children 0-17 years old (median age 2.4 months), none of the 20 virulence factors tested by PCR was found to correlate with sepsis severity [26]. The Burden of Antibiotic Resistance in Neonates from Developing Societies (BARNARDS) study was conducted to assess the burden of AMR in neonates in seven low-middle income countries [6]. In this study, Gram-negative (GN) pathogens from neonatal sepsis were isolated and characterized through whole genome sequencing (WGS) for resistance and virulence genes. The number of virulence genes carried by each isolate through a virulence score was used. The results obtained suggested that yersiniabactin and/or aerobactin/salmochelin virulence genes may be involved in a more rapid onset and mortality. However, the inability to follow up all neonates and additional local factors likely to contribute to patient's death hampered the authors' capacity to attribute mortality singularly to the presence/absence of genomic traits.
Our results showed a strong positive correlation between the number of resistance and virulence genes carried by each isolate. Many studies have been conducted on either AMR or virulence. However, the biological effect and connection between these two factors are of particular importance [40,41]. Indeed, a negative or positive relationship can be found among them. Enhanced virulence or AMR frequently has been reported to have a fitness cost on bacteria but their relationship changes according to different bacterial species, the resistance and virulence genes involved and the host's immune system [42,43]. Some antibiotics, such as ceftazidime, cefotaxime and quinolones, have been reported to enhance the increase of the deletion and transposition of DNA regions that are specific for VFs [42]. By contrast, a positive correlation has been shown between AMR and virulence with the use of other antibiotics. In particular, uropathogenic strains of E. coli carrying the blaCTX-M-15 resistance gene also harbored more colV, colE2-E9, colIa-Ib, hlyA, and csgA genes as well as the blaOXA-2 beta lactamase was correlated with increased expression of colM, colB, colE, and crl genes [44]. Prophages are another mechanism that has been shown to be involved in both virulence and resistance diffusion through the spread of toxins to other resistant strains [45]. Porins and biofilm also play an important role in the relationship between virulence and resistance, with the first acting as a channel controlling the entrance of both antibiotics and VFs into the pathogens (i.e., the OmpC gene) and the second favoring antimicrobial treatment tolerance and infection persistence at the same time also increasing the transfer of resistance and virulence genes among the cells [46]. Lastly, an enhancement in both resistance and virulence characteristics of the pathogens can occur through mobile genetic elements like plasmids. These self-replicating extra-chromosomal elements are capable of transferring among different bacteria while carrying virulence and resistance genes [47]. This mechanism is independent of any antibiotic pressure.
There are several limitations in the study design, due to both confounding factors and heterogeneity of our sample. Firstly, the study is not powered to demonstrate a significant correlation between the presence of single resistance/virulence genes and neonates' clinical outcome; this required us to analyze resistome and virulome by grouping genes according to the class. Moreover, we are well aware of the rapid dynamics of bacterial genetics, and the selection of a single time point isolate can cause bias and affect the results. At the same time, considering the possible implication of heteroresistance, the selection of certain colonies in the first place could have potentially missed relevant isolates. We found a significantly higher number of virulence genes in E. coli isolates compared with other species; this could be due to an over representation of E. coli genes in the virulence database (VFDB). Different pathogens can have different impacts on neonates with sepsis, and pooling data on multiple strains could alter the results. Lastly, neonates and infants represent a broadly heterogeneous population, being characterized by different gestational ages, birth weight, underlying conditions, and risk factors. Despite these limitations, this was the first study to correlate virulence factors of non-E. coli Enterobacterales with mortality in septic neonates. Although a number of studies have tried to correlate single virulence genes with Gram-negative sepsis outcome, virtually none of them have investigated the role of the entire virulome in causing mortality. Considering the huge number of virulence genes involved and the limited sample size, we presented a potential way to use the entire body of information gained from the WGS by grouping genes as the number of classes of virulence genes involved. Having access to a very complete dataset including patient-level treatment, outcome data, and isolates available for a comprehensive genetic characterization, we believe that our data provide new and relevant information on the molecular picture of GN pathogens causing neonatal BSIs.

Study Design and Data Source
A case series of neonates with clinical sepsis and microbiologically confirmed GN-BSIs was constructed from three separate studies: the NeonIN [30], the NeoMero clinical trial [31], and the CLAHRC [32]. The NeonIN is a multinational network which prospectively collect data on neonatal infections from neonatal units. The detailed data procedures have been previously described [48]. NeoMero was a European-based randomized controlled trial to compare the efficacy of meropenem with standard of care (ampicillin + gentamicin or cefotaxime + gentamicin) in the treatment of late onset neonatal sepsis. CLAHRC is a United Kingdom-based prospective cohort study to collect data for patients with GN-BSIs in three hospitals in South London which aims to characterize clinical management of patients with GN-BSIs in all ages and identify potential risk factors associated with 28-day treatment outcomes.

Selection Criteria, Available Data and Definitions
Patients with a microbiologically-confirmed diagnosis of GN-BSIs due to Enterobacterales were selected from the above studies among European Neonatal Intensive Care Units between 2010-2015. Inclusion criteria were (I) age between 0-90 days, (II) Enterobacterales blood isolate available for the sequencing, (III) patient-level clinical and treatment data, and (IV) 28-day clinical outcome data.
Data available for the analysis included demographic, risk factors, clinical characteristics, antibiotic treatment, susceptibility results, resistome and virulome of the isolated bacteria, and 28-day treatment outcome.
The date of the blood culture was considered as the day of the sepsis onset. The MDR of the isolates was defined according to Magiorakos A.-P. et al. as acquired nonsusceptibility to at least one agent in three or more antimicrobial categories [49]. The first 48-h antibiotic treatment was defined as concordant or discordant based on the inclusion of at least an active drug against the blood isolate. The evaluated outcome was defined as 28-day case-fatality (patient alive/dead).

Microbiological Methods
Bacterial isolates from blood were cultured from frozen stocks (−80 • C) on blood agar plates. Species were identified by MALDI-TOF mass spectrometry (Bruker, Karlsruhe, Germany). Antibiotic susceptibility profiles were obtained according to the European Committee on Antimicrobial Susceptibility Testing (EUCAST) 2019 Clinical Breakpoints with disk diffusion tests for the following antibiotics: amikacin, ampicillin, amoxicillinclavulanate, ceftriaxone, ceftazidime, aztreonam, ciprofloxacin, gentamicin, meropenem, piperacillin-tazobactam, and trimethoprim-sulphametoxazole. To facilitate the analyses, isolates that were defined as having increased exposure were classified as non-susceptible. The isolates included in the study were subjected to WGS using the Illumina MiSeq platform (Illumina, San Diego, CA, USA), with paired-end runs of 2 × 300 bp, after Nextera XT library preparation. The obtained reads were assembled using SPAdes [50]. For each genome, we determined the ST using an in-house script (available upon request) and the Multilocus sequence typing (MLST) schemes and gene alleles sequences available on PubMLST (pubmlst.org). The isolates were further characterized at the genomic level with the identification of resistant and virulence genes using ABRicate (Seemann T, Abricate, Github https://github.com/tseemann/abricate, accessed on: 11 November 2019) and the following databases: The Comprehensive Antibiotic Resistance Database (CARD) [51] and Resfinder [52] for the resistance genes and VFDB [53] for the virulence genes.

Statistical Analysis
Qualitative variables were summarized by absolute frequencies and percentages, and quantitative variables by median and IQR. A descriptive analysis was conducted with the potential association between variables and outcome of interest evaluated by chi-squared or Fisher's exact test as more appropriate for qualitative variables, and Mann-Whitney or t-test as appropriate for quantitative variables. A Spearman's rank correlation was calculated to evaluate the correlation between the number of resistance genes and virulence genes carried by each isolate. We performed a survival analysis to investigate the 28-day mortality predictors using the Cox regression model (primary endpoint), after evaluated the proportional hazard (PH) assumption. In case of non-PH assumption, the weighted Cox regression model was performed [54]. A bivariate analysis (univariate) was carried out and the variables for which the p-value was <0.10 in univariate analysis were included in the multivariate model. All variables entered as covariates were evaluated at the baseline. As secondary endpoint, a survival analysis was conducted using the Kaplan-Meier curves to assess the probability of death at 28 days. To facilitate the analysis, the resistome and virulome data obtained in sequencing were categorized as number of classes of virulence and resistance genes carried by each isolate. Because of the huge heterogeneity of treatment regimens among the included patients, antibiotics were coded according to the WHO ATC/DDD Index 2020 at the 4th level [55]. p values <0.05 were considered as statistically significant.
All statistical analyses were performed using R Statistical Software (version 4.0.2) [56].

Ethics
The source studies were approved by the Ethical Committees of the participating institutions, and all enrolled patients' legal guardians provided informed consent. Given the retrospective nature of the present study, ethical approval for this analysis was not necessary.

Conclusions
To conclude, this pilot study demonstrated the feasibility of investigating the association between neonatal sepsis mortality and the causative Enterobacterales isolates virulome. The limited sample size of our cohort did not allow us to determine the role of single virulence genes in neonatal GN-BSIs but grouping genes as the number of classes involved allowed us to investigate the impact of the entire virulome in neonatal sepsis outcome. This knowledge may be useful for predicting clinical outcomes, detecting virulent strains, and helping with vaccine development. Expanding research on anti-virulence molecules together with the development of new antibiotics could be crucial to improving the management of these fragile patients. Further research would be advisable to elucidate the correlation with the timing of sepsis onset to personalize the clinical approach. These findings need further exploration in larger global studies, ideally including host immunopathological response, in order to develop a tailor-made therapeutic strategy.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/antibiotics10060706/s1, Table S1: Isolate-specific accession numbers. Funding: This study and the APC were funded by Penta Foundation.

Institutional Review Board Statement:
The source studies were approved by the Ethical Committees of the participating institutions. Given the retrospective nature of the present study, ethical approval for this analysis was not necessary.
Informed Consent Statement: All enrolled patients' legal guardians provided informed consent.

Data Availability Statement:
The corresponding author confirms that he had full access to all the data in the study and had final responsibility for the decision to submit for publication. Wholegenome shotgun projects have been deposited in GenBank (project code PRJEB44870). The isolatespecific accession numbers are indicated in Supplementary Table S1.