Nasopharyngeal virome analysis of COVID-19 patients during three different waves in Campania region of Italy.

From December 2019, severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) infection has spread rapidly, leading to a global pandemic. Little is known about possible relationships between SARS-CoV-2 and other viruses in the respiratory system affecting patient prognosis and outcomes. This study aims to characterize respiratory virome profiles in association with SARS-CoV-2 infection and disease severity, through the analysis in 89 nasopharyngeal swabs collected in a patient's cohort from the Campania region (Southern Italy). Results show coinfections with viral species belonging to Coronaviridae, Retroviridae, Herpesviridae, Poxviridae, Pneumoviridae, Pandoraviridae, and Anelloviridae families and only 2% of the cases (2/89) identified respiratory viruses.

INTRODUCTION

In December 2019, a few cases of “pneumonia of unknown aetiology” were reported in Wuhan (Hubei region, China). On January 7, 2020, a new coronavirus (CoV), highly related to bats’ SARS‐like virus, was isolated and identified. Not even a decade since the Middle East respiratory syndrome‐related CoV and 15 years after the severe acute respiratory syndrome coronavirus 1 (SARS‐CoV‐1) epidemic, another menace, posing a new unimaginable challenge.

Since the advent of this pandemic, in October 2021, more than 240 million confirmed cases and over 4.9 million deaths have been reported (World Health Organization data), and a variable case‐fatality rate estimated to be slightly below 3%, influenced by several factors such as patient's age and comorbidities, healthcare setting, geography, and epidemic phase.

Most of the infected individuals have mild clinical symptoms, while only about 20% of positive patients progress to severe disease. During COVID‐19 treatment, many factors have been shown to affect patient prognosis and one of those is respiratory coinfections with other viruses. The COVID‐19 outbreak occurred first during the winter season, with a high incidence of other respiratory viruses such as influenza viruses. Indeed, SARS‐CoV‐2 coinfections with entero/rhinovirus, human metapneumovirus, respiratory syncytial virus (RSV), other coronaviruses (non‐SARS‐CoV‐2), and influenza A virus have been reported in several studies. , In these situations, infected patients are more prone to develop acute respiratory distress syndrome, which makes their management more challenging. Still, little is known about the impact of coinfections in COVID‐19 patients and a proper virome profiling of SARS‐CoV‐2 infection sites is lacking.

Our study is a descriptive analysis aiming at the identification of SARS‐CoV‐2 and other viral coinfections, to assess the possible association of such coinfections with disease severity. Here, nasopharyngeal swabs samples were collected from a Campania region (Southern Italy) cohort of 89 patients that have been diagnosed with COVID‐19. If and how SARS‐CoV‐2 infections can influence the composition of the upper respiratory tract remains unclear. Therefore, the differentiation between SARS‐CoV‐2 single‐infection and SARS‐CoV‐2 coinfection with other pathogens, and in particular other viruses, is of huge importance for clinical management.

METHODS

Samples cohort

The cohort of SARS‐CoV‐2 infection cases from the Campania region (Southern Italy) consists of 89 patients. Nasopharyngeal swabs were collected during the three main COVID‐19 waves in Italy: first wave (March–May 2020); second wave (September–November 2020); and third wave (January–February 2021). Infections were then confirmed through a positive molecular test. In total, 27 positive cases from the first period were included in this study, as well as 43 from the second and 19 from the third period. Forty‐six percent of patients were female (n  =  41) and 54% were male (n  =  48) with a median age (interquartile range) of 55 years, ranging from 3 to 99 years (ethical approved number 1316, November 23, 2020). Patients were distributed on the basis of symptom severity as previously described into nonsevere (total: n = 49; asymptomatic: n = 26; and mild: n = 13 cases), moderate (n = 6), severe (n = 10 included 3 deceased), and unknown groups (n = 34). Patient data are summarized in Table 1.

Table 1

Epidemiological features of the 89 cohort members between the three collection periods

	Mar–May 2020 (n = 27)	Sep–Nov 2020 (n = 43)	Jan–Feb 2021 (n = 19)
Age (years)
0–20	3 (11%)	3 (7%)	5 (26%)
21–40	4 (15%)	8 (19%)	1 (5%)
41–60	4 (15%)	15 (35%)	–
61–80	10 (37%)	13 (30%)	6 (32%)
>80	3 (11%)	4 (9%)	7 (37%)
Unknown	3 (11%)	–
Gender
Male (%)	16 (59%)	27 (63%)	5 (26%)
Female (%)	8 (30%)	16 (37%)	14 (74%)
Unknown	3 (11%)	–	–
Disease severity (%)
Asymptomatic	11 (41%)	15 (35%)	–
Mild	7 (26%)	6 (14%)	–
Moderate	2 (7%)	–	4 (21%)
Severe	5 (2 dead) (19%)	4 (1 dead) (9%)	1 (5%)
Unknown	2 (7%)	18 (42%)	14 (74%)

Note: The values shown in this table are expressed in the format of number (percentage).

Library preparation, sequencing, and bioinformatics analysis

RNA was extracted from 200 µl of 89 nasopharyngeal swabs using ELITeInGenius fully automated system (ELITechGroup) and ELITeInGenius SP RNA cartridge (ELITechGroup), which exploits a magnetic bead technology, eluting in 100 µl. Extracted RNAs were retro‐transcribed using SensiFAST™ cDNA Synthesis Kit (meridian BIOSCIENCE). The viral load of each sample was assessed by real‐time polymerase chain reaction (RT‐PCR), targeting the Sars‐CoV‐2 viral nucleoprotein gene (forward primer: GGGGAACTTCTCCTGCTAGAAT, reverse primer: CAGACATTTTGCTCTCAAGCTG). RNA concentration was quantified using a Qubit RNA HS Assay Kit (Thermo Fisher Scientific). Libraries were made starting from 100 ng of RNA extract and using the TruSeq Stranded Total RNA Kit (Illumina) according to the manufacturer's guidelines. Briefly, RNA was depleted for ribosomal RNA, fragmented, and first‐strand complementary DNA was synthesized. The following synthesis of the second strand was performed using dUTPs instead of dTTP to quench the amplification of the second strand during the PCR amplification step. After adenylation of double‐strand DNA (dsDNA) fragments, indexed adapters were ligated and DNA fragments containing adapter molecules were enriched by 15 cycles of PCR. Final library concentration was assessed using the Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific), while library size was verified by Agilent 4200 Tapestation System (Agilent), showing an average size of 400 bp. Equimolar pools of the samples were prepared and sequenced on the NextSeq 500 (Illumina) in 2 × 75paired‐end mode at a final concentration of 1.7 pMol or on NovaSeq 6000 (Illumina) in 2 × 100 bp paired‐end mode at a final concentration of 250 pMol. The sequencing runs generated 57.6 Gbp of data, with 91.8% of passing filter reads and 94.5% of reads with a quality ≥Q30 for NextSeq and 880.7 Gbp of data, with 82.45% of passing filter reads and 93.3% of reads with a quality ≥Q30 for NovaSeq. Raw sequencing data were analyzed with the HOME‐BIO pipeline. Host‐related sequences were filtered out by mapping on the human reference genome (GRCh38.p13 release 37) and viral taxonomy assignation was obtained with default parameters by querying RefSeq complete viral genomes/proteins database. Classification data were then imported in R software (version 3.6.3) and normalized in reads per million (RPM) values (RPM mapped on the viral database).

RESULTS AND DISCUSSION

A total of 9.7 billion raw reads were obtained, with an average of 109.4 million reads per sample (range 26 853.188–240 193.698 reads). For the entire dataset, 609 million reads were mapped on the virus database, with an average of 6.8 million viral reads per sample (range 1718–52 613.608 reads). Sequences related to viruses and their targeted natural hosts were identified. Those specific hosts include invertebrates, plants, fungi, protozoa, and bacteria. For further analysis, reads related to bacteria and phages have not been considered in this study. We focused our attention on eukaryotic viruses, for which the reads relative abundance per sample varied from a minimum of 0.01% of the total viral reads.

As expected, SARS‐CoV‐2 (Coronaviridae family), is the most abundant viral species identified with an average RPM of 703 555 (range 1582–993 313 RPM). In addition, six other viral families were detected during the three different waves: Retroviridae, Herpesviridae, Poxviridae, Pneumoviridae, Pandoraviridae, and Anelloviridae (Figure 1A–C and Table 2). The Retroviridae family was identified in 76% (68/89) of the samples, 6 females and 14 males belonging to the first wave, 10 females as well as 22 males from the second, and 8 females and 5 males from the third (Figure 1D). No direct association with the disease severity seems to be revealed. Amongst those 68 patients, 1 was infected by Lentivirus of human immunodeficiency virus‐1 species, and the others by human endogenous retroviruses K (HERV‐K) species. The patient positive for human immunodeficiency virus‐1 was a female infected by Sars‐CoV‐2 during the first wave, which outcome was fatal. It was previously shown that, in HIV patients, the mortality associated with COVID‐19 disease is higher. Regarding the patients infected by HERV‐K, 10 were characterized by severe outcomes (Figure 1D). Souza et al.  recently reported (preprint version) that the presence of Retroviridae HERV‐K in the lower respiratory tract of severe COVID‐19 patients is associated with early mortality. These results are pointing out its possible association with disease severity.

Figure 1

Mean of RPM reads related to detected families among first (A), second (B), and third waves (C). Values are reported as the percentage of all RPM assigned to detected families in the considered period. (D) Heatmap reporting Log10 RPM values of species on entire dataset. RPM, reads per millionl SARS‐COV‐2, severe acute respiratory syndrome coronavirus 2

Table 2

Overview of the most abundant viruses (family and genus) detected in SARS‐COV‐2 positive samples from three different waves in the Campania Region

Mar–May 2020 (10 genus from 7 families)		Sep–Nov 2020 (6 genus from 4 families)		Jan–Feb 2021 (5 genus from 5 families)
Family (genus)	Reads number (RPM)	Family (genus)	Reads number (RPM)	Family (genus)	Reads number (RPM)
Coronaviridae (Betacoronavirus)	517.032 516.869	Coronaviridae (Betacoronavirus)	765.266 764.777	Coronaviridae (Betacoronavirus)	867.010 866.518
Pneumoviridae (Orthopneumovirus)	694 1.225	–	–	–	–
Retroviridae (Lentivirus, human immunodeficiency virus‐1)	71.339 86	Retroviridae Human endogenous retroviruses	6.532 6.529	Retroviridae Human endogenous retroviruses	6.010 6.009
Human endogenous retroviruses	71.233
Herpesviridae (Simplexvirus) (Lymphocryptovirus)	19.160 18.539 138	Herpesviridae (Simplexvirus) (Roseolovirus) (Lymphocryptovirus)	6.704 6.623 8 36	Herpesviridae (Simplexvirus)	849 815
Poxviridae Chordopoxvirinae (subfamily) Orthopoxvirus BeAn58058 virus	7.793 7.793	Poxviridae Orthopoxvirus (BeAn58058virus)	680 680	Poxviridae Orthopoxvirus (BeAn58058virus)	543 543
Anelloviridae (Alphatorquevirus)	28 28	–	–	Anelloviridae Alphatorquevirus	14 14
Pandoraviridae (Pandoravirus)	71	–	–	–	–

Note: The viral reads (expressed as a mean of the reads/total samples for each period) identified from samples collected during the first, second, and third waves were classified into 7, 4, and 5 families, respectively.

Abbreviations: RPM, reads per million; SARS‐COV‐2, severe acute respiratory syndrome coronavirus 2.

Another highly represented viral family was the Herpesviridae, which was detected in 21% of the patients (20/89). The specific species found included human‐alfa‐herpesvirus 1 in 9% (8/89) of our samples, human gamma‐herpesvirus 4 (Epstein–Barr virus) in 4.5% (4/89) and human‐betaherpesvirus 6A as well as human‐betaherpesvirus 7A found in, respectively, in 1% and 2.2% (1/89 and 2/89).

In particular, human‐alfa‐herpesvirus 1 was found in five males and in three females distributed along the three different waves. Human gamma‐herpesvirus 4 was found in two males and two females from the first and second waves, while human‐betaherpesvirus 6A was identified in a male from the second wave. Human‐betaherpesvirus 7A was detected in a male and a female during the first sampling campaign. In our data, the human gamma‐herpesvirus 4 was detected in patients with mild to severe/deadly outcomes (Figure 1D).

Human‐alfa‐herpesvirus 1 was also found in samples with severity ranging from nonsevere (n = 1), moderate (n = 2), to severe (n = 5) (Figure 1D). Consequently, it seems that those two species are linked with poorer outcomes in our cohort of COVID‐19 patients. Interestingly, Katz et al. observed that human‐alfa‐herpesvirus 1 reactivation occurs more frequently in COVID‐19 patients than in the normal population. Also, it was already noted, in previous studies, that the human‐alfa‐herpesvirus 1 presence in the lungs of pneumonia patients was associated with worse outcomes. , A reason for Herpesviridae detection in the more severe COVID‐19 cases might be SARS‐CoV‐2 advanced infection association with immunosuppression.

Amongst reads attributed to Poxviridae family, BeAn 58058 virus was detected in 32.5% (29/89) of our subjects. Positive samples for the BeAn virus are spread over the three waves (Figure 1D). BeAn 58058 is a zoonotic orthopoxvirus able to infect a wide range of hosts, both wild and domestic animals as well as humans. BeAn 58058 has previously been identified in postmortem Covid‐19 patients as a frequently nonpathogenic detected species.

Reads matching to the Anelloviridae family were detected as well and belonged to the Alphatorquevirus genus. In our samples, Torquetenovirus18 was found in a 76‐year‐old female patient, from the first wave, with a nonsevere (mild) outcome (Table 3) and in a 69‐year‐old male (third wave, moderate outcome, Table S2). This was codetected with Herpesviridae species (Figure 1D). It is noteworthy that, even though anelloviruses are not known to be pathogenic, they are considered possible markers of immunosuppression. In our study, Anelloviridae reads were detected at low abundance only in two patients with nonsevere (mild) and moderate outcomes. This result is most likely related to a technical limitation of RNA‐seq. Indeed, it has been shown that RNA sequencing represents a challenge for detection and quantitation of DNA virus, such as the Anelloviridae family, in biological samples as this method was not specifically designed for genomes with such complexity. Furthermore, due to the intrinsic design of the RNA‐seq assays, the low abundance of detected reads relates more to low viral RNA expression levels than to the poor representativeness of these DNA viruses in the samples.

Table 3

Metagenomic detection of viruses from human nasal‐throat swab samples SARS‐CoV‐2 positive

Sample ID	Clinical outcome	C t value	SARS‐CoV‐2 Reads	Other virus detected (reads)	Sample ID	Clinical outcome	C t value	SARS‐CoV‐2 Reads	Other virus detected (reads)
3_CA44	Asymptomatic	33.52	1.186	HERV‐ K113 (105179)	SA49	Asymptomatic	27.23	128.385	HERV‐K113 (58450)
4‐CA04	Mild	27.62	27.068	HERV‐ K113 (607115) Human beta herpesvirus7 (773) Human gamma herpesvirus4 (2320) BeAn58058virus (107501) Torquetenovirus18 (773) Pandoravirus (773)	SA56	Asymptomatic	25.90	122.347	HERV‐K113 (40449) BeAn 58058 virus (8426)
5_COS41	Dead	29.35	53.429	HERV‐ K113 (127591)	SA04	Asymptomatic	21.2	751.085	Respiratory syncytial virus (1440) HERV‐K113 (515)
6‐E6‐NA	Severe	16.79	17.4587	HERV‐ K113 (12785) Human alpha herpesvirus1 (310554) Human gamma herpesvirus4 (969)	7‐H6‐NA	Severe	24.75	553.888	HERV‐K113 (147339) Human alpha herpesvirus1 (909)
SA47	Asymptomatic	22.81	71.9161	HERV‐K113 (5650)	SA06	Severe	21.685	682.398	HERV‐K113 (35971)
8‐BE14	Dead	34.96	3.492	HERV‐K113 (562281) Human immunodeficiency virus 1 (2328) Human alpha herpesvirus 1 (15133) Pandoravirus (1164)	SA12	Moderate	27.24	11.012	HERV‐K113 (25860)
SA68	Asymptomatic	25.48	104.904	HERV‐K113 (53134) BeAn 58058 virus (5136)	p34‐A7	Mild	15.45	55.620	HERV‐K113 (24953)
SA73	Asymptomatic	22.94	704.533	HERV‐K113 (18507) Human alpha herpesvirus 1 (156540)	SA16	Asymptomatic	26.33	54.533	Human coronavirus HKU1 (61377) HERV‐K113 (92600) BeAn 58058 virus (31437)

Pandoravirus genus, and more specifically the Pandoravirus neocaledonia species, was found in two females from the first wave (Figure 1D). One of those patients had a fatal outcome and human alpha‐herpesvirus 1 was codetected (Figure 1D). The other female, the same one that presented the Torquetenovirus18, had a nonsevere outcome (mild) and presented human gamma‐herpesvirus 4 sequences. In both patients of them were codetected a high number of reads matched the HERV‐K as well as the BeAn 58058 virus (Figure 1D). Pandoraviruses are typical giant viruses of amoebas and are often detected in environmental samples, insects, and simian bushmeat. However, data showed that these giant viruses are present in humans as well when looking into various body parts of both healthy and sick individuals. This kind of virus is found in intensive care units, in patients suffering from pneumonia, and seems to be associated with ventilator use. We, unfortunately, do not know if the patient with a fatal outcome was indeed ventilated.

The Pneumoviridae family reads were found in a young (asymptomatic) child from the first wave (Figure 1D and Table 3). More specifically, they matched on the RSV species. RSV is known to cause bronchiolitis and lower respiratory tract infection in young children that can rarely progress into pneumonia.

In addition, we also detected the HKU1 coronavirus (Coronaviridae family) in an 86‐year‐old female from the first wave (Figure 1D and Table 3). Human coronaviruses, such as HKU1, generally cause mild upper‐respiratory tract illness and are responsible for common colds in human adults, however severe lower respiratory tract infections can sometimes occur in elderly people, infants, or immunocompromised patients. It is known to coinfect patients with other respiratory viruses, including other Coronaviridae pathogens.

Starting from our descriptive analysis, respiratory viral coinfections seem to be not closely associated with SARS‐CoV‐2 infection or disease severity, period of diagnosis, and gender.

Unlike several other studies reporting influenza virus coinfection with SARS‐CoV‐2, we found no presence of these viruses infection in our samples (in line with the evidence of an unusual global low circulation of influenza viruses during the pandemic period). In the patient cohort, described here, besides asymptomatic RSV detection in an 8‐year‐old child, we also found HKU1 in an elderly female. Both were diagnosed during the first wave before the Italian measures to wear face masks (introduced on April 26, 2020), social distancing, and other measures intended to stop disease spread.

A study conducted by the Icahn School at Mount Sinai, New York, reports that coinfection with other respiratory viruses appears to be rare in patients with SARS‐CoV‐2 infection. Some viruses, such as rhinoviruses, are known to interfere with the ability of other viruses to establish an infection. Hence, in our samples, the lack, or a low presence of other respiratory viruses, can be analyzed in light of these studies. Different mechanisms of the interference have been suggested, including alteration of cell surface viral receptor, cell death, or the host interferon responses. The protective antibody‐driven interferences have also been proposed for the conflict of genetically close viruses such as parainfluenza, metapneumovirus, and RSV. The immune response can be triggered by a virus, through different mechanisms, and their interactions can determine an advantage concerning competition between coinfecting viruses. From these considerations, we can speculate that competitive advantage may play a role in SARS‐CoV‐2 interaction with other respiratory viruses during coinfection, and thus could be one of the reasons why the coinfection rate in SARS‐CoV‐2 patients we analyzed is low. Factors other than viral interference could determine low virus co‐detection rates, such as variations in virus seasonality based on environmental factors and or differences in virus‐host range (e.g., range of cell types, viruses preferentially infect different age groups).

Interestingly, as others have reported before, there seems to be an association between Herpesvididae and SARS‐CoV‐2 infections , (Figure 1D). Nonetheless, in our data, the samples presenting human‐alfa‐herpesvirus 1 and human gamma‐herpesvirus 4 are mostly patients with severe outcomes and is in agreement with other observations. Additionally, we detected the contemporary presence, in two different SARS‐CoV‐2's positive samples (Figure 1D, Table 3, and Table S2), of members of the Herpesviridae and Anelloviridae family. This kind of coinfection is considered worthy of study by the scientific literature. In fact, Mallet et al. analyzed in a recent paper the association of virological markers, as the presence of herpesvirus and anellovirus with clinical outcomes and various immunological parameters to better define the causes and consequences of viral reactivation in 377 patients admitted to the Intensive Care Unit. The concomitant presence of herpesvirus and anellovirus (detected also from us) may have important clinical implications. Between potential coinfector, influencing the SARS‐CoV‐2's disease severity, find in HIV a valuable candidate. HIV and SARS‐CoV‐2 infections were found to be a dire combination and, indeed, the only HIV‐positive patient from our cohort had a fatal outcome (Figure 1D and Table 3). The purpose of this study is to characterize virome composition in COVID‐19 patients.

This study presents some limitations: it is only involving a single COVID‐19 patient cohort from the Campania region including only a single time‐point for each case. Detection of viruses employing supplementary specimen types such as oropharyngeal and broncho‐alveolar lavage fluids could also provide important additional information. Despite the intrinsic exploratory purpose of this study, it lifted up questions about whether some viruses with uncertain pathogenicity could be contributing to symptoms manifestation, complicating clinical outcome, or might be possible biomarkers of infection or host response.

CONFLICT OF INTERESTS

The authors declare that there are no conflict of interests.

AUTHOR CONTRIBUTIONS

Francesca Rizzo and Gianluigi Franci: Conceptualization. Roberta Astorri, Pasquale Pagliano, and Giorgio Giurato:Software analysis.Teresa Rocco, Ylenia D'Agostino, and Jessica Lamberti:Investigation.Carlo Ferravante, Giuseppina Sanna, Viola Melone, Aldo Manzin, Gianluigi Franci, and Giovanni Pecoraro: Data curation.Giuseppina Sanna and Carlo Ferravante: writing – original draft preparation. Gianluigi Franci, Francesca Rizzo, Giuseppina Sanna, Carlo Ferravante, and Giorgio Giurato: Writing – review and editing. Francesca Rizzo, Gianluigi Franci, Alessandro Weisz, Aldo Manzin, and Massimiliano Galdiero: Supervision.Giorgio Giurato, Massimiliano Galdiero, and Alessandro Weisz: Funding acquisition. All authors have read and agreed to the published version of the manuscript.

Supporting information.

Click here for additional data file.

Supporting information.

Click here for additional data file.