Genetic Characteristics of Anelloviruses Detected in Individual Viromes of Children with Acute Respiratory Infection Using the Metagenomic Approach

1. Introduction

Metagenomic studies, especially in the field of virome, have shown significant growth in recent years due to their advantages investigation of viruses that are inaccessible to classical methods, characterizing viral diversity in various environments, detecting of viruses in samples with low viral load, identifying novel strains and study viral evolution. Thus, this led to detailed studies on the human respiratory virome [1,2,3]. It has been demonstrated that one of the ubiquitous members of the human respiratory virome is the enigmatic Anelloviridae family viruses [1]. It is presented in up to 90% of individuals and acquired in childhood [4,5]. Various transmission routes are described such as vertical, fecal-oral, parenteral, sexual [6,7,8]. With age, the “anellome” (the composition of anello-viruses (AVs) in the body), representing the diversity of AVs within an individual, be-comes stable; however, in childhood, the composition is variable with a predominance of TTMV (Torque teno mini virus) and TTMDV (Torque teno midi virus) [9,10].

The Anelloviridae family comprises 30 genera, four of which are detected in humans: Alphatorquevirus (Torque teno virus, TTV), Betatorquevirus (TTMV, previously known as TTV-like Mini virus, TLMV), Gammatorquevirus (TTMDV), and the recently discovered Hetorquevirus (Torque teno hominid virus) [11]. This group includes non-enveloped viruses with single-stranded negative-sense DNA circular genomes ranging in size from 1.6 to 3.9 kb [5,11]. AVs exhibit a high degree of heterogeneity among DNA viruses [9,12,13]. It is proposed to consider the species demarcation threshold as nucleotide sequence identity of ORF1 at 31% [11]. This means that if a new AV shows 69% or higher pairwise similarity with a classified member of the species, it is classified within that species.

While it is not definitively established whether AVs cause diseases, it is currently considered that they are commensal viruses. However, their involvement in the development of a wide range of pathologies, including respiratory infections, hepatitis, multiple sclerosis, lymphoma, autoimmune diseases, and others, is also suspected [3,4,6,14,15]. It is reliably noted that an increase in viral load occurs in conditions of compromised immune function (such as HIV/AIDS) [16,17], and replication is controlled by the immune system [18]. Moreover, there is a growing number of studies investigating whether AVs are associated with respiratory diseases [1,6,8,19,20,21,22,23]. In this study, we employed viral metagenomics to investigate the presence, diversity, phylogenetic relationship of AVs among pediatric patients exhibiting acute respiratory symptoms.

2. Materials and Methods

2.1. Sample Collection

The study was approved by the Ethics Committee of the Federal Research Center for Fundamental and Translational Medicine. Nasal and throat swabs were collected from 1310 pediatric patients hospitalized with respiratory symptoms (from October 2022 to May 2023) in Novosibirsk, Russia. Clinical diagnoses and demographic data of patients were taken from medical records. Samples were placed in tubes with transport medium (Dulbecco’s modified Eagle’s medium (Capricorn Scientific, Ebsdorfergrund, Germany) with 0.5% bovine serum albumin, 100 μg/mL of gentamicin sulfate (BioloT, Saint-Petersburg, Russia), and 50 units/mL of amphotericin B (BioloT, Saint-Petersburg, Russia). Tubes were stored in liquid nitrogen immediately. After transport to the laboratory, all samples were stored at -80 °C for future studies.

2.2. Virus Detection

All samples were tested for the presence of major respiratory viruses. RNA was extracted using the RIBO-sorb kit (Interlabservice, Moscow, Russia). Reverse transcription was performed using the REVERTA-L kit (Interlabservice, Moscow, Russia). The resulting cDNA was used to detect respiratory syncytial virus (HRSV); alphacoronaviruses and betacoronaviruses (HCoV); parainfluenza virus (HPIV); metapneumovirus (HMPV); rhinovirus (HRV); adenovirus (HAdV) and bocavirus (HBoV) using the AmpliSens ARVI-screen-FL kit (Interlabservice, Moscow, Russia).

Detection of influenza virus A, influenza virus B and SARS-CoV-2 was carried out in a one-step real-time PCR reaction using the AmpliPrime Influenza SARS-CoV-2/Flu(A/B/H1pdm09) kit (Nextbio, Moscow, Russia).

2.3. Bacteria Detection

Additionally, bacterial culture for opportunistic and pathogenic organisms (including Acinetobacter baumannii, Candida albicans, Candida lusitaniae, Candida parapsilosis, Enterobacter cloacae, Escherichia coli, Klebsiella oxitoca, Klebsiella pneumoniae, Moraxella catarrhalis, Neisseria spp., Pseudomonas aeruginosa, Staphylococcus aureus, Streptococcus pneumoniae, Streptococcus viridans) was performed in the clinical diagnostics laboratory of the Municipal Children’s Infectious Diseases Hospital, Novosibirsk, Russia. During the microbiological study, the primary sowing of clinical material on nutrient media was carried out, followed by the re-sowing of grown colonies of microorganisms on selective media. The species identification of microorganisms was carried out using the MALDI-TOF mass spectrometry method.

Four individual samples from children aged 2-4 years diagnosed with pneumonia or bronchitis were selected for further viral metagenomic analysis. We selected three samples in which respiratory viruses were present as a monoinfection and one sample without respiratory pathogens (for more details, refer to the Table 1).

2.4. Sample Preparation for Metagenomics

2.4.1. Enrichment of Samples with Virus-Like Particles

The sample preparation for viromic analysis was conducted in accordance to the NetoVIR protocol, with modifications [24]. Test tubes with transport medium and viscose swabs were mixed using a vortex for 3 minutes. An aliquot of each sample was then transferred to a clean tube and centrifuged for 3 min at 17,000 × g. The supernatant was filtered on a 0.8 μm filter (Sartorius Vivaclear, Göttingen, Germany) for 1 min at 17,000 × g. 20× enzyme buffer was added to the filtrate. Then benzonase and 1 μL of microccocal nuclease were added and incubated for 2 h at 37 °C. The reaction was stopped with 0.2 M EDTA and extraction began immediately.

2.4.2. RNA Extraction, Reverse Transcription and Amplification of Double-Stranded cDNA

Nucleic acids were extracted using a column-based extraction kit (Biolabmix, Novosibirsk, Russia). Complementary DNA was obtained using Sequence-Independent, Single-Primer-Amplification (SISPA) according to the protocol [25]. Briefly, first-strand cDNA was prepared using K-8N primer (5’-GACCATCTAGCGACCTCCACNNNNNNNN-3’). To do this, 50 pmol of the K-8N primer was added to the viral RNA and incubated for 2 min at 70°C. Then RT buffer, M-MuLV reverse transcriptase, and water from the M-MuLV-RH kit (Biolabmix, Novosibirsk, Russia) were added. The reaction mixture was incubated at 25 °C for 10 min, 42 °C for 60 min and 70 °C for 10 min.

To synthesize the second strand of cDNA, 10x Klenov Bufer, 20 pmol K-8N primer and 250 mM of each dNTP were added to the entire volume of the previous reaction. The reaction mixture was incubated at 95 °C for 3 min. Then cooled to 4 °C, added Klenow fragment (SibEnzyme, Novosibirsk, Russia) and incubated at 37 °C for 60 min. The reaction products were purified using Agencourt AMPure XP beads (Beckman Coulter, Brea, CA, USA).

Double-Stranded cDNA (ds-cDNA) amplification was carried out using the Encyclo Plus PCR kit (Evrogen, Moscow, Russia). To 10 μl of ds-cDNA, 5 μl of 10x Encyclo buffer, 1 μl of dNTP mix (10 mM each), 5 μl of K-primer (5’-GACCATCTAGCGACCTCCAC-3’) 100 pmol/μl, 1 μl of 50x Encyclo polymerase mix and 28 μl of water were added. The reaction mixture was incubated at 95 °C for 4 min, followed by 35 cycles at 95 °C for 15 s, 55 °C for 30 s, 72 °C for 2 min, and a final extension at 72 °C for 10 min. The reaction products were purified by Reaction mixtures DNA isolation kit (Biolabmix, Novosibirsk, Russia) and quantified through a Qubit 2.0 Fluorometer (Thermo Fisher Scientific, USA) with a Spectra Q HS Kit (Sesana, Moscow, Russia).

2.5. Library Preparation and Sequencing

Preparation of whole transcriptome libraries was carried out using the KAPA HyperPrep Kit (Roche, Switzerland). Briefly, Briefly, the size of cDNA fragments was assessed using a 4150 TapeStation automated gel electrophoresis platform (Agilent, USA). Then, the terminal sections of the cDNA fragments were repaired and the 3′ ends were adenylated. After this, the adapter was ligated and the target cDNA fragments were selected using KAPA Pure Beads (Roche, Switzerland). Next, the libraries were amplified with indexing primers and target fragments were selected using KAPA Pure Beads (Roche, Switzerland). Finally, the quality and quantity of whole-transcriptome cDNA libraries were assessed on a Qubit 4.0 (Thermo Scientific, Waltham, MA, USA) fluorometer using the dsDNA HS kit (Thermo Scientific, USA) and an automatic gel electrophoresis platform 4150 TapeStation (Agilent, Santa Clara, CA) using the HS D1000 ScreenTape kit.

Sequencing was performed on the GenoLab M platform (GeneMind, Shenzhen, China). The finished whole transcriptome libraries were diluted to a concentration of 4 nM and pooled according to the operating manual of the GenoLab M sequencing system (GeneMind, Shenzhen, China). The library pool was denatured and diluted to a final concentration of 2.8 pM, and the volume fraction of PhiX in the pool was 1%. Datasets is deposited in the Sequence Read Archive (SRA) (accession number are available upon request).

2.6. Bioinformatics Processing

Sequencing was conducted on the GenoLab M (GeneMind, Shenzhen, China) platform and approximately 13-17 million 150 bp paired-end reads per sample were achieved. Raw reads were trimmed to remove low-quality bases, ambiguous bases and adapter sequences using Trimmomatic v0.39 [26]. Trimmed reads have been decontaminated from the human genome using Bowtie2 [27]. To reduce the deduplication level, we utilized fastp v0.23.4 [28], Clumpify (from the BBMap v39.06 package [29]), and FastUniq v1.1 [30]. These reads were processed using the Virome Paired-End Reads (ViPER) pipeline (https://github.com/Matthijnssenslab/ViPER). For de novo assembly, we utilized metaSPAdes v3.15.5 [31] and filtered the assembled contigs based on coverage (≥10x) and length (≥300 nt). The annotation and classification of the assembled contigs were performed using DIAMOND v2.0.15 [32] and KronaTools v2.8.1 [33]. Contigs from Anelloviridae were filtered out and handled individually. To assemble according to the reference, the contigs of AVs were mapped against the nr/nt database BLAST to obtain the closest full genomic (or complete CDS) sequence, followed by consensus building from the libraries using BWA-MEM [34], SAMtools [35]. Sequences (> 400 nt) obtained were analyzed by BLASTn and classified by Anelloviridae genus. Seven sequences of AVs from phylogenetic analysis are available in GenBank (accession numbers are available upon request).

2.7. Phylogenetic Analysis

Phylogenetic analysis was performed based on the complete (or almost complete) amino acid ORF1 sequences of the seven AVs (five TTMVs, two TTMDVs) acquired in this study, their closest relatives based on BLASTx and other representative virus strains. MAFFT v7.520 [36] was employed for multiple sequence alignment using the mafft-linsi option and a maximum of 100 iterations separately for TTMVs and TTMDVs. Subsequently, the alignments were merged (“--merge”) and trimmed using TrimAl v1.2 with the “--gappyout” option [37]. Maximum likelihood (ML) phylogenetic trees were constructed for amino acid dataset using MEGA 11.0 under the best-fit substitution model (GTR + F + G + I) [38]. Bootstrap with 1,000 samples were performed to assess the robustness of tree topologies.

3. Results

We examined four nasal and throat swabs from children aged 2-4 years with respiratory symptoms from Novosibirsk, Russia, collected between November and December 2022. Three swabs showed the presence of respiratory pathogens, while sample E3 exhibited no detection of either bacterial or viral respiratory pathogens (Table 1). Despite the absence of major respiratory pathogens, the patient exhibited pneumonia symptoms, fever, malaise, nasal congestion, and cough.

To obtain comprehensive information on the prevalence of viruses in children with respiratory symptoms, a virome study was performed. The metagenomic analysis results were consistent with the qPCR diagnostic data.

3.1. Characteristics of the Obtained Reads and Contigs

Following sequencing, there were 13-17 million reads obtained. After trimming, deduplication, and removal of host reads, there were approximately 2 million reads remaining in samples E1, E2, and E3, and 4 million reads in sample E4 (Table 2). Through de novo assembly, we obtained 967 contigs in E1, 497 in E2, 944 in E3, and 962 in E4. Among these, 3, 1, 25, and 4 contigs belonged to AVs in each respective sample. After reference-based assembly using these contigs, we obtained 13 sequences of AVs (length 1,165-2,832 nt). In addition, 3 contigs (1 from E1, 2 from E4) remained unchanged in the attempt to assemble them based on the reference. Thus, a total of 16 unique sequences of AVs were obtained, 8 of which were > 2,000 nt (1 each in E1, E2 and 6 in E3), and 8 were < 2,000 nt (2 in E1, 2 in E3, and 4 in E4).

3.2. Prevalence of Viral Families by Samples

The prevailing number of viral reads were assigned to bacteriophages (Caudoviricetes, Crassvirales, Siphoviridae, Myoviridae, Inoviridae, Podoviridae, Herelleviridae), comprising 97-99% in PCR-positive samples (E1, E2, E4) and 89% in PCR-negative sample (E3) (Figure 1A). Among eukaryotic viruses (excluding bacteriophages), reads mapping to the family Anelloviridae predominated in sample E1 and Е4. Meanwhile in sample E2 the family Orthomyxoviridae predominated. Additionally, in the samples E1 and E4 we found reads mapping to the Human bocavirus and Pneumoviridae respectively (Figure 1B), which is corresponding to the respiratory viruses detected in the PCR analysis. It is worth mentioning that in sample E4, contigs belonging to the family Astroviridae were identified, specifically Astrovirus VA3. This finding is noteworthy as it is usually associated with gastrointestinal disorders [39,40], however there are also some studies reporting the presence of this virus in nasopharyngeal swabs from febrile children [41,42,43].

3.3. Characterization of AVs Sequences

Sequences of AVs were detected in all tested samples, with sample E3 exhibiting a significantly higher number of assembled contigs of Anelloviridae (p <0.05) compared to the others (Figure 1A, Table 2). While this abundance in E3 may be related to gender, as male typically have a richer anellome [44], the possibility that AVs may have played a role in the development of inflammation cannot be excluded, as the patient exhibited respiratory symptoms.

The overwhelming majority of AVs sequences identified across all patient samples belonged to the Betatorquevirus genus (TTMV) (11 out of 16, 68.8%). Meanwhile, three sequences (18.8%) were attributed to the Gammatorquevirus genus (TTMDV) and two sequences (12.5%) to the Alphatorquevirus genus (TTV). The higher content of TTMVs and TTMDVs in children is in line with studies reporting an increased prevalence of these AVs in younger age groups [10]. In our study, all individuals had TTMVs. However, E2 and E4 exhibited monoinfections (solely TTMV), E1 showed double coinfections (TTM + TTMV), and E3 presented with triple coinfections (TTM + TTMV + TTMDV) (Figure 2).

3.4. Variability of AVs

To establish a genetic relationship between AVs strains, a phylogenetic analysis based on the amino acid sequence of ORF1 was conducted for 7 assembled full-length or almost full-length genomes (the smallest being 2,520 nt ORF1, 94% of the complete genome) (Figure 3). The remaining 9 sequences had lengths less than 75% of ORF1, hence they were not included in the phylogenetic analysis.

All our sequences collected from pediatric patients show similarity to genomes obtained in a large viromic study of blood from febrile Tanzanian children [45] (Figure 3, evolutionary distances). These genomes display significant differences from previously well-described species in NCBI and do not cluster with them. This feature is particularly evident in the genus Gammatorquevirus (Figure 3), where the branch with reference strains is separate from the group of our sequences and those from the virome study of febrile children. This may indicate a significant gap in the detection and well description of TTMDV sequences, which could potentially fill the gap between these two branches.

Studied TTMV genomes differentiated into distant clades from each other and demonstrated a high degree of divergence. The exceptions are NODEA508_E2 and NODE99_E3, which exhibit a high degree of identity (99% amino acid similarity of ORF1, Supplementary Table S1) despite being collected from different children. Due to the highly divergent TTMV sequences, it can be hypothesized that patient E3 was infected with multiple strains of Betatorquevirus that evaded host immune system eradication. Two sequences of TTMDVs NODEA98_E3 and NODEA297_E3 are very similar (100% amino acid identity of ORF1, Supplementary Table S1) and cluster together on the phylogenetic tree.

Additionally, a plot was constructed showing the similarity of ORF1 for TTMV sequences covering over 80% of the complete genome (Figure 4). The alignment was performed against the prototype species for the TTMV genus - Torque teno mini virus 1-CBD279. The greatest distance from the prototype strain was demonstrated by NODE_28, NODE_58, NODE_98 from sample E3.

4. Discussion

Thus, the current metagenomic study provided viral profiles for children with respiratory symptoms. The analysis of the obtained data confirmed the results of PCR diagnostics for major respiratory pathogens. Additionally, it provided information on AVs in the examined pediatric patients. The research findings indicate that different AVs are presented in the virome of children with acute respiratory symptoms. The presence of viruses from Anelloviridae family in all samples is accordance with data on the prevalence of AV-infections in children with respiratory diseases [6,48].

Interestingly, in the symptomatically diseased patient without respiratory viruses, there was the highest representation of AVs in terms of both read numbers and genetic diversity. It is possible that this representation of AVs can activate immune and inflammatory responses, influencing the occurrence and maintenance of respiratory symptoms. Although a similar phenomenon has so far been well described only for TTV infection, in which TTV interacts with PAMPs (Pathogen-associated molecular patterns) receptors, modulating inflammation [1,8,19,48,49,50]

The low representation of AVs in patients with respiratory viruses (E1, E2, E4 as opposed to E3) prompts speculation about the conditions under which AVs have an opportunity for replication. It is likely that the immune system, activated to combat pathogens, suppresses the replication of AVs. After the elimination of the main pathogen, AVs increase to their initial levels, and the anellome stabilizes. Similar observations in fluctuations of AV levels have been made in patients with COVID-19, where the anellovirus load significantly decreased in the first weeks of SARS-CoV-2 infection [51].

TTVs were detected only in samples E1 and E3, possibly due to a tendency for anellome to shift toward the presence of Alphatorquevirus with age [10,52]. TTMDVs were identified only in the E3 sample, which exhibited the highest diversity of annelloviruses.

All patients exhibited respiratory symptoms and TTMVs were present in the anellome. However, contrary data can be found that the human AV TTMV is not associated with fever and severe infection [6]. Unfortunately, the widespread presence of AVs complicates the search for a definitive answer to whether this group of viruses truly worsens the clinical course of the disease.

5. Limitations

A technical negative control sample (i.e., a sample with transport medium and probe) that underwent metagenomic analysis was missing in our study. However, we compared the assembled contigs of respiratory viruses and found no identical contigs across different samples, indicating the absence of cross-contamination.