Bacterial ribosomal RNA detection in cerebrospinal fluid using a viromics approach

Patients participated in the MeninGene study, a nationwide prospective cohort study of community-acquired bacterial meningitis in the Netherlands, methods of which have been described elsewhere [19, 20]. Briefly, patients with a positive CSF culture were identified by the Netherlands Reference Laboratory for Bacterial Meningitis (NRLBM), which receives the cultured pathogen from 85% of bacterial meningitis patients in the Netherlands. The NRLBM notified the researchers, who contacted the treating physician, who subsequently informed patients or their legal representative about the study. Patients could also be included by their treating physician without notification by the NRLBM. All patients or representatives gave written informed consent, and the study was approved by the Medical Ethics Committee of the Amsterdam UMC (METC2013_043). Clinical data were collected using an online case record form and patient outcome was recorded using the Glasgow Outcome Scale [21]. Leftover CSF was stored at treatment centers at − 80 °C and transferred to the Amsterdam UMC biobank facility. For this study, 89 CSF samples with sufficient residual material were selected. Researchers performing library preparation, sequencing and metagenomic analysis were blinded to patient clinical information and the diagnosed pathogen. Subsequently, data were unblinded for optimization of threshold parameters. For controls, previously generated [22] sequencing data from 74 patient CSF samples tested negative in culture for bacteria were included in analysis, approved by a separate decision of the Medical Ethics Committee of the Amsterdam UMC (METC2014_290). These samples were from patients undergoing lumbar puncture for suspected CNS infection, and either viral CNS infection was diagnosed or CNS infection was eventually ruled out.

VIDISCA library preparation was performed on CSF as previously described [14]. Briefly, 110 µl of CSF was centrifuged for 10 min at 5000 g, and supernatant was treated with TURBO DNase (Thermo Fisher Scientific) for 30 min at 37 °C. Nucleic acids were extracted using the Boom method [23], followed by reverse transcription primed with non-ribosomal hexamers [17] and second strand synthesis using Klenow fragment (3′ → 5′ Exo-, NEB). After clean-up by phenol/chloroform extraction and ethanol precipitation, dsDNA was digested with Mse1 (NEB) and ligated to sample specific adapters. Size selection with AMPure XP beads (Agencourt) was done to remove small DNA fragments. After a 28-cycle PCR, further size selection was done to retain fragments 200–600 bp long. Library concentrations were analysed using the Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific), and pooled at equimolar concentration. Pool concentration and fragment length distribution were analysed using the Qubit and Bioanalyzer (High Sensitivity Kit, Agilent Genomics) instruments respectively. Sequencing was carried out on the Ion S5 System with the Ion 510 Chip Kit (Thermo Fisher Scientific).

For bacterial rRNA identification, reads were mapped to the SILVA 138.1 SSU and LSU NR99 rRNA databases [24] using BWA MEM v0.7.17-r1188 [25]. Outputs were processed using the PathoID module of PathoScope v2.0.7 [26], and hits to phylum Chordata were removed. Remaining reads were realigned to the GenBank nt database using BLASTn [27], and hits to the five bacterial pathogens included in this study were counted. Blinded diagnostic predictive ability was explored by selecting the pathogen with the highest read count per sample as the predicted species, and those with equal counts as indeterminate. Unblinded detection performance per pathogen was then measured via sensitivity and specificity calculations. This was repeated varying several threshold parameters to understand their impact and optimize detection performance; parameters were minimum pathogen specific read count (≥ 1 read versus ≥ 10 reads), pathogen read identity to reference (≥ 97% versus 100%), and sample sequencing depth (all samples versus samples ≥ 10,000 reads). The level of human background per sample was estimated by mapping reads to a human rRNA database, subsetted from the aforementioned SILVA database. Hits were realigned to the same database using BLASTn, with reads retained and counted if they matched with 100% nucleotide identity for at least 100 bp.

Sources of bacterial reads (rRNA versus non-rRNA) were assessed by mapping them to genome assemblies of Streptococcus pneumoniae (GCF_002076835.1), Neisseria meningitidis (GCF_008330805.1), Listeria monocytogenes (GCF_000196035.1), Haemophilus influenzae (GCF_004802225.1), and Klebsiella pneumoniae (GCF_000240185.1). This was done first using original assemblies, and then versions with rRNA genes masked by RepeatMasker v4.1.1 [28]. Reads were curated by realignment to the respective original reference using BLASTn (requiring 100% identity for ≥ 100 bp). The rRNA count was calculated by subtracting the masked assembly read count from the original assembly read count. To determine if non-rRNA was predominately from mRNA or intact gDNA, the proportion mapping to reference coding sequences from each pathogen was calculated, with quality filtration as above. The proportion of coding sequence reads was compared with the gene density of each genome (coding sequence length/total genome length), since this represents the expected proportion of coding sequence reads derived from randomly sequenced pure gDNA.

Detection of viruses was done using a simplified version of a previously published workflow [29]. Reads were aligned to a database of viral proteins, with hits realigned to the GenBank nt database using BLASTn. Those aligning to non-viral sequences were removed as false positives, while those aligning either to viruses or to no reference were retained for manual examination. Reads from viruses of interest were aligned to respective reference genomes (Enterovirus D68: AY426531.1, HIV-1: AF286365.1) in Geneious v4.8.5 (https://​www.​geneious.​com) to visualize genomic position and coverage. HIV-1 reads were subtyped using the Los Alamos National Laboratory HIV BLAST tool (www.​hiv.​lanl.​gov).

We explored whether clinical characteristics were associated with detection of bacteria using Fisher’s exact tests for dichotomous data or Mann–Whitney U tests for continuous data. Sequencing depth was correlated to clinical characteristics using Spearman’s rank correlation. A p-value of < 0.05 was considered statistically significant.

Between 2006 and 2022, 2705 patients with bacterial meningitis were included in the MeninGene study. We selected 89 patient samples from this cohort for VIDISCA analysis (Additional file 2: Table S1). Culture based pathogen diagnoses were S. pneumoniae (n = 24), N. meningitidis (n = 22), L. monocytogenes (n = 24), H. influenzae (n = 18), and K. pneumoniae (n = 1). Thirty-nine patients were female (44%) and the median age was 58 years (interquartile range (IQR) 34–67). Nine patients (10%) were being treated with antibiotics at the moment of presentation. The most frequent symptoms on presentation were fever in 72 patients from 81 case record reports (89%) and headache in 65 patients from 81 reports (80%). An altered mental status, defined as a Glasgow Coma Scale (GCS) score < 14, was seen in 47 patients (53%). Ten patients (11%) were in a comatose state (GCS < 8) upon presentation. Median number of CSF leukocytes was 2560/mm³ (IQR 768–5680) with a median CSF total protein concentration of 3.1 g/L (IQR 1.8–5.5). Seventy-three patients (82%) had a favorable outcome, defined as a score of five on the Glasgow Outcome Scale at discharge, and six died (7%). The negative controls consisted of 38 patients with viral CNS infection and 36 patients with initial suspicion of CNS infection, eventually ruled out.

We next explored the impact of threshold parameters on overall and per pathogen diagnostic performance. Alignment identity cut-off, sequencing depth, and pathogen specific read count impacted test diagnostic accuracy, with overall sensitivity to any pathogen ranging from 40% (30–51% CI) to 69% (56–79% CI), and specificity from 87% (82–90% CI) to 99% (97–100% CI; Fig. 1, Additional file 2: Table S2). For individual pathogens, there was high sensitivity to S. pneumoniae, N. meningitidis, and H. influenzae across parameters, and poor sensitivity to L. monocytogenes (Fig. 1, Table 1). The single K. pneumoniae positive sample was successfully detected across all parameter thresholds. Specificity was high across pathogens, with the exception of certain parameter thresholds for S. pneumoniae. Diagnostic performance was higher when only samples with ≥ 10,000 reads were considered, driven by increased sensitivity at low specificity cost, though the impact varied from minimal (e.g. N. meningitidis) to large (e.g. S. pneumoniae) (Fig. 1). For minimum pathogen read count, requiring ≥ 10 reads reduced false positive rates for S. pneumoniae from ~ 50% to ~ 10% when compared to ≥ 1 read (Fig. 1). The impact was lower for L. monocytogenes, H. influenzae, and N. meningitidis, which already had low false positivity rates at the ≥ 1 read threshold (≤ 11%). A ≥ 1 read threshold increased sensitivity for all pathogens compared with ≥ 10 reads, particularly for L. monocytogenes and H. influenzae. Alignment identity requirement had minimal overall impact, though in some cases the ≥ 97% cut-off increased sensitivity compared to the 100% cut-off at low to no specificity cost (Fig. 1), leading us to select this as the universal cut-off.

To further understand the differences in diagnostic performance between pathogens, we produced a scatterplot of individual rRNA read counts for each, including all samples (Fig. 2). This highlighted a universally low pathogen read count for L. monocytogenes with a maximum of 15 reads, which likely contributes to the 21% sensitivity we observed (5 of 24 positive samples detected). Notably, 37% of L. monocytogenes culture positive samples were low-depth (< 10,000 reads), compared to a cohort average of 25%. Both L. monocytogenes and H. influenzae had zero false positive detections, while both K. pneumoniae and N. meningitidis had five, all with < 10 reads (Fig. 2). S. pneumoniae was frequently detected in samples from patients with bacterial meningitis caused by different pathogens, although 26 of 31 (84%) of these had < 10 reads and all were < 250 reads. Analysis of control CSF samples revealed similar false positive patterns (Additional file 1: Figure S1), with no L. monocytogenes false positives, though this time four false H. influenzae detections were made, all with < 10 reads. Six and two N. meningitis and K. pneumoniae false positives were found respectively (all < 10 reads). Again S. pneumoniae carried the highest false positive risk, with 28 false positives, 26 of which had < 10 reads, and all of which had < 19.

From exploratory testing, we found detection of the CSF culture diagnosed bacterium by VIDISCA was associated with no antibiotic treatment prior to presentation [49 of 78 (63%) versus 2 of 9 (22%); p = 0.03]. Detection by VIDISCA was lower in patients on immunosuppressive therapy (2 of 12 (17%) versus 48 of 79 (61%); p = 0.004), and higher in patients presenting with an altered mental status [GCS score < 14; 32 of 47 (68%) versus 19 of 42 (45%)] or a comatose state [GCS score < 8; 9 of 10 (90%) versus 42 of 79 (53%); p = 0.04]. GCS score was lower in cases accurately detected by VIDISCA, with 12 (IQR 9–15) versus 15 (IQR 12–15; p = 0.05). The level of C-reactive protein (CRP) was higher in cases detected by VIDISCA, with 206 mg/L (IQR 101–353) versus 91 mg/L (IQR 41–148; p = 0.003). The number of CSF leukocytes did not differ between groups, but CSF protein was higher in VIDISCA detected cases, with 4.1 g/L (IQR 2.9–6.0) versus 2.1 g/L (0.9–3.7; p < 0.001). No correlations between raw read count and clinical variables were found. Raw read count also did not correlate with human rRNA read count as a percentage of the total (Spearman’s rho = 0.16, p = 0.14), but the latter was weakly correlated with CSF leukocyte count (Spearman’s rho = 0.29; p = 0.006).

We hypothesized bacterial reads would primarily derive from rRNA and not residual gDNA, due to our library preparation methods. Generally, rRNA did make up the major fraction of detected pathogen reads (Fig. 3), especially for L. monocytogenes, N. meningitidis, and S. pneumoniae samples (all with median values ≥ 87%). S. pneumoniae was also notable in that some samples had very high rRNA counts, which were not seen in samples containing other bacteria. Across species however, substantial read fractions from bacterial-non-rRNA were also found. In particular, H. influenzae had a low median value of 36% rRNA, and the single K. pneumoniae sample had only 10% rRNA, with a high overall bacterial-read count. To determine the likeliest source of non-rRNA reads (mRNA versus residual gDNA) we determined the proportion of non-rRNA aligning to coding sequences of respective pathogen genomes. In most cases this proportion was consistent with a predominately mRNA source, being higher than the expected value for randomly sequenced gDNA (Additional file 1: Figure S2). For H. influenzae however, the median proportion was precisely the expected value for gDNA, suggesting DNase treatment or centrifugation may have underperformed in CSF samples containing this species.

Viral metagenomic analysis identified six CSF samples positive for at least one human virus (Table 2). These belonged to four families (Flaviviridae, Anelloviridae, Picornaviridae, and Retroviridae). Reads from Enterovirus D68 (EV-D68, family Picornaviridae) were identified in one sample, from three regions of the viral genome (Additional file 1: Figure S3A). The patient, who was taking prednisone, presented reporting gastrointestinal symptoms for several days and confusion on the day of presentation, followed by an epileptic seizure, and was diagnosed with L. monocytogenes meningitis. Separately, two distinct HIV-1 reads were detected in the CSF sample of a patient with meningococcal meningitis, both overlapping regions of the Env gene (Additional file 1: Figure S3B). Alignment to subtyped reference genomes showed the reads belonged to HIV-1 subtype B. HIV-1 infection in this patient had been discovered just prior to their admission due to meningitis, and antiretroviral therapy had not yet been started. The patient’s CD4 count was 160 cells/mm³, with a serum viral load of 150,000 copies/mL.