# Phylogeography and phylogeny of Rhinoviruses collected from Severe Acute Respiratory Infection (SARI) cases over successive epidemic periods in Tunisia

Sondes Haddad-Boubaker, Cherif Ben Hamda, Kais Ghedira, Khaoula Mefteh, Aida Bouafsoun, Ilhem Boutiba-Ben Boubaker, Amin Slim, Khaled Menif, Henda Triki, Mohamed Ali Ben Hadj Kacem, et al.

► To cite this version:

Sondes Haddad-Boubaker, Cherif Ben Hamda, Kais Ghedira, Khaoula Mefteh, Aida Bouafsoun, et al.. Phylogeography and phylogeny of Rhinoviruses collected from Severe Acute Respiratory Infection (SARI) cases over successive epidemic periods in Tunisia. PLoS ONE, 2021, 16 (11), pp.e0259859. 10.1371/journal.pone.0259859. pasteur-03546852

HAL Id: pasteur-03546852

https://riip.hal.science/pasteur-03546852

Submitted on 28 Jan 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
RESEARCH ARTICLE

# Phylogeography and phylogeny of Rhinoviruses collected from Severe Acute Respiratory Infection (SARI) cases over successive epidemic periods in Tunisia

Sondes Haddad-Boubaker<sup>1,nab,*</sup>, Cherif Ben Hamda<sup>2</sup>, Kais Ghedira<sup>2</sup>, Khaoula Mefteh<sup>1,3</sup>, Aida Bouafsoun<sup>1,3</sup>, Ilhem Boutiba-Ben Boubaker<sup>4,5</sup>, Amin Slim<sup>4</sup>, Khaled Menif<sup>6</sup>, Henda Triki<sup>nab</sup>, Mohamed Ali Ben Hadj Kacem<sup>6</sup>, Hanen Smaoui<sup>1,3</sup>

1 Laboratory of Microbiology, Bechir Hamza Children's Hospital, Bab-Saadoun Square, Tunis, Tunisia,

2 Laboratory of Bioinformatics, Biomathematics and Biostatistics, Institut Pasteur de Tunis, University of Tunis El-Manar, Tunis, Tunisia, 3 Microbiology of Children and Immunocompromised, Faculty of Medicine of Tunis, University of Tunis El-Manar, Tunis, Tunisia, 4 Laboratory of Microbiology, Charles Nicolle Hospital, Tunis, Tunisia, 5 Laboratory Research "Antimicrobial Resistance", Faculty of Medicine of Tunis University of Tunis El-Manar, Tunis, Tunisia, 6 Pediatric Intensive Care Unit, Bechir Hamza Children's Hospital in Tunis, Bab-Saadoun Square, Tunis Tunisia

nab Current address: Laboratory of Clinical Virology, WHO Regional Reference Laboratory for Poliomyelitis and Measles for the EMR, Institut Pasteur de Tunis, University of Tunis El-Manar, Tunis, Tunisia

nabb Current address: Laboratory of Viruses, Hosts and Vectors, Institut Pasteur de Tunis, University of Tunis El Manar, Tunis, Tunisia

* Sondeshaddadboubaker@gmail.com

<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

## OPEN ACCESS

**Citation:** Haddad-Boubaker S, Ben Hamda C, Ghedira K, Mefteh K, Bouafsoun A, Boutiba-Ben Boubaker I, et al. (2021) Phylogeography and phylogeny of Rhinoviruses collected from Severe Acute Respiratory Infection (SARI) cases over successive epidemic periods in Tunisia. PLoS ONE 16(11): e0259859. https://doi.org/10.1371/journal.pone.0259859

**Editor:** Ahmed S. Abdel-Moneim, Taif University, SAUDI ARABIA

**Received:** June 11, 2021

**Accepted:** October 27, 2021

**Published:** November 22, 2021

**Copyright:** © 2021 Haddad-Boubaker et al. This is an open access article distributed under the terms of the [Creative Commons Attribution License](https://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** All relevant data are within the paper and its [Supporting Information](#) files.

**Funding:** This work was supported by the Research and High Education Ministry of Tunisia [Research Laboratory LR18ES39]. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

## Abstract

Rhinoviruses (RV) are a major cause of Severe Acute Respiratory Infection (SARI) in children, with high genotypic diversity in different regions. However, RV type diversity remains unknown in several regions of the world. In this study, the genetic variability of the frequently circulating RV types in Northern Tunisia was investigated, using phylogenetic and phylogeographic analyses with a specific focus on the most frequent RV types: RV-A101 and RV-C45. This study concerned 13 RV types frequently circulating in Northern Tunisia. They were obtained from respiratory samples collected in 271 pediatric SARI cases, between September 2015 and November 2017. A total of 37 RV VP4-VP2 sequences, selected among a total of 49 generated sequences, was compared to 359 sequences from different regions of the world. Evolutionary analysis of RV-A101 and RV-C45 showed high genetic relationship between different Tunisian strains and Malaysian strains. RV-A101 and C45 progenitor viruses' dates were estimated in 1981 and 1995, respectively. Since the early 2000s, the two types had a wide spread throughout the world. Phylogenetic analyses of other frequently circulating strains showed significant homology of Tunisian strains from the same epidemic period, in contrast with earlier strains. The genetic relatedness of RV-A101 and RV-C45 might result from an introduction of viruses from different clades followed by local dissemination rather than a local persistence of an endemic clades along seasons. International traffic may play a key role in the spread of RV-A101, RV-C45, and other RVs.
**Competing interests:** The authors have declared that no competing interests exist.

**Introduction**

Rhinoviruses (RVs) are a common cause of upper respiratory diseases. In early childhood, they constitute a leading cause of severe lower respiratory infections causing important rates of hospitalization and mortality all around the world [1]. Clinical manifestations may vary from bronchiolitis and wheezing illnesses to pneumonia, and exacerbation of asthma [2, 3]. Besides, the infection may cause airway damage leading to subsequent asthma development [4].

RVs are non-enveloped, positive-sense RNA viruses, classified within the *Enterovirus* genus of the family *Picornaviridae* [5, 6]. They are characterized by a wide genetic diversity. Initially, they were classified into 101 distinct serotypes subdivided into two species, RV-A (76 serotypes) and RV-B (25 serotypes) [7]. Later, with the improvement of molecular-based diagnostic methods, new RVs were discovered [8] and RVs are now classified into three species, RV-A, RV-B, and RV-C, including 169 types [5]. In general, type identification is based on the VP4-VP2 coding region pairwise distances [9–11]. Nevertheless, using this approach, a significant number of RVs are still provisionally attributed. Thus, for accurate type identification, further genetic analysis using VP1 genomic coding region was also recommended [10, 11].

High genotypic diversity of RVs was reported in various regions of the world with co-circulation of different types during long periods or with turn-over of same types [12–14]. However, while RVs are highly variable viruses through mutations and recombination, leading to the emergence of divergent clades [15], studies on the genetic relationships between strains belonging to the same types remain scarce. It is not clear whether viruses from a RV clade can persist and establish continuous circulation and evolution in a geographic region or if RV circulation corresponds to the regular introduction of strains from different clades.

In Tunisia, we have recently described the circulation of 27 types of RV-A, B, and C, responsible for children's SARI cases, between September 2015 and December 2017 [16]. Some of the circulating viruses were sporadic, while others circulated for long periods. In the present study, lineage variability and genetic relationships among the frequently circulating RV types in Tunisia were investigated for a better comprehension of RV molecular epidemiology and circulation dynamics. A particular concern was assigned to the most frequent types RV-A101 and RV-C45, using phylogeographic analyses.

**Material and methods**

**Ethics statement**

The study was approved by the local Ethical Committee of Bechir Hamza Children Hospital of Tunis, Tunisia. It was performed under ethical standards according to the 1964 Declaration of Helsinki and its later amendments. The sequences investigated in this study were obtained from samples collected for diagnostic purposes, as part of routine laboratory tests for a panel of respiratory viruses. The Ethical Committee waived the need for the consent of parents or guardians of the minors included in the study, since samples were obtained for diagnostic purposes and only viral sequences were investigated. The specimens were used after de-identification of the samples with respect to patient anonymity and patient data protection.

**Study population**

Between September 2015 and November 2017, nasopharyngeal samples (n = 271) were collected from Tunisian children presenting Severe Acute Respiratory Infection (SARI) and admitted in the Pediatric Intensive Care Unit (PICU) of Bechir Hamza Children's Hospital in Tunis. This hospital receives patients from Tunis (the capital), around 11 districts from
northern Tunisia as well as complicated cases from other regions of the country. Among the 271 investigated samples, 57 were RV positive, using an in-house real-time PCR [16, 17]. Among the RV positive samples, 49 had sufficient residual samples to be used for VP4-VP2 PCR amplification and sequencing [7] and the 49 VP4-VP2 sequences were generated [16].

## RV VP4-VP2 sequences

The present study focused on the most frequently detected RV types in Tunisia, in the aim to evaluate the contribution of endemic circulation and virus importation from a region to another in the circulation dynamics of RVs. Thus, among the 49 generated sequences, 37 sequences, RV-A (n = 23); RV-B (n = 3) and RV-C (n = 11) were included. As we have previously demonstrated, they correspond to the most frequent RV types A and C in Northern Tunisian between September 2015 and November 2017: RV-A (8 types) and RV-C (3 types), and all RV-B circulating in the same period (3 types) (Table 1) [16]. Viral RNA was extracted using the QIAamp MinElute Virus Spin Kit (Qiagen, GmbH, Hilden, Germany) according to the manufacturer's recommendations. The VP4-VP2 genomic region was amplified using the primer pair 9565-reverse and 9895-forward, producing a fragment of approximately 549 nucleotides [7]. The reaction was performed using the One-step RT-PCR kit from QIAGEN (GmbH, Hilden, Germany) in a total PCR mixture of 25µl containing 2.5µl of extracted RNA, 5µl of RT-PCR buffer, 400µM of each dNTP, 1µM of each primer, 1µl of enzyme Mix and RNase-free water. The thermal cycling conditions were as follows: 20mn at 50°C for reverse transcription, 15mn at 95°C for initial denaturation, 35 cycles of 1mn at 95°C for denaturation, 1mn at 60°C for annealing, 1mn at 72°C for the extension, and 7mn at 72°C for a final extension step.

The VP4-VP2 amplicons were initially visualized in ethidium-bromide containing 1% agarose gel. The amplicons of the expected size were then purified by the Qiaquick PCR purification kit from Qiagen (Gmbh, Hilden, Germany). Sequences were obtained by automated sequencing using the Big Dye terminator chemistry, according to the manufacturer's protocol (Applied Biosystems) in an automated sequencer (ABI 3130). Accession numbers of investigated sequences are indicated in Table 1.

## Data collection

In this study, a total of 359 VP4-VP2 RV sequences from several countries of the world, collected from the NCBI GenBank database (http://www.ncbi.nlm.nih.gov/nucleotide were investigated, in addition to the sequences obtained as part of the present study). Only well-annotated sequences, including the time of isolation and the geographic origin were selected. Duplicate sequences from same country same year were removed and a unique sequence of duplicates was retained. Phylogeographic analyses used 97 and 30 sequences of RV-A101 and RV-C45, respectively (S1 and S2 Tables). Phylogenetic analyses used 232 sequences of different RV types: RV-A12 (n = 27), RV-A78 (n = 24), RV-A32 (n = 21), RV-A89 (n = 14), RV-A15 (n = 21), RV-A16 (n = 11), RV-A19 (n = 24), RV-B72 (n = 22), RV-A48 (n = 13), RV-A42 (n = 15), RV-C53 (n = 18) and RV-C43 (n = 22). The accession numbers of the used sequences are indicated in the drawn trees.

## Bayesian phylogeography reconstruction

Bayesian analyses concerned the most circulating RV types in Tunisia, RV-A101 and RV-C45. Sequences were first aligned using Muscle 3.6 software [18]. Time-scaled phylogenies were inferred by Bayesian Markov Chain Monte Carlo (MCMC) sampling using BEAST v1.10.4 [19]. In a first step, the best substitution models were identified using the likelihood-based
<table><caption>Table 1. Description of Tunisian sequences investigated.</caption><thead><tr><td>Sequence name</td><td>Date of isolation</td><td>Accession number</td><td>Type</td></tr></thead><tbody><tr><td>HRV38.Tun.02.2016</td><td>Feb/2016</td><td>MN583163</td><td>RV-A101</td></tr><tr><td>HRV37.Tun.03.2016</td><td>Mar/2016</td><td>MN583162</td><td></td></tr><tr><td>HRV04.Tun.03.2016</td><td>Mar/2016</td><td>MN583129</td><td></td></tr><tr><td>HRV06.Tun.08.2016</td><td>Aug/2016</td><td>MN583131</td><td></td></tr><tr><td>HRV25.Tun.09.2016</td><td>Sep/2016</td><td>MN583150</td><td></td></tr><tr><td>HRV26.Tun.10.2016</td><td>Nov/2016</td><td>MN583151</td><td></td></tr><tr><td>HRV29.Tun.01.2017</td><td>Jan/2017</td><td>MN583154</td><td></td></tr><tr><td>HRV30.Tun.12.2016</td><td>Dec/2016</td><td>MN583155</td><td>RV-A32</td></tr><tr><td>HRV10.Tun.01.2017</td><td>Jan/2017</td><td>MN583135</td><td></td></tr><tr><td>HRV14.Tun.04.2017</td><td>Apr/2017</td><td>MN583139</td><td></td></tr><tr><td>HRV35.Tun.09.2015</td><td>Sep/2015</td><td>MN583160</td><td>RV-A12</td></tr><tr><td>HRV36.Tun.01.2017</td><td>Jan/2017</td><td>MN583161</td><td></td></tr><tr><td>HRV11.Tun.02.2017</td><td>Feb/2017</td><td>MN583136</td><td></td></tr><tr><td>HRV01.Tun.09.2015</td><td>Sep/2015</td><td>MN583126</td><td>RV-A15</td></tr><tr><td>HRV02.Tun.10.2015</td><td>Oct/2015</td><td>MN583127</td><td></td></tr><tr><td>HRV03.Tun.11.2015</td><td>Nov/2015</td><td>MN583128</td><td>RV-A78</td></tr><tr><td>HRV45.Tun.10.2017</td><td>Oct/2017</td><td>MN583170</td><td></td></tr><tr><td>HRV17.Tun.07.2017</td><td>Jul/2017</td><td>MN583142</td><td>RV-A16</td></tr><tr><td>HRV34.Tun.11.2017</td><td>Nov/2017</td><td>MN583159</td><td></td></tr><tr><td>HRV20.Tun.02.2016</td><td>Feb/2016</td><td>MN583145</td><td>RV-A19</td></tr><tr><td>HRV24.Tun.08.2016</td><td>Aug/2016</td><td>MN583149</td><td></td></tr><tr><td>HRV28.Tun.10.2015</td><td>Oct/2015</td><td>MN583153</td><td>RV-A89</td></tr><tr><td>HRV22.Tun.04.2016</td><td>Apr/2016</td><td>MN583147</td><td></td></tr><tr><td>HRV13.Tun.03.2017</td><td>Mar/2017</td><td>MN583138</td><td>RV-C45</td></tr><tr><td>HRV42.Tun.03.2017</td><td>Mar/2017</td><td>MN583167</td><td></td></tr><tr><td>HRV48.Tun.03.2017</td><td>Mar/2017</td><td>MN583173</td><td></td></tr><tr><td>HRV16.Tun.05.2017</td><td>May/2017</td><td>MN583141</td><td></td></tr><tr><td>HRV44.Tun.06.2017</td><td>Jun/2017</td><td>MN583169</td><td></td></tr><tr><td>HRV39.Tun.06.2017</td><td>Jun/2017</td><td>MN583164</td><td></td></tr><tr><td>HRV08.Tun.10.2016</td><td>Oct/2016</td><td>MN583133</td><td>RV-C43</td></tr><tr><td>HRV27.Tun.10.2016</td><td>Oct/2016</td><td>MN583152</td><td></td></tr><tr><td>HRV46.Tun.10.2016</td><td>Oct/2016</td><td>MN583171</td><td></td></tr><tr><td>HRV07.Tun.09.2016</td><td>Sep/2016</td><td>MN583132</td><td>RV-C53</td></tr><tr><td>HRV43.Tun.05.2017</td><td>May/2017</td><td>MN583168</td><td></td></tr><tr><td>HRV32.Tun.08.2016</td><td>Jul/2016</td><td>MN583157</td><td>RV-B72</td></tr><tr><td>HRV31.Tun.12.2016</td><td>Sep/2016</td><td>MN583156</td><td>RV-B42</td></tr><tr><td>HRV49.Tun.03.2017</td><td>Mar/2017</td><td>MN583174</td><td>RV-B48</td></tr></tbody></table>

https://doi.org/10.1371/journal.pone.0259859.t001

criteria (AIC) (S1 Fig), and the Smart Selection Tool (SMS) [20]. The time to the most recent common ancestor (tMRCA) was evaluated using the TN93 model with gamma invariant sites for RV-A101 and using the GTR model with gamma invariant sites for RV-C45 based on the SMS tool prediction. Five parametric models of coalescent population growth were then used: Bayesian skyline, constant population size, exponential growth, expansion growth and logistic growth under the strict, uncorrelated relaxed and random local clocks. The MCMC chains were run separately for each data set (RV-A101 and RV-C45). A total of 2 × 10⁷ MCMC iterations were run, with a burn-in of 2 × 10⁶ iterations and a sampling frequency of 6,000. In total, fifteen independent runs (five models under three clocks) were run. Convergence and effective
<table><caption>Table 2. Fitting models bayes factors and HPD95% related to RV-A101.</caption><thead><tr><th rowspan="2">Clock/Tree prior</th><th>Constant</th><th>Exponential</th><th>Expansion</th><th>Bayesian</th><th>Logistic</th></tr><tr><th>Size</th><th>Growth</th><th>Growth</th><th>Skyline</th><th>Growth</th></tr></thead><tbody><tr><td rowspan="2">Random local</td><td>-2 048,50</td><td>-2 049,50</td><td>-2 048,52</td><td>-2 038,93</td><td>-1 513,96</td></tr><tr><td>1980(1977,1982)</td><td>1980(1977,1982)</td><td>1980(1977,1982)</td><td>1980(1976,1982)</td><td>1981(1977,1982)</td></tr><tr><td rowspan="2">Strict</td><td>-2 048,11</td><td>-2 049,32</td><td>-2 049,10</td><td>-2 041,37</td><td>-1 555,08</td></tr><tr><td>1980(1976,1982)</td><td>1980(1977,1982)</td><td>1980(1977,1982)</td><td>1980(1976,1982)</td><td>1980(1977,1982)</td></tr><tr><td rowspan="2">Uncorrelated Relaxed</td><td>-2 051,79</td><td>-2 055,36</td><td>-2 055,15</td><td>-2 042,35</td><td>-1 650,94</td></tr><tr><td>1980(1977,1982)</td><td>1980(1977,1982)</td><td>1980(1977,1982)</td><td>1980(1977,1982)</td><td>1980(1977,1982)</td></tr></tbody></table>

https://doi.org/10.1371/journal.pone.0259859.t002

sampling size (ESS) of estimates were assessed by visual inspection using Tracer v1.7 [21]. The best fitting models were then selected by the calculation of a Bayes Factor (BF) [22], using Marginal Likelihoods Estimation (MLE) mean of path sampling (PS) [23] and stepping-stone sampling (SS) (Tables 2 and 3) [24]. Finally, the Maximum Clade Credibility (MCC) trees with temporal and spatial annotation were constructed with a 10% burn-in removed using Tree Annotator being part of the BEAST package. Phylogenetic trees were generated using FigTree v1.4.2 (http://beast.bio.ed.ac.uk/figtree). Processed data sets were finally visualized using the “SPREAD 3” program on Google Earth Pro version 7.3 (http://www.google.com/earth/download/ge/) to produce a graphical animation of the continuous movement dynamics of RV-A101 and RV-C45 investigated sequences in the key markup language (kml) file format. These kml files contain information on the routes and times of virus movements, that can be displayed using Google Earth.

## Phylogenetic analyses

In the aim to study genetic diversity and to identify the epidemiological origin of the other circulating RV types in Northern Tunisia, phylogenetic analyses were conducted using Mega version 6.0. (https://www.megasoftware.net/). The sequences were first aligned using Muscle 3.6 software. Dendrograms were drawn using the Neighbor-joining method and confirmed with 1000 bootstrap replicates.

## Results

### Epidemic history of RV-A101 and RV-C45

**Evolutionary analysis.** Fig 1 displays the maximum clade credibility trees (MCC) of RV-A101 (Fig 1A) and RV-C45 (Fig 1B) for 100 and 36 sequences, respectively, from several countries. The tree in (Fig 1A) shows that at least four RV-A101 lineages are circulating in the world. However, no consistent geographical or time-related cluster could be identified.

<table><caption>Table 3. Fitting models bayes factors and HPD95% related to RV-C45.</caption><thead><tr><th rowspan="2">Clock/Tree prior</th><th>Constant</th><th>Exponential</th><th>Expansion</th><th>Bayesian</th><th>Logistic</th></tr><tr><th>Size</th><th>Growth</th><th>Growth</th><th>Skyline</th><th>Growth</th></tr></thead><tbody><tr><td rowspan="2">Random local</td><td>-1 785,53</td><td>-1 788,06</td><td>-1 785,80</td><td>-1 784,33</td><td>-1 788,40</td></tr><tr><td>1995(1989,1996)</td><td>1979(1961,1991)</td><td>1979(1961,1991)</td><td>1982(1942, 1996)</td><td>1982(1963,1999)</td></tr><tr><td rowspan="2">Strict</td><td>-1 792,07</td><td>-1 794,26</td><td>-1 794,06</td><td>-1 795,69</td><td>-1 791,73</td></tr><tr><td>1943(1976,1977)</td><td>1947(1913,1975)</td><td>1944(1895,1971)</td><td>1950(1767, 1977)</td><td>1949(1928,1977)</td></tr><tr><td rowspan="2">Uncorrelated Relaxed</td><td>-1 781,67</td><td>-1 783,83</td><td>-1 781,15</td><td>-1 778,40</td><td>-1 783,60</td></tr><tr><td>1993(1957,1996)</td><td>1993(1996,1987)</td><td>1993(1996,1987)</td><td>1995(1976, 1996)</td><td>1995(1992,1996)</td></tr></tbody></table>

https://doi.org/10.1371/journal.pone.0259859.t003
<figure><img src="image_4.png" /><figcaption>**Fig 1.** Maximum clade credibility (MCC) trees for RV-A101 (A) and RV-C45 (B). The branches are colored according to the most probable ancestor location of their descendant nodes. Tunisian RV-A101 and RV-C45 strains are highlighted in a Cyan and Purple rectangular, respectively.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0259859.g001

Tunisian sequences are highlighted with a cyan background color. The tree shows that almost all lineages are circulating in the USA. The Tunisian strains seem profoundly related to Malaysian ones. Likewise, several clades for RV-C45, could also be distinguished (Fig 1B) forming four potential monophyletic groups. The Tunisian sequences are highlighted with a purple background color. Tunisian strains seem extremely related to Mongolian, Malaysian, and Indian strains.

**Bayesian skyline plots.** We performed Bayesian coalescent analysis on the RV-A101 and RV-C45 sequences from Tunisia and other countries using five parametric models of coalescent population growth under the strict, uncorrelated relaxed, and random local clocks. For each demographic and molecular clock model, chain lengths of 20 million were used and sampled every 6,000 states. The data were analyzed using the Bayesian Skyline Plot (BSP) that represents the estimated change in the effective number of infected individuals over time for the RV-A101 and RV-C45 viruses. Fig 2 displays both Bayesian skyline plots corresponding to RV-A101 (Fig 2A) and RV-C45 (Fig 2B), respectively.

Fig 2A shows that the RV-A101 epidemic history was characterized by four phases of epidemic population growth including an initial period of relatively constant population size from 1980 to 2001, followed by a brief decline between 2001 and 2003, a third period of exponential growth until 2005 and a final phase with a significant drop from 2005 to 2019. Fig 2B

<figure><img src="image_5.png" /><figcaption>**Fig 2.** Bayesian skyline plot of RV-A101 (A) and RV-C45 (B). Bayesian skyline plots of the RV-A101 (A) and RV-C45 (B) genes show the changes in effective population size (genetic diversity) through time. The thick solid line indicates the median value, and the blue area is the 95% Highest Posterior Density (HPD) of the genetic diversity estimates.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0259859.g002
shows that the RV-C45 epidemic history was characterized by two phases including an epi-
demic population growth period between 1980 and 2008 followed by a second period of a
decline from 2008 to 2019. The waves on population size could be a consequence of an increase
of international trades and exchanges during the last decades.

**Most Recent Common Ancestor (tMRCA).** The progenitor virus of the RV-A101 that ultimately gave rise to all RV-A101 sequences was estimated around 1981 with an HPD (Highest Posterior Density) 95% [1977, 1982] and the one that eventually gave rise to all RV-C45 circulating in the world was estimated around 1995 with an HPD95% [1976, 1996], respectively.

**Phylogeographic analysis.** Phylogeography generally describes the geographical distribution of lineages and has been used to reconstruct the geospatial dynamics of disease spread and the diffusion process of RV-A101 and RV-C45 in Tunisia and in the world by mapping the spatial estimates annotated in the MCC trees on Google Earth Pro version 7.3 (http://www.google.com/earth/download/ge/). This mapping allows visualizing the virus's geographic spread process over time. The links between different geographic regions represent branches in the MCC tree on which virus migration occurs and circle areas reflect the number of branches maintaining a particular location at that time point. The panels in Fig 3 show the temporal dynamics of RV-A101(Fig 3A) and RV-C45 (Fig 3B) spatial dispersal processes in the world.

Our results suggest that the most ancient RV-A101 virus dated around 1981 in the world. Based on the investigated RV-A101 sequences, the phylodynamic and phylogeographic analysis revealed that the virus started its migration from the USA around 1981 to Australia in 1988, Asia (the Philippines and South Korea) starting from 2003, European countries including United Kingdom, Germany, Spain, and Italy from 2005 to 2009 than East Africa around 2010. The USA would have been the origin of the spread to other countries worldwide. The results showed that RV-A101 isolates were introduced to the USA several times.

The virus spread in South Korea and reached India, Singapore, and China in 2011. It migrated into Colombia through Spain around 2013. Concerning North Africa, a migration of the virus from China to Egypt would have occurred in 2011 and then to Taiwan in 2012. From the USA, the virus would have reached Malaysia in 2011 where it spreads from 2011 to 2015, then migrated to Tunisia around 2016–2017.

For the RV-C45, the virus seems to have migrated from Finland around 1998 to the USA in 2001 where it started to spread then continued its migration to reach Japan in 2004 and from there to Thailand, around 2006, where it spread until 2007. During this period, the virus migrated from Japan to Argentina and Egypt around 2012 and from Thailand to China around 2010. In 2011, the virus migrated again from Thailand to Japan and reached Cyprus in 2014. Meanwhile, in 2013 RV-C45 migrated from the USA to Mongolia, Japan, and Thailand. RV-C45 then migrates from Mongolia to reach Tunisia in 2016 and Malaysia around 2015. It spreads in Tunisia until 2017 then migrates into India in 2018.

**Phylogenetic analysis of other circulating RV-A**

In the aim to analyze genetic variability and identify the epidemiological origin of RV strains circulating in Tunisia, during the study period, we assessed phylogenetic analyses of strains circulating during at least 5 months (RV-A12, RV-A78, RV-A32, RV-A89, RV-A15, RV-A16, RV-A19, RV-A18) (Fig 4A).

Three strains of RV-A12 were detected in 2015 and 2017. The two strains circulating in 2017 were identical and slightly different from the one that circulated in 2015 (96.1% identity). With other Kenyan strains, they seem to constitute an independent clade confirmed by a high bootstrap value (82%). The most related strain was a Kenyan strain (MH459725) detected in 2016 (98.1–98.8% identity) (Fig 4A.1).
<figure><img src="image_6.png" /><figcaption>**Fig 3.** Temporal dynamics of RV-A101 (Panel A) and RV-C45 (Panel B) and spatial dispersal in the world. Google Earth (https://earth.google.com). The figure is similar but not identical to the original image and is therefore for illustrative purposes only.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0259859.g003

Two RV-A78 Tunisian strains were detected in 2015 and 2017, the other was genetically quite different (90.2% identity) and belong to two different clades confirmed by high bootstrap values (77% and 91%): the one detected in 2015 was highly related to strains that originated from the Americas and especially to a strain isolated in 2014 (MF160344) with 96.2% identity rate. The strain detected in 2017 clustered with African strains detected in Egypt (97.3% identity rate) and Kenya (97.1% identity rate) (Fig 4A.2).

Regarding RV-A32, three Tunisian strains were detected in 2016 and 2017. The two strains from 2017 shared 99.6% nucleotide identity; and had 99.3%, 98.9% identity rates with the strain from 2016. These strains were deeply related to a strain from Hong Kong (HQ444790) with identity rates ranging from 98% to 99.1%. They constitute an independent clade supported by a high bootstrap value (96%). Another Tunisian strain (RV-A32) obtained in 2010, was previously published in GenBank. However, it is slightly different from the 2016 and 2017 Tunisian sequences (Identity rates ranging from 97.2% to 98%) (Fig 4A.3).
<figure><img src="image_7.png" /></figure>

A.1

<figure><img src="image_8.png" /></figure>

A.2

<figure><img src="image_9.png" /></figure>

A.3

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

A.7

<figure><img src="image_14.png" /></figure>

A.6

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

B.3

<figure><img src="image_17.png" /></figure>

C.1

<figure><img src="image_18.png" /></figure>

C.2
**Fig 4.** Molecular typing of Tunisian RV strains obtained during the study period (2015–2017). **A.** Molecular typing of RV-A strains: A1. RV-A12, A2. RV-A78, A3. RV-A32, A4. RV-A89, A5. RV-A15, A6. RV-A16, A7. RV-A19; **B.** Molecular typing of RV-B strains: B1. RV-B72, B2. RV-B48, B3. RV-B42; **C.** Molecular typing of RV-C strains: C1. RV-C53, C2. RV-C43. The evolutionary history was inferred using the Neighbor joining method. The percentage of replicate trees in which the associated taxa clustered together in the bootstrap test (1000 replicates) is shown next to the branches. Evolutionary analyses were conducted in MEGA6. Investigated Tunisian sequences are indicated in bold and Red; The most related sequences to investigated ones are indicated in bold and Green; Previously available Tunisian are indicated in bold in Blue.

https://doi.org/10.1371/journal.pone.0259859.g004

For RV-A89, two Tunisian strains were detected in 2015 and 2016; they are identical and closely related to Pacific Asiatic strains that circulated between 2010 and 2014 (98–98.4% identity rate) (Fig 4A.4).

Two strains of RV-A15 were detected, in Tunisia, in 2015 and 2016. Both are closely related to each other (99% identity rate) and with strains from the USA (KU237003) and Spain (KF957929) with 98.3% and 99.3% identity rates, respectively (Fig 4A.5).

The two RV-A16 Tunisian strains here investigated were slightly different from each other (96% identity rate). and clustered with pacific Asian strains. The HRV34.TUN.11.2017 strain was closely related to a strain obtained in Thailand in 2014 (98.6% identity) while the HRV17. TUN.07.2017 strain was slightly different and presented approximately the same identity rates with all strains included in this clade (95.8–96% identity rates) (Fig 4A.6).

For RV-A19, the phylogenetic tree grouped all strains into two clusters supported by high bootstrap values (73 and 85%). The newly investigated Tunisian strains are different from each other. Indeed, the strain HRV24.TUN.08.2016 seems to be genetically independent but at the same time, it shares an identity rate of 98% with a Peruvian strain (JX129441) taken in 2010. However, the strain HRV20.TUN.02.2016 is related to a Kenyan strain (MH459624) with an identity rate of 99.2%. Other Tunisian strains, obtained in 2011, were available in GenBank. They seem to belong to different clades. The first strain (KM234342) is identical to an Egyptian (KT286895), Spanish (JN855961), and Netherlander strains (KP736959). The second Tunisian strain (KM234342) was related to another Netherlander (JQ042514) and Mongolian (462714) strains (99.1 and 98.9% identity rates respectively) (Fig 4A.7).

## Phylogenetic analysis of RV-B

All three detected RV-B viruses were investigated (Fig 4B). Phylogenetic analysis of RV-B72 demonstrates that the Tunisian strain is strongly related to strains from Singapore (MH648104) obtained in 2009 and Malaysia (KY09355) in 2012, with an identity rate of 99.1 and 98.9%, respectively (Fig 4B.1). Regarding RV-B48, the Tunisian strain is related to Kenyan (KX348029) and Singaporean (MH648106) strains, obtained in 2010 with an identity rate of 97.4% (Fig 4B.2). Also, for RV-B42, the Tunisian strain was deeply related to a Chinese (GU568096) strain obtained in 2007 with a 97.3% identity rate (Fig 4B.3).

## Phylogenetic analysis of other circulating RV-C

Regarding RV-C, we are interested in the most circulating types, either RV-C45: RV-C53; RV-C43 (Fig 4C).

Two Tunisian strains of RV-C53 were circulating in Tunisia in 2016 and 2017. Both strains were closely related to each other (99.3% identity rate) and with a strain from Kenya (MH459850) (97.8 and 98.6% identity rates respectively) (Fig 4C.1).

Regarding RV-C43, three Tunisian strains were detected during the study period. Among them, HRV08.Tun.10.2016 and HRV46.Tun.10.2016 were closely linked to each other (99.1% identity rate) and are related to the Kenyan strains (MH460082, MH460081) with 97.7 and
98.6% identity rates respectively. In the other hand, the HRV27.Tun.10.2016 strain was closely related to a strain obtained in Thailand (KR054525) with 98% identity rate (Fig 4C.2).

Discussion

RVs represent one of the leading causes of hospitalization and mortality in neonates and young children. Thus, their molecular characterization is of major interest to understand their epidemiology and genetic evolution. Previous studies focused mainly on the prevalence and type diversity of circulating RVs, with few data on the phylogeny and phylogeography of the different types [12–14]. In-depth molecular studies of circulating strains can help understanding the genetic diversity within each type, identifying endemic and persistent clades in a region and evaluating the contribution of virus importation and exchange from a region to another in the virus spread. To the best of our knowledge, this is the first study describing the phylogeny and phylogeography of the RV types circulating in a region.

In this study, phylogeographic analyses were conducted on the most frequently detected RV types in Northern Tunisia (RV-A101 and RV-C45), from September 2015 to November 2017 [16]. Both types belong to the most frequently circulating RV species in the world (A and C); they were detected in Europe, the USA, Asia, and Africa [11]. Tunisian strains were closely related to each other and to other strains from Asia suggesting importation events from Asiatic countries followed by the dissemination of different strains from the same clade during 11 and 4 months, respectively. After January and June 2017, and during the study period, the circulation of both types seems interrupted and, likely, replaced by other RV types. The global topology of the phylogenetic trees showed at least four potential monophyletic groups of RV-A101 and RV-C45 that circulated in the world. However, the evolution and spread of virus clusters do not respect any geographic distribution; they seem disseminated throughout the world without any geographic restriction. This finding may indicate a regular exchange of strains from one region to another and supports the hypothesis that circulation of RV-A101 and RV-C45 occurs through regular introductions of viruses from different clades rather than a local persistence of the same clade from a season to another. Unfortunately, a very restricted number of RV sequences from several seasonal epidemics in Tunisia are available. Nevertheless, considering the available RV sequences, our phylogenic analysis of other RV-A and RV-C types displayed a remarkable homogeneity among strains obtained during the same seasonal period; while comparing to strains from 2011, they seem different and related to strains from distinct geographic origins. The same result was also obtained for RV-A78 and RV-A12 sequences where strains obtained in 2015 were related to strains from different geographic origins compared to strains obtained in the following epidemic period in 2017. This suggests that, in general, the genetic diversity of RV types results from regular introductions of new viruses rather than a continuous local evolution of specific genotypes. Long-term surveillance of RV circulating types during multiple seasonal periods may further support this finding. Also, investigation of another genetic region, such as the VP1 region or full genomic RV sequences may support our hypothesis. In our study, given the limited number of samples obtained from young children and the multiple bacterial and virologic investigation conducted for diagnosis, the residual amounts of samples were only sufficient for sequencing the VP2-VP4 region. Other studies concerning other respiratory viruses such as the influenzae virus revealed that the genetic diversity is also a result of frequent introductions of genetically divergent viruses from different clades during seasonal epidemics [25–27]. For instance, genetic analysis of the H3N2 virus over 15 years in China did not demonstrate any continuous virus clade persistence over seasonal epidemics [28].

On the other hand, phylogeographic analyses highlights the significant spread of both virus types from one country to another. The present work suggests that RV-A101 and C45 progenitor
viruses rise approximately in the eighteenth and nineteenth, respectively. Since the early 2000s, they knew a big spread all over the world and were detected in The United States of America, Australia, Asian countries, European countries, and some African countries, despite the restricted number of available sequences from some regions. The level of international viral traffic was high and rapid which suggests travel-related dissemination of RV-A101 and C45. Given the mild clinical presentation of RV infection in adults, travelers may significantly contribute to the spread of RVs, globally. Furthermore, phylogenetic analysis of other detected RVs circulating in Tunisia revealed strong epidemiological links with strains from the USA, Asian and African countries. During our study period, different virus types from different geographic origins co-circulated simultaneously, with extinction and replacement of some strains by other strains from the same type or other types but in general from different geographic origins.

## Conclusions

In this study, the genetic diversity of the most frequently circulating RVs was investigated, during two successive SARI epidemic periods in Northern Tunisia, using phylogeny and phytogeography. It highlights the high genetic homology among viruses from the same type during a seasonal epidemic period with lower homology rates with strains from another period. The study suggests high international spread of RV-A101 and RV-C45. Thus, circulation of RV-A101 and RV-C45 seems to result from regular introductions of viruses from different clades followed by their local dissemination, rather than a local persistence of viruses from the same clade along seasons. International traffic may play a pivotal role in the spread of RV-A101, RV-C45, and other RV types. Long-term surveillance of RV circulating types during multiple seasonal periods may further support these findings.

## Supporting information

S1 Table. Description of RV-A101 sequences included in this study.
(DOCX)

S2 Table. Description of RV-C45 sequences included in this study.
(DOCX)

S1 Fig. Smart selection tool output of RV-A101 and RV-C45.
(JPG)

## Acknowledgments

The authors are grateful to Pr Amel Kechrid, Pr Nejla Ben Jeballah and the staff of the Pediatric Intensive Care Unit (PICU) of Bechir Hamza Children's Hospital in Tunis, Tunisia, for their help in collecting samples and patient's information.

## Author Contributions

**Conceptualization:** Sondes Haddad-Boubaker, Kais Ghedira, Mohamed Ali Ben Hadj Kacem, Hanen Smaoui.

**Data curation:** Sondes Haddad-Boubaker, Cherif Ben Hamda, Kais Ghedira, Khaoula Mefteh, Khaled Menif.

**Formal analysis:** Sondes Haddad-Boubaker.

**Investigation:** Sondes Haddad-Boubaker.
**Methodology:** Sondes Haddad-Boubaker.

**Resources:** Khaoula Mefteh, Aida Bouafsoun, Khaled Menif.

**Supervision:** Sondes Haddad-Boubaker, Hanen Smaoui.

**Validation:** Sondes Haddad-Boubaker.

**Visualization:** Sondes Haddad-Boubaker.

**Writing – original draft:** Sondes Haddad-Boubaker, Cherif Ben Hamda.

**Writing – review & editing:** Sondes Haddad-Boubaker, Cherif Ben Hamda, Kais Ghedira, Ilhem Boutiba-Ben Boubaker, Amin Slim, Khaled Menif, Henda Triki, Mohamed Ali Ben Hadj Kacem, Hanen Smaoui.

## References

1. Vandini S, Biagi C, Fischer M, Lanari M. Impact of Rhinovirus Infections in Children. Viruses. 2019; 11 (6): 521. https://doi.org/10.3390/v11060521 PMID: 31195744

2. Cheuk DK, Tang IW, Chan KH, Woo PC, Peiris MJ, Chiu SS. Rhinovirus infection in hospitalized children in Hong Kong: a prospective study. Pediatr Infect Dis J. 2007; 26: 995–1000. https://doi.org/10.1097/INF.0b013e3181586b63 PMID: 17984805

3. Miller EK, Lu X, Erdman DD, Poehling KA, Zhu Y, Griffin MR, et al. Rhinovirus associated hospitalizations in young children. J Infect Dis. 2007; 195: 773–781. https://doi.org/10.1086/511821 PMID: 17299706

4. Jackson DJ, Gangnon RE, Evans MD, Roberg KA, Anderson EL, Pappas TE, et al. Wheezing rhinovirus illnesses in early lifepredict asthma development in high-risk children. Am J Respir Crit Care Med. 2008; 178: 667–72. https://doi.org/10.1164/rccm.200802-309OC PMID: 18565953

5. International Committee on Taxonomy of Viruses (ICTV), The online 10th report of the International Committee on Taxonomy of Viruses. 2019. Available: https://ictv.global/taxonomy. Accessed: 03/02/2021

6. Simmonds P, Gorbalenya AE, Harvala H, Hovi T, Knowles NJ, Lindberg A M, et al. Recommendations for the nomenclature of enteroviruses and rhinoviruses. Archives of Virology. 2020; 165: 793–797. https://doi.org/10.1007/s00705-019-04520-6 PMID: 31980941

7. Savolainen C, Blomqvist S, Mulders MN, Savolainen C, Blomqvist S, Mulders MN, et al. Genetic clustering of all 102 rhinovirus prototype strains: serotype 87 is close to enterovirus 70. J Gen Virol. 2002; 83: 333–340. https://doi.org/10.1099/0022-1317-83-2-333 PMID: 11807226

8. Lamson D, Renwick N, Kapoor V, Liu Z, Palacios G, Ju J, et al. MassTag polymerase-chain-reaction detection of respiratory pathogens, including a new rhinovirus genotype, that caused influenza-like illness in New York State during 2004–2005. J Infect Dis. 2006; 194: 1398–1402. https://doi.org/10.1086/508551 PMID: 17054069

9. McIntyre CL, McWilliam Leitch EC, Savolainen-Kopra C, Hovi T, Simmonds P. Analysis of genetic diversity and sites of recombination in rhinovirus species C. J Virol. 2010; 84: 10297–310. https://doi.org/10.1128/JVI.00962-10 PMID: 20668080

10. Simmonds P, McIntyre C, Savolainen-Kopra C, Tapparel C, Mackay IM, Hovi T. Proposals for the classification of rhinovirus species C into genotypically assigned types. J Gen Virol. 2010; 91: 2409–2419. https://doi.org/10.1099/vir.0.023994-0 PMID: 20610666

11. McIntyre CL, Knowles NJ, Simmonds P. Proposals for the classification of rhinovirus species A, B and C into genotypically assigned types. J Gen Virol. 2013; 94: 1791–1806. https://doi.org/10.1099/vir.0.053686-0 PMID: 23677786

12. Zhao Y, Shen J, Wu B, Liu G, Lu R, Tan W. Genotypic diversity and epidemiology of Rhinovirus among children with Severe Acute Respiratory Tract Infection in Shanghai, 2013–2015. Front Microbiol. 2018; 9: 1836. https://doi.org/10.3389/fmicb.2018.01836 PMID: 30131797

13. Hung HM, Yang SL, Chen CJ, Chiu CH, Kuo CY, Huang KA, et al. Molecular epidemiology and clinical features of rhinovirus infections among hospitalized patients in a medical center in Taiwan. J Microbiol Immunol Infect. 2019; 52: 233–241. https://doi.org/10.1016/j.jmii.2018.08.009 PMID: 30201131

14. Lu QB, Wo Y, Wang LY, Wang HY, Huang DD, Zhang XA, et al. Molecular epidemiology of rhinovirus in children with acute respiratory diseases in Chongqing, China. Sci Rep. 2014; 4: 6686. https://doi.org/10.1038/srep06686 PMID: 25328042
15. Rathe JA, Liu X, Tallon LJ, Gern JE, Liggett SB. Full-genome sequence and analysis of a novel rhinovirus strain within a divergent RV-A clade. Arch Virol. 2010; 155: 83-7. https://doi.org/10.1007/s00705-009-0549-8 PMID: 19936613

16. Haddad-Boubaker S, Mefteh K, Mejri C, Bouaffsoun A, El Moussi A, Boutiba I, et al. High genotypic diversity of Rhinoviruses obtained from Tunisian children with severe acute respiratory infection. J Infect Dev Ctriest. 2021; 15: 726-735. https://doi.org/10.3855/jidc.12880 PMID: 34106898

17. Volle R, Nourrisson C, Mirand A, Regagnon C, Chambon M, Henquell C, et al. Quantitative real-time RT-PCR assay for research studies on enterovirus infections in the central nervous system. J Virol Methods. 2012; 185: 142-8. https://doi.org/10.1016/j.jviromet.2012.06.019 PMID: 22766179

18. Edgar RC. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Research. 2004; 32: 1792-97. https://doi.org/10.1093/nar/gkh340 PMID: 15034147

19. Drummond AJ, Rambaut A. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol Biol. 2007; 7: 214. https://doi.org/10.1186/1471-2148-7-214 PMID: 17996036

20. Lefort V, Longueville JE, Gascuel O. SMS: Smart Model Selection in PhyML. Molecular biology and evolution. 2017; 34: 2422-2424. https://doi.org/10.1093/molbev/msx149 PMID: 28472384

21. Rambaut A, Suchard MA, XD& DA. TRACER [Internet]. 2014. Available: http://tree.bio.ed.ac.uk/software/tracer/ Accessed: 03/02/2021

22. Kass RE, Raftery AE, Bayes factors. Journal of the American Statistical Association.1995; 90: 773-795. https://doi.org/10.1080/01621459.1995.10476572

23. Gelman A, Meng XL. Simulating normalizing constants: from importance sampling to bridge sampling to path sampling. Statistical Science. 1998; 13: 163-185. https://doi.org/10.1214/ss/1028905934

24. Xie W, Lewis PO, Fan Y, Kuo L, Chen MH. Improving marginal likelihood estimation for bayesian phylogenetic model selection. Systematic Biology. 2011; 60: 150-160. https://doi.org/10.1093/sysbio/syq085 PMID: 21187451

25. Pollett S, Nelson MI, Kasper M, Tinoco Y, Simons M, Romero C, et al. Phylogeography of Influenza A (H3N2) Virus in Peru, 2010-2012. Emerg Infect Dis. 2015; 21: 1330-1338. https://doi.org/10.3201/eid2108.150084 PMID: 26196599

26. Rambaut A, Pybus OG, Nelson MI, Viboud C, Taubenberger JK, Holmes EC. The genomic and epidemiological dynamics of influenza A virus. Nature. 2008; 453: 615-9. https://doi.org/10.1038/nature06945 PMID: 18418375

27. Nelson MI, Simonsen L, Viboud C, Miller MA, Holmes EC. Phylogenetic analysis reveals the global migration of seasonal influenza A viruses. PLoS Pathog. 2007; 3: 1220-8. https://doi.org/10.1371/journal.ppat.0030131 PMID: 17941707

28. Cheng X, Tan Y, He M, Lam TT, Lu X, Viboud C, et al. Epidemiological dynamics and phylogeography of influenza virus in southern China., Infect Dis. 2013; 207: 106-14. https://doi.org/10.1093/infdis/jis526 PMID: 22930808