Deep sequencing applied to the analysis of viromes in patients with beta-thalassemia.

To date, blood banks apply routine diagnosis to a specific spectrum of transfusion-transmitted viruses. Even though this measure is considered highly efficient to control their transmission, the threat imposed by emerging viruses is increasing globally, which can impact transfusion safety, especially in the light of the accelerated viral discovery by novel sequencing technologies. One of the most important groups of patients, who may indicate the presence of emerging viruses in the field of blood transfusion, is the group of individuals who receive multiple transfusions due to hereditary hemoglobinopathies. It is possible that they harbor unknown or unsuspected parenterally-transmitted viruses. In order to elucidate this, nucleic acids from 30 patients with beta-thalassemia were analyzed by Illumina next-generation sequencing and bioinformatics analysis. Three major viral families: Anelloviridae, Flaviviridae and Hepadnaviridae were identified. Among them, anelloviruses were the most representative, being detected with high number of reads in all tested samples. Human Pegivirus 1 (HPgV-1, or GBV-C), Hepatitis B Virus (HBV) and Hepatitis C Virus (HCV) were also identified. HBV and HCV detection was expected due to the high seroprevalence in patients with beta thalassemia. Our results do not confirm the presence of emerging or unsuspected viruses threatening the transfusion safety at present, but can be used to actively search for viruses that threaten blood transfusion safety. We believe that the application of viral metagenomics in multiple-transfused patients is highly useful to monitor possible viral transfusion threats and for the annotation of their virome composition.

INTRODUCTION

Blood Transfusion Services worldwide apply highly sensitive diagnostic tests to detect Human Immunodeficiency virus (HIV), Hepatitis B Virus (HBV), Hepatitis C virus (HCV), syphilis, and in some countries Human T-cell lymphotropic virus (HTLV) and Chagas disease as well. Therefore, hemotherapeutic procedures including blood transfusion are considered safe. Although the transfusion risk of transmission of routinely tested viral agents has been dramatically reduced, even the most sensitive detection techniques cannot guarantee a “zero risk”. Therefore, occasional transfusion transmission of routinely tested agents is observed mainly due to circulation of novel mutation variants and/or the presence of very low viral load levels such as the ones reported in occult HBV infection , , HCV and HIV .

Apart from the risk of transfusion transmission of the routinely tested viruses, there is a risk of parenteral transmission of virtually any viral agent that will thereafter establish viremia. In this respect, emerging viruses which can threaten the safety of blood transfusion are especially important. The list of these viral agents is extensive and includes typical arboviral agents with short-termed viremia like dengue (DENV) , and Zika viruses (ZIKV) and others with prolonged viremic phases like Hepatitis E virus (HEV) , and Parvovirus type 4 (PARV4) .

Multiple-transfused patients with hereditary hemoglobinopathies are at high risk of acquiring parenterally-transmitted viruses. Historically, these patients demonstrate higher prevalences of HIV, HCV and HBV compared to the general population. Moreover, due to frequent transfusions, one may infer that unsuspected viral infections which are not routinely tested in Blood Transfusion Services can also be present.

Considering that it is impossible to test all possible infections that can threaten transfusion safety, a suitable approach would be the virome examination of multiple -transfused patients through metagenomic High Throughput Sequencing (mHTS). This approach gives a detailed information regarding the overall viral abundance and diversity in any type of clinical sample, thus speeding up the discovery of emerging viruses , . Therefore, the use of mHTS for evaluating the virome in patients receiving multiple transfusions can provide important information on the presence of unsuspected viral infections that can threaten blood transfusion.

The aim of our study was to annotate the virome of patients with beta-thalassemia major in order to search for viral agents that can threat blood transfusion safety. The choice of this group was based on the following reasons: (i) Historically, patients with beta thalassemia present with higher indexes of transfusion-transmitted diseases, which have been reduced nowadays, but continue elevated compared to the general population. This demonstrates that patients with beta-thalassemia are under a constant threat of acquiring parenterally-transmitted infections, especially the emerging ones. (ii) Presence of endemic viral diseases in the examined region like dengue and Zika infections that are transmitted by arthropod vectors but might also be transmitted by blood transfusion.

MATERIALS AND METHODS

Study population

In the present study, we evaluated the virome of 30 patients with beta-thalassemia, who sought medical attention at the hemoglobinopathies outpatient clinic of the Blood Center of Ribeirao Preto, Faculty of Medicine of Ribeirao Preto, University of Sao Paulo. From these patients, 20 were male (66.7 %) and 10 were female (33.3%) with a mean age of 26.6 years (range 6-47 years). The serological screening for the routinely tested viruses transmitted by blood transfusion revealed a seroprevalence of 23.3% (n=7/30) for anti-HCV IgG and 16.7% of anti-Hbs Ag (n=7/30).

Between September-October, 2018, blood samples were collected individually from each patient during monthly examinations at the outpatient clinic. All patients included in this study were transfusion-dependent and represented a total of 73.2% of all thalassemia patients treated in our institution. Approximately 6 mL of whole blood were collected in a sterile vacutainer tube (Beckton-Dickinson, NJ, USA), after individual informed consent. The samples were initially centrifuged at low speed for plasma separation (1,931 × g, 10 min) and if not immediately processed they were stored at −80 °C. The study was approved by the Institutional Ethics Committee of the University Hospital at the Faculty of Medicine of Ribeirao Preto, University of Sao Paulo, Brazil (process N° HCRP-12196/2018).

Sample preparation and Illumina sequencing

Initially, 600 μL of plasma obtained from each patient were pre-treated with Turbo DNase (Turbo DNA-free kit, Thermo Fisher Scientific, Waltham, MA, USA) in order to remove free host/bacterial DNA in the analyzed samples (DNAse incubation at 37 °C for 30 min). After DNase inactivation, five individual samples were mixed in a single pool. The total pool volume (3 mL) was extracted using the High Pure Viral Nucleic Acid Large Volume Kit (Roche, Sao Paulo, Brazil) following the manufacturer's instructions with minor modifications: the polyA carrier RNA from the kit was replaced by the neutral carrier i.e. GenElute Linear Polyacrylamide carrier (LPA) (Merck & Co., Kenilworth, New Jersey, USA) and Isopropanol was added to the Binding buffer before the filtration step. After column extraction, viral nucleic acids were submitted to reverse transcription using the Superscript III First-Strand Synthesis System (ThermoFisher Scientific, Waltham, MA, USA). Amplification of synthesized cDNA was performed using the QuantiTect Whole Transcriptome Kit (QIAGEN, Hilden, Germany). The preparation of libraries for sequencing was performed using the Nextera DNA Flex Library Preparation Kit (Illumina, San Diego, CA, USA) and the Nextera DNA CD Indexes kit following the manufacturer's instructions. The sequencing the dual-indexed libraries was carried out with an Illumina NextSeq 550 sequencer using the NextSeq High Output Kit v.2.5, 300 cycles (Illumina, San Diego, CA, USA) following the manufacturer's instructions. The sensitivity of viral identification sequencing was evaluated using pools containing clinically relevant DNA and RNA viral agents with varying levels of viral load, as described previously .

Bioinformatics processing of the obtained sequencing data

Sequence reads were initially processed for quality examination using FastQC (v. 0.11.8, Babraham Institute, Cambridge, UK). The low-quality reads, sequences with ambiguous bases and duplicated reads were removed using PRINSEQ (v. 0.20.4, SDSU, San Diego, CA, USA). To remove any host-related sequences and to retain as much viral sequences as possible, we used DeconSeq (v. 0.4.3, SDSU, San Diego, CA, USA). The remaining reads were taxonomically analyzed using Kraken (v. 2, Johns Hopkins University, Baltimore, MD, USA). The unclassified sequences were analyzed using Blastn and Blastx. Assembly was performed by SPAdes (v. 3.13.0, St. Petersburg State University, CAB, St. Petersburg, Russia). All generated contigs were analyzed using Blastn with nucleotide database obtained by the NCBI and an E-value cutoff of 10-3. For protein similarity evaluation, we used DIAMOND (v.2.0.6, University of Tübingen, Tübingen, Germany). The presence of environmental or other contaminants/possible misclassifications was evaluated using the P-DiP pipeline as previously established .

Phylogenetic analysis of the most important viral contigs

Phylogenetic analysis of the consensus contigs was performed using datasets of complete or partial genomes obtained from NCBI (National Center for Biotechnology Information). Multiple alignment was performed using MAFFT (v.7.450, Kyoto University, Kyoto, Japan) and the phylogenetic signal was verified using TREE-PUZZLE v. 5.3 (Max-Planck-Institute, Berlin, Germany). The reconstruction of the trees was performed with IQ-TREE (v. 1.6.8, CIBV, Vienna, Austria) applying the Maximum Likelihood method with a statistical support of 10,000 bootstrap replicates. The phylogenetic trees were visualized with FigTree (v. 1.4, IEB, University of Edinburgh, Edinburgh, UK). Only bootstrap values of above 75% were considered statistically significant and were represented in important phylogenetic branches.

RESULTS

The Illumina sequencing generated a total of 808,427,182 reads in plasma samples obtained from patients with beta thalassemia. We determined an average of 134 million reads with nucleotide lengths of 150 bp for each plasma pool (range: 119-160×106 reads/pool) (Table 1). After trimming, an average of 13.5 x 106 reads (10.48%) were eliminated from five pools, while one pool (number 6) showed removal of 82 x 106 million reads (50.9% from its total reads). After taxonomic classification, the classified sequences (range: 73-119 x 106 reads/pool) were kept for de novo assembly and Blastn/Blastx analysis. On average, viral sequences corresponded to 0.023 % of reads per pool (range: 0.002 to 0.11%). The unclassified sequences corresponded to a medium of 5.82 % of the total reads (range: 5.19% - 6.24%). From the classified viral sequences, the maximum contig length was 3,786 bp and belonged to HCV (start position 3,441 bp ends 7,158 bp when aligned with the reference HCV genome NC_038882. This region encompassed the entire NS3, NS4A and NS4B genes, and partially to NS5A gene). We observed that ~ 93% of the classified contigs were assigned as human sequences. From the known viruses, we detected three major viral families: Anelloviridae, Flaviviridae and Hepadnaviridae. Most of the viral reads corresponded to anelloviruses and were mainly assigned to Torque teno virus (TTV) (> 95%) from the Alphatorquetenovirus genus, detected in each pool. Other Anelloviridae family genus including Gammatorquevirus (Torque teno midi virus, TTMDV) and Betatorquevirus (Torque teno mini virus, TTMV), were also represented (1.5-16.7% of all anelloviruses present). Human pegivirus-1 (HPgV-1, formerly known as GB virus C or Hepatitis G virus) was also found in three out of six plasma pools with a particularly high number of reads in pool number 2 (158,790 reads) (Table 2). The performed phylogenetic reconstruction of the largest HPgV-1 contig (1,360 bp, sequence starting at position 2,228 bp and ending at position at 3,554 bp when aligned with HPgV-1 genome NC_001710 corresponding to the p7-NS2 gene), demonstrated that it belongs to genotype 2 (subgenotype 2A) (Figure 1).

Table 1

Workflow of the sequences analysis, from raw data to assignment.

Pool number	Total reads	Reads after trimming	Classified reads	Viral reads	Viral Reads (%)	Unclassified reads
Pool 1	126,241,206	112,680,449	106,827,862	2,642	0.002%	5,852,587
Pool 2	140,189,184	120,822,545	113,278,558	141,841	0.11%	7,543,987
Pool 3	123,039,880	109,175,346	102,493,334	26,556	0.02%	6,682,012
Pool 4	119,051,092	107,960,664	101,267,675	3,414	0.003%	6,692,989
Pool 5	139,336,406	128,896,726	119,770,867	3,673	0.003%	7,125,859
Pool 6	160,569,414	78,111,788	73,753,074	3,555	0.002%	4,358,714

Table 2

Number of reads of the main viral families detected by the applied metagenomic pipeline in plasma samples from patients with beta-thalassemia.

Samples	Pool 1	Pool 2	Pool 3	Pool 4	Pool 5	Pool 6
Anelloviridae	1,794	13,126	42,954	4,480	4,081	3,088
Alphatorquevirus	1,724	11,384	41,587	3,912	3,516	3,012
Gamatorquevirus	52	1,232	827	408	313	23
Betatorquevirus	5	251	82	42	132	17
HPgV-1 (GBV-C)	-	158,790	5	26	-	-
Hepatitis B	-	-	-	8	-	-
Hepatitis C	-	-	-	-	-	1,390

Figure 1

Phylogenetic trees of the most abundant viruses detected by viral metagenomic analysis in patients with β-thalassemia major. The phylogenetic trees were reconstructed using the General Time Reversible (GTR) and Transversion Nucleotide Substitution model (TVM) selected by the ModelFinder through Bayesian Information Criteria (BIC), as implemented by IQtree software v. 1.6.2. The Maximum-Likelihood method was used for tree topology reconstruction and for statistical support we used 10,000 ultrafast bootstrap replicates. Only bootstrap values above 75% were maintained at important tree branches: A) Phylogenetic tree of the contig of Human Pegivirus-1, which was classified as subgenotype 2A of the main genotype 2 that is commonly found in patients with high risk of parenteral exposure. The GTR+F+R4 nucleotide substitution model was used for tree reconstruction; B) Taxonomic classification of the HBV contig, which clustered within genotype A that is widely spread in Brazil. For tree reconstruction the TVM+F+I+G4 nucleotide substitution model was used; C) Phylogenetic tree of HCV. The arrow shows the position of the assembled contig within the subgenotype 1A of the main HCV genotype 1. For the tree reconstruction, the GTR+F+R5 nucleotide substitution model was used.

We have also detected a large number of sequences (1,390 reads) in one of the pools (pool 6, Table 2) corresponding to HCV (family Flaviviridae). The detection of HCV nucleic acids was related to the presence of HCV viremia in one of the samples of the pool (viral load of HCV, 56,793 IU/mL). The performed phylogenetic analysis of the obtained 3,786 bp contigs demonstrated that it belongs to HCV genotype 1A (Figure 1). The established bioinformatics pipeline was also able to detect 8 reads belonging to HBV (Hepadnaviridae family) in one pool (pool 4). The phylogenetic analysis of the obtained HBV contigs (244 bp, sequence starting at position 1,610 bp and ending at position 1,853 bp aligned with reference genome NC_003977.2 and corresponding partially to the genomic region of the X protein) demonstrated that it belongs to genotype A. The consensus sequences used in the phylogenetic analyses were deposited under the following numbers in the NCBI i.e., HCV (MT596912), HPgV-1 (MT596914) and HBV (MT596913).

DISCUSSION

In this study, performed among chronically transfused Brazilian patients with beta-thalassemia, we evaluated their virome by mHTS. By applying this approach, we tried to identify emerging or unsuspected viral agents which can threaten the transfusion safety. The main reasons to perform this study are the high-risk of the thalassemic patients to acquire parenteral infections and the endemicity of the studied region for arboviruses and other tropical infections which can be transmitted via blood transfusion , .

In all of the screened pools we found viruses showing persistent viremia, thus demonstrating that transmission can occur parenterally. Among them, the most abundant were the members of the Anelloviridae family, represented by the three major genera: Alpha-, Beta- and Gammatorquetenvirus. By far, this is the most abundant group of viruses found in blood samples analyzed by viral metagenomics – . The obtained results are in agreement with the fact that anelloviruses are characterized by a high level of genome diversity and high prevalences in the human population which in some cases can reach 100% , . Despite the high anellovirus prevalence, the way these viruses interact with the hosts is not elucidated.

Another widely distributed virus, which was found by the metagenome pipeline was the HPgV-1, showing an elevated number of reads in one of the tested pools. Further phylogenetic analysis of the major contigs confirmed the presence of HPgV-1 genotype 2 (subgenotype 2A), which is the most prevalent genotype in Brazil , that is efficiently transmitted by blood transfusion. Similar to the anelloviruses, HPgV-1 is highly prevalent in different patient groups, especially in patients with high parenteral exposure – . This is in accordance with our findings of a high number of HPgV-1 contigs detected among patients with beta-thalassemia. Our results are similar to data from the literature which demonstrated that HPgV-1 and anelloviruses are the most commonly found commensal viruses detected in human plasma . In this respect, these commensal viruses were also most frequently identified by viral metagenomic approaches in blood donations obtained from healthy donors , . In general, healthy blood donors are a suitable control group that can also be used to compare our results in order to evaluate the impact of these commensal viruses on blood transfusion. Although we did not implement such a control, we believe that the composition of the commensal virome between patients and donors might be quite similar based on similarities obtained from the literature data , .

On the other hand, important parenterally-transmitted viruses like HCV and HBV were also identified which differs from the normal blood donor virome. The detection of HCV reads in one of the tested pools was due to the high seroprevalence of anti-HCV IgG in the tested group of patients (23.3%). The phylogenetic analysis of the largest HCV contig (3,786 bp) demonstrated that it belongs to genotype 1A. This finding is in accordance with the literature, pointing out that genotype 1A is the most prevalent HCV genotype in Brazil . We also detected a few reads belonging to HBV. The HBV contig was composed of 244 bp and although relatively small, it was further classified as belonging to HBV genotype A, which is also a common Brazilian genotype . The identification of few HBV reads probably reflects low viral loads and that the sensitivity of the applied sequencing is adequate for the identification of low viral load infections. This was supported on one hand by our previous studies which sensitively identify underrepresented viruses in artificially spiked pools , and on the other by the extraordinary deepness of the applied sequencing (around 200 million reads for each tested pool) which possibly enabled the detection of infections with low viral loads. In addition, the obtained results show that the presence of HCV-RNA and HBV-DNA in multiple - transfused patients with beta-thalassemia major is in accordance with the genotype distribution of these two viruses in Brazil, and that the used HTS pipeline was capable of determining the presence of active viremia of chronic viral infections.

A shortcoming was that annotated viral sequences were not confirmed by direct molecular methods like PCR. However, metagenomics has a high diagnostic potential due to the correlation between the percentage of viral reads of a given virus and its viral load in the examined sample , . Nevertheless, sequence reads of only highly represented agents were taken into account, except for HBV, which was previously confirmed by serological analysis. According to some authors, metagenomics could be considered even superior to direct molecular methods when working with highly diverse viral taxa such as the Anellovirus family for which exact primer matching and efficient amplification are difficult to be achieved – .

In this study, we did not characterize the sequence of any new virus that can be regarded as an emerging threat to multiple-transfused patients with beta-thalassemia. In general, the identification of emerging viruses is dependent on nucleotide alignments with known viruses and it is difficult to taxonomically classify phylogenetically distant viral agents . However, we identified some sequences with low nucleotide sequence coverage (20-30%) belonging to the Marseilleviridae, Phycodnaviridae and Mimiviridae families which represent giant viruses some of which have been putatively transmitted by blood transfusion . Due to the low percentage of genomic coverage, we could not attribute with certainty that these viral agents are emerging viruses threatening the transfusion policies in Brazil. Therefore, we believe that there is a need for improvement of the current bioinformatics approach to be able to classify and identify sequences of phylogenetically distant viruses .

CONCLUSION

In conclusion, this study examined the viral landscape in patients with beta-thalassemia. We detected only known commensal and parenterally-transmitted viruses belonging to different viral families in high-risk patients with beta-thalassemia major. Our study demonstrates that viral metagenomic approaches represent an interesting tool for the evaluation of viral diversity in patients subjected to multiple transfusions, reinforcing its relevance for the discovery of viral agents that can potentially threaten the transfusion safety.