Longitudinal monitoring of DNA viral loads in transplant patients using quantitative metagenomic next-generation sequencing

Patients and sample selection

Six immunocompromised patients (one allogeneic stem cell transplant patient, four kidney transplant patients, one haematological patient) were selected based on available follow-up EDTA plasma samples that previously tested positive for one or more transplantation-related DNA viruses. Samples had been previously (July 2008 – December 2019) sent to the Clinical Microbiological Laboratory (CML) of the Leiden University Medical Center (LUMC, The Netherlands) for viral load monitoring. The samples were respectively qPCR positive for ADV, BKV, CMV, EBV, B19V, and TTV, with a wide range of viral loads. Patient plasma samples were stored at -80°C until mNGS analysis.

Ethical approval

Approval was obtained from the ethical committee from the LUMC (P11.165 NL 37682.058.11, and Biobank Infectious Diseases protocol 2020-03 & 2020-04 B20.002).

Extraction of nucleic acids, internal controls

Patient plasma samples were spiked with an internal control (baculovirus, Arc Bio, LLC) prior to extraction. Nucleic acids were extracted from 200 μl plasma using the MagNApure 96 DNA and Viral NA Small volume extraction kit on the MagNAPure 96 system (Roche Diagnostics, Almere, The Netherlands) with 100 μL output eluate. The eluate was concentrated using vacuum centrifugation by a SpeedVac vacuum concentrator (Thermo Scientific) to a volume of 26 µl.

Library preparation and sequencing

Sequence libraries were prepared using the Galileo Viral Panel sequencing kit (Arc Bio, LLC., Cambridge, MA, USA) according to the manufacturer’s instructions. The protocol is based on enzymatic fragmentation at 37°C for 5 minutes, followed by end repair and A-tailing at 65°C for 30 minutes. Subsequently, fragments were ligated using unique dual-index adapters at 20°C for 15 min and purified using magnetic Kapa Pure Beads (Roche). Human DNA was depleted using human depletion reagents at 45°C for 2 hours followed by 45 °C for 15 minutes, after which libraries were amplified using library amplification primers for 45 °C for 30 seconds, by 14 cycles of 98°C for 10 seconds and 65°C for 75 seconds and 65°C for 5 minutes. The final library preparation products were purified using magnetic Kapa Pure Beads (Roche) and quantified using a Qubit fluorometer (Thermo Fisher) followed by equally pooling using the Arc Bio calculation pooling tool. After a final quantity and quality check using a Bioanalyzer (Agilent), samples were sequenced using the NovaSeq6000 sequencing system (Illumina, San Diego, CA, USA) at GenomeScan B.V. (Leiden, the Netherlands) aiming at 10 million reads per library.

Calibration samples

Initial calibration runs were performed testing the multi-analyte mixture (MAM) of whole-virus particles at viral loads of 0, 1,000, 5,000, 10,000, and 100,000 copies/ml or IU/ml plasma, in quintuple (Arc Bio, LLC) for the following 10 viruses: ADV, BKV, CMV, EBV, HHV-6A, HHV6B, HSV-1, HSV-2, JCV, and VZV. For TTV and B19V, no Arc Bio calibrator panels were available, and therefore the Galileo Signal values were plotted against the calibrator plot of other viruses that showed optimal agreement with the viral load (respectively JCV and VZV), representing a semi-quantitative result.

Bioinformatic analysis

After demultiplexing of the sequence reads using bcl2fastq (version), FASTQ files were uploaded to the Galileo Analytics web application [13][15] which automatically processes data for quality assessment and pathogen detection using a custom database of DNA viruses involved in transplant-associated infections: ADV, CMV, EBV, HHV-6A, HHV-6B, HSV-1, HSV-2, JCV, VZV, B19V and TTV.

Galileo Analytics web application aligns sequence reads to the genomes of the DNA viruses in their calibration kit and scores these read alignments based on complexity, uniqueness and alignment scores and reports this in a signal value. The signal value is normalized for read counts across libraries, correcting for differences in genome lengths and technical bias, based on the spiked-in normalization controls. The Signals reported are related to the genomic depth and the observed amount of viral DNA being present in a sample, belonging to non-confounding genomic regions [13]. The sample Signals were visualized in linear calibration curves (Supplementary Figure 1).

• Download figure
• Open in new tab

Supplementary Figure 1.

Calibration graphs of the six viruses in six patients in this study with associated slope, intercepts and R² values. Concentrations are expressed in log₁₀ copies or IU/ml. Calibrator samples are shown in black dots, clinical samples in orange.

Analysis of performance, and additional findings

Performance of the metagenomic Galileo Viral Panel assay was assessed in comparison with routine qPCR, analyzing both qualitative and quantitative detection. Additional findings by mNGS were confirmed by additional qPCR analysis. In case no remaining sample was available, the Galileo Analytics software results were compared with results from analysis using alternative bioinformatic tools: metagenomic taxonomic classifier Centrifuge (1.0.4-beta) [16] and de novo assembly-based viral metagenomic analysis software Genome Detective [17].

Calibration curves

After metagenomic sequencing, the viral loads were calculated for each virus by the Galileo Analytics web application. Signals of both the calibrators and patient plasma samples were plotted in load graphs (Supplementary Figure 1) and the corresponding viral load of the patient samples was extrapolated. Since no calibrator panels for B19V and TTV virus were available, these signals were plotted against other calibration curves of viruses that showed the optimal agreement with the known viral load for semi-quantitative detection. All calibration sample signals correlated well with the titer (R² range 0.84-0.92).

Viral load by mNGS versus qPCR

In total six patients were tested by qPCR and mNGS for quantification of different viruses at subsequent time points. Agreement between the methods for qualitative detection was 100% for the viruses targeted by PCR. Quantitative results per patient are shown in Table 1, and Figure 1 depicts viral loads by mNGS versus qPCR per target virus. CMV and EBV viral loads showed highest agreement, with a maximum difference in viral load of 0.70 log₁₀ IU/ml. Mean differences in viral loads were 0.43 for CMV and 0.19 log₁₀ IU/ml for EBV. For ADV, viral loads were higher when quantified with mNGS with a mean difference of 0.90 log₁₀ c/ml. For BKV, viral loads by mNGS were lower in comparison with qPCR, with a mean difference of 1.32 log₁₀ c/ml. When taking into account viral loads measured above the limit of quantification of 2.5 log₁₀ c/ml, as applied in our diagnostic qPCR for BKV, the mean difference is 0.62 log₁₀ c/ml and a trend towards better agreement with higher viral loads could be observed. Semi-quantitative detection of B19V and TTV viruses by mNGS resulted in mean differences of, respectively, 0.39 log₁₀ IU/ml and 3.0 log₁₀ c/ml in comparison with qPCR.

• Download figure
• Open in new tab

Figure 1.

Viral loads as predicted by Galileo Viral Panel mNGS versus qPCR (copies/ml for ADV, BK, and TTV, and IU/ml for CMV, EBV and B19V). B19V and TTV results were considered semi-quantitative since no Galileo calibration panels were available for these targets.

Longitudinal patient follow-up and clinical decision-making

Furthermore, for each patient the viral loads over time were plotted in graphs with clinical information about symptomatology and treatment (Figure 2). For CMV, EBV and BKV, in our clinical practice, specific viral load thresholds are used to decide whether immunosuppression should be tapered and/or antiviral therapy should be administered. Viral load quantification around these thresholds showed good agreement in identifying these clinical decision-making breakpoints. In Patient 3, antiviral treatment with foscarnet was started for CMV-reactivation when viral load measured by qPCR exceeded 4.0 log₁₀ IU/ml. By mNGS, this critical threshold for treatment initiation was correctly identified with a viral load by mNGS of 5.44 log₁₀ IU/ml. In the same patient, rituximab was administered when the EBV load by qPCR was repeatedly above the threshold of 4.0 log₁₀ IU/ml, consistently quantified thrice above 4.0 log₁₀ IU/ml before administration of rituximab both by qPCR and mNGS.

• Download figure
• Open in new tab

Figure 2.

Longitudinal follow-up of DNA viral loads in immunosuppressed patients over time, as predicted by mNGS (Galileo Viral Panel, Arc Bio) versus qPCR. Clinical information and therapeutic agents are included.

For B19V, ADV and TTV, no predefined thresholds were used for changing the treatment regimen. For all viruses, the observed trends in load over time in each patient were comparable for qPCR and mNGS, despite the semi-quantitative nature of the B19V mNGS assay. Effect of treatment (anti-viral drugs, immunoglobulins and/or tapering of immunosuppressive drugs) in patients was estimated by follow-up of viral loads by qPCR. For B19V in Patient 5 and 6, the effect of intravenous immunoglobulins (IVIG) could be assessed by the decreasing viral load in the weeks after administration, as observed by mNGS. For ADV, in patient 1, antiviral therapy with cidofovir was started when increasing viral loads were consistently detected, both by qPCR and mNGS.

Additional findings

For some samples additional viral reads were detected in the pathogenic mNGS reports, of viruses that were not initially tested for by qPCR (Supplementary Table 1). Most additional findings were supported by a secondary bioinformatic analysis using Centrifuge and Genome Detective: BK (1 patient), CMV (1 patient), HHV-6B (1 patient), and TTV (4 patients). In a few cases, additional findings were not confirmed by a second analysis leaving some low mNGS signals for CMV, EBV and HSV. JCV was detected by mNGS in a sample with a high concentration of BKV, possibly indicating forced alignment contamination due to high sequence homology between JCV and BKV [13][14].