The Emerging Role of Rhinoviruses in Lower Respiratory Tract Infections in Children – Clinical and Molecular Epidemiological Study From Croatia, 2017–2019

Rhinoviruses (RVs) are increasingly implicated not only in mild upper respiratory tract infections, but also in more severe lower respiratory tract infections; however, little is known about species diversity and viral epidemiology of RVs among the infected children. Therefore, we investigated the rhinovirus (RV) infection prevalence over a 2-year period, compared it with prevalence patterns of other common respiratory viruses, and explored clinical and molecular epidemiology of RV infections among 590 children hospitalized with acute respiratory infection in north-western and central parts of Croatia. For respiratory virus detection, nasopharyngeal and pharyngeal flocked swabs were taken from each patient and subsequently analyzed with multiplex RT-PCR. To determine the RV species in a subset of positive children, 5′UTR in RV-positive samples has been sequenced. Nucleotide sequences of referent RV strains were retrieved by searching the database with Basic Local Alignment Tool, and used to construct alignments and phylogenetic trees using MAFFT multiple sequence alignment tool and the maximum likelihood method, respectively. In our study population RV was the most frequently detected virus, diagnosed in 197 patients (33.4%), of which 60.4% was detected as a monoinfection. Median age of RV-infected children was 2.25 years, and more than half of children infected with RV (55.8%) presented with lower respiratory tract infections. Most RV cases were detected from September to December, and all three species co-circulated during the analyzed period (2017–2019). Sequence analysis based on 5′UTR region yielded 69 distinct strains; the most prevalent was RV-C (47.4%) followed by RV-A (44.7%) and RV-B (7.9%). Most of RV-A sequences formed a distinct phylogenetic group; only strains RI/HR409-18 (along with a reference strain MF978777) clustered with RV-C strains. Strains belonging to the group C were the most diverse (41.6% identity among strains), while group B was the most conserved (71.5% identity among strains). Despite such differences in strain groups (hitherto undescribed in Croatia), clinical presentation of infected children was rather similar. Our results are consistent with newer studies that investigated the etiology of acute respiratory infections, especially those focused on children with lower respiratory tract infections, where RVs should always be considered as potentially serious pathogens.

Rhinoviruses (RVs) are increasingly implicated not only in mild upper respiratory tract infections, but also in more severe lower respiratory tract infections; however, little is known about species diversity and viral epidemiology of RVs among the infected children. Therefore, we investigated the rhinovirus (RV) infection prevalence over a 2-year period, compared it with prevalence patterns of other common respiratory viruses, and explored clinical and molecular epidemiology of RV infections among 590 children hospitalized with acute respiratory infection in north-western and central parts of Croatia. For respiratory virus detection, nasopharyngeal and pharyngeal flocked swabs were taken from each patient and subsequently analyzed with multiplex RT-PCR.
To determine the RV species in a subset of positive children, 5 UTR in RV-positive samples has been sequenced. Nucleotide sequences of referent RV strains were retrieved by searching the database with Basic Local Alignment Tool, and used to construct alignments and phylogenetic trees using MAFFT multiple sequence alignment tool and the maximum likelihood method, respectively. In our study population RV was the most frequently detected virus, diagnosed in 197 patients (33.4%), of which 60.4% was detected as a monoinfection. Median age of RV-infected children was 2.25 years, and more than half of children infected with RV (55.8%) presented with lower respiratory tract infections. Most RV cases were detected from September to December, and all three species co-circulated during the analyzed period (2017)(2018)(2019). Sequence analysis based on 5 UTR region yielded 69 distinct strains; the most prevalent was RV-C (47.4%) followed by RV-A (44.7%) and RV-B (7.9%). Most of RV-A sequences formed INTRODUCTION Rhinoviruses (RVs) are small, non-enveloped viruses that belong to the family Picornaviridae, genus Enterovirus. To date, there are 171 rhinovirus (RV) genotypes recognized and classified into the three species as RV-A (83 types), RV-B (32 types), and RV-C (56 types) (Royston and Tapparel, 2016;Pan et al., 2018). RV-A and RV-B were discovered by isolation on monkey kidney cells in 1950s (Price, 1956) while RV-C genotypes, which are not cultivable using ordinary culture methods, have been identified decades after following the rise of molecular techniques (Lamson et al., 2006;Lau et al., 2007). There are also differences between RV species in utilization of cell entry receptor: a majority of RV-A and RV-B attach to the intercellular adhesion molecule (ICAM)-1 (classified as the major receptor group) and the others alternatively bind low density lipoprotein receptor (LDL-R) (minor receptor group), whereas RV-C utilizes human cadherin-related family member 3 (CDHR3) (Bochkov et al., 2015;Royston and Tapparel, 2016).
Rhinovirus genome is a 7.2-kb single-stranded, positive-sense RNA with a single open reading frame (ORF) joined to a 5 untranslated region (5 UTR) and a short viral priming protein (VPg) (Jacobs et al., 2013). ORF encodes a poly-protein which is cleaved by virally encoded proteases in 11 proteins. Four proteins -VP1, VP2, VP3, and VP4 -make up the viral capsid and account for the virus' antigenic diversity, while the remaining non-structural proteins are involved in viral genome replication and assembly. A rather conserved 5 UTR region that harbors internal ribosomal entry site (IRES) is usually utilized for RV detection from clinical samples, while more precise genotyping is based on VP4/VP2 or VP1 sequence analysis.
Rhinovirus have been neglected for decades, primarily because they were considered less virulent and only capable of causing mild common cold, and were not recognized, until recently, as important causative agents of lower respiratory tract infections (LRTIs) and severe respiratory disease . Introduction of sensitive and technically simple molecular detection assays, especially multiplex PCR, enabled affordable RV detection in line with other common respiratory viruses in clinical samples. Together with coronaviruses, RVs are indeed responsible for majority of upper respiratory tract infection (URTI), but also for substantial rate of LRTIs in all age groups (Lauinger et al., 2013;Zhao et al., 2018Zhao et al., ,Čivljak et al., 2019. Many studies showed that RV is one of the leading causes of pneumonia, bronchiolitis and other form of severe respiratory disease (Zhao et al., 2018), standing side by side with respiratory syncytial virus (RSV) in children and influenza in elderly (Esposito et al., 2013;Ning et al., 2017;Čivljak et al., 2019). There are several studies that report increased rates of asthma exacerbations and LRTIs among children with RV-C when compared to those infected with RV-A and RV-B (Bizzintino et al., 2011;Lauinger et al., 2013); however, recent studies did not find any relationship between a specific RV species and the severity of clinical presentation (Jacobs et al., 2015;van der Linden et al., 2016;Zhao et al., 2018).
Data on RV prevalence in Croatia are scarce, and molecular epidemiology was thus far not described. This study aims to determine the RV prevalence, compare it with prevalence patterns of other common respiratory viruses, as well as to explore clinical and molecular epidemiological features of RV infections among hospitalized children with acute respiratory infection.

Ethical and Safety Issues
The research was performed in accordance with relevant guidelines/regulations and in line with the Declaration of Helsinki, as revised in 2013. Written informed consent was obtained from all participants (children's parents or their legal guardians). The study was approved by the Ethics Committee of the Dr. Andrija Štampar Teaching Institute of Public Health and conducted as part of the Croatian Science Foundation project entitled "New and neglected respiratory viruses in vulnerable groups of patients" (No. IP-2016-06-7556). All standard biosecurity and institutional safety procedures have been adhered to.

Patients and Sample Collection
From May 2017 to April 2019, a total of 590 patients were included from hospitals located in north-western and central part of Croatia: Clinical Hospital Zagreb and General Hospital Karlovac, respectively. Inclusion criteria were: age lower than 18 years, a clinical diagnosis of acute respiratory tract infection (ARI), and need for hospitalization [on ward (≥1 day) or day hospital (for more than 6 but less than 24 h)]. Exclusion criteria were presumed bacterial respiratory infection -including otitis, sinusitis and bacterial pneumonia, healthcare-associated infection, and ambulatory treated patients.
Patients were categorized into the four groups according to their respective age (i.e., <1, 1-2.99, 3-4.99 and ≥5 years of age), and two groups according to the localization of infection in those with upper respiratory tract infection (URTI) and lower respiratory tract infection (LRTI). URTI was defined by symptoms of the common cold, coryza, cough, and hoarseness with or without fever, so clinical syndromes of respiratory catarrh, rhinitis and/or pharyngitis represented URTI category. LRTI was defined by symptoms of tachypnea, wheezing, severe cough, breathlessness, and by specific clinical signs such as nasal flaring, jugular and intercostal retractions, cyanosis (in rare instances), as well as wheezing, crackles and inspiratory rhonchi or generally reduced breath sounds during auscultation. Clinical syndromes of bronchitis, bronchiolitis and pneumonia were included in LRTI category (Tregoning and Schwarze, 2010;Ljubin-Sternak et al., 2016). To avoid unnecessary X-ray exposure, chest radiographs were taken only for some of the patients in order to exclude or confirm bacterial pneumonia.
For respiratory virus detection, nasopharyngeal and pharyngeal flocked swabs from each patient were collected, combined, and placed in viral transport medium (UTM TM , Copan, Italy). Specimens were immediately transported to the Molecular microbiology laboratory at the Public Health Institute where they were stored at −80 • C until tested. The results of virology testing were released to the physicians periodically, approximately once per week. As a part of routine care, nasopharyngeal, pharyngeal swabs, and blood cultures were taken from hospitalized patients and submitted for bacterial diagnostics using standard cultivation methods. Demographic and clinical data, antimicrobial use records, and the results of routine bacterial studies were collected by a retrospective review of patient charts.

Multiplex RT-PCR
To isolate viral DNA and RNA from viral transport medium, 300 µL has been extracted according to the manufacturer's protocol using Ribospin TM vRD Kit (GeneAll Biotechnology, Seoul, Korea). Multiplex RT-PCR for 15 respiratory viruses using Seeplex R RV15 detection kit (Seegene Inc., Seoul, Korea) was performed. Briefly, multiplex PCR and cDNA synthesis as onestep reaction was performed and set up in three different tubes with three sets of primers. More specifically, "A set" contained primers for simultaneous amplification of target sequences of adenovirus (AdV), human coronavirus (HCoV) 229E/Nl63, parainfluenza virus (PIV) types 1-3, and PCR internal control to check for the presence of substances that may interfere with amplification; "B set" contained primers for HCoV OC43, RV groups A/B/C, RSV type A and B, influenza (Flu) type A, and PCR internal control; and "C set" contained primers for human bocavirus (HBoV), Flu type B, human metapneumovirus (HMPV), PIV type 4, human enterovirus (HEV), and the overall process control (human RNase P was included throughout the entire process as a control from nucleic acid extraction to amplification). Amplification was performed on thermal cycler GeneAmp R 9700 PCR System (Applied Biosystems, Foster City, United States). Detection of PCR products was done by microchip electrophoresis on the MCE R -202 MultiNA device (Shimadzu, Kyoto, Japan) -including software analysis showing results in the form of electropherograms and virtual image gels.

Rhinovirus Typing -Reverse Transcription, PCR and Sequencing
To determine RV species, 5 UTR of the genome in RV-positive samples was sequenced. Total RNA was extracted from 500 µL of RV-positive samples by the method reported by Chomczynski and Mackey (1998). Reverse transcription was performed at 42 • C for 60 min, in a reaction mix containing 10 µL of isolated RNA, 1 × PCR buffer (GE Healthcare, United Kingdom), 0.1 mM of each dNTP, 20 U of RNase inhibitor (Thermo Fisher Scientific, United States), 1.25 mM MgCl2, 2.5 mM of random hexanucleotide primers and 50 U of MuLV reverse transcriptase (Thermo Fisher Scientific, United States) in a final volume of 20 µL.
Sequencing reactions were set up with purified DNA, one of the specific primers used for amplification, and a BigDye Terminator v3.1 Cycle Sequencing Kit (Thermo Fisher Scientific, United States) according to the manufacturer's protocol. Sequencing and sequence analysis were performed on a 3130 Genetic Analyzer (Thermo Fisher Scientific, United States).

Phylogenetic Analysis
Nucleotide sequences of referent RV strains were retrieved by searching the database with BLAST (Basic Local Alignment Tool) 1 and used to construct alignments and phylogenetic trees. Alignments were performed using MAFFT multiple sequence alignment tool available at the EMBL-EBI website 2 , and edited in AliView v1.23 (Larsson, 2014). Phylogenetic trees were generated using the maximum likelihood method with Molecular Evolutionary Genetics Analyses (MEGA) software v6.06 (Tamura et al., 2013), under the most appropriate model of nt substitution determined with jModeltest v2.1.4 (Darriba et al., 2012). Bootstrap probabilities for 1,000 iterations were calculated to evaluate confidence estimates. Sequence conservation (defined as percentage of genomic positions identical in all strains; gaps were ignored during computing) and evolutionary distances (p-distances) within and between groups have been calculated using MEGA 6.06. The sequences of HRSV strains obtained in this study were deposited in the GenBank under acc. nos. MN369460 -MN369528.

Statistical Analysis
Data analysis was performed using Stata/MP (Ver.15.1; StataCorp LLC, College Station, United States). Age was presented using medians with stated interquartile range (IQR). Demographic and clinical parameters were compared by χ 2 -test or Fisher's exact test for categorical variables and by Kruskal-Wallis test for continuous variables. To assess the strength of association between dependent variable (RV positive patients) and age we used univariate logistic regression. P-value was set to <0.05.
According to the clinical presentation, there was a significant difference in the proportion of LRTIs between the type of the respiratory viruses (P = 0.002) (Figure 2). More than half of children infected with RV (110; 55.8%) presented with LRTI; nonetheless, this was not significantly different from the proportion of RV positive children with URTI (87; 44.2%) (P = 0.336), while children with RSV infection significantly more often presented with LRTI (P < 0.001).
Of the three respective groups, strains belonging to the group C were the most diverse, with identity of 41.6% (142 of 341 identical positions), while group B was the most conserved with 71.5% identity among strains (241/337). Group A comprised strains which shared an overall 54.3% identical positions (183/337). Calculated p-distances between groups showed group A and C are more closely related (p-distance 0.234) than to group B. Similar p-value was calculated between group B and group A (p-distance 0.31) or group C strains (p-distance 0.33), respectively.
No significant difference has been demonstrated in clinical symptoms based on the RV species, with the exception of increased frequency of antibiotic treatment in those infected with RV-A species (P = 0.012) ( Table 2). Most RV cases were detected from September to December, and all three species co-circulated during the analyzed period (Figure 4).

DISCUSSION
In this study we initially evaluated the prevalence of RV and other common respiratory viruses, as well as RV species distribution in hospitalized children with symptoms of ARI over a period of 2 years. Our results are consistent with previously published studies that investigated the etiology of ARI, especially those focused on children with LRTIs (Chen et al., 2015;Ning et al., 2017). Indeed, RV is right next to RSV when addressing the most common causes of bronchiolitis in hospitalized children, as demonstrated in very low-birth-weight infants from Argentina and in one multicentre prospective study from the United States, respectively (Mansbach et al., 2012;Miller et al., 2012). Although RV can be found both in upper and lower respiratory tract, the potential for spread and the pathophysiology of infection in those two regions differs (as evidenced by studies conducted on RSV) (Kim et al., 2016;González-Parra and Dobrovolny, 2019). Possible causes of such disparity are fundamental differences in the immune responses and virus-cell interactions between these two anatomical regions, resulting in altered disease manifestation and spread (González-Parra and Dobrovolny, 2019). There is also evidence that mucus velocity is decreased in small children (akin to elderly individuals) (Grubb et al., 2016), making them more susceptible to LRTI and creating in turn a potential niche for more detrimental effect of RV infection.
When age is concerned, our study has showed that RV holds a middle ground, not affecting very young children as is the case with RSV. This is comparable with the results found in the study by Cebey-López et al. (2015), where RSV infection was also seen in the youngest age groups (less than 1 year of age), RV was present in somewhat older children (mean between 14.4 and 40.9 months), while AdV predominated in patients older than 50 months. Conversely, some other author groups pointed out how RV can predominate in practically all age groups (Tsagarakis et al., 2018).
In order to determine the RV species, we decided to sequence and analyze 5 UTR, which has been established as a relatively simple and rapid technique for identifying the RV serotypes in clinical samples. Moreover, it has been shown that 5 UTR RT-PCR demonstrated greater sensitivity than VP4-VP2 PCR, as reflected by the higher positivity rate in amplification of clinical isolates, and further, no need for a nested PCR or multiple primer pairs -reducing in turn the contamination rate, cost and turnaround time (Kiang et al., 2008). Despite using primers which target the region widely used for typing purposes, more than half of the samples produced no amplicons suitable for sequencing. Such low rate of successful RV sequencing from clinical samples may be a consequence of mutations in primer regions or low viral  load in original sample, but most probably arises due to RNA degradation during freezing/thawing cycles. Phylogenetic groups were readily distinguished and all three RV groups were detected, with group C and A predominating. The proportion of the three RV species revealed in this study (RV-A 44.7%, RV-B 7.9%, and RV-C, 47.4%) is consistent with prior studies worldwide (RV-A, 35.9-67.7%; RV-B, 1.5-13%; RV-C, 23-59.3%) (Lauinger et al., 2013;Rahamat-Langendoen et al., 2013;Jacobs et al., 2015;Tsatsral et al., 2015;Ratnamohan et al., 2016;van der Linden et al., 2016;Zhao et al., 2018).
Intermixing of groups was not observed, except for a single group A sequence (RI/HR409-18) which clustered with group C sequences, indicating a recombination event or co-infection as was proposed by Richter et al. (2015). Furthermore, group B strains were detected sporadically during the analyzed period, which is in accordance with other studies (Lauinger et al., 2013;Launes et al., 2015;Richter et al., 2015). These results represent the first report on RV diversity in Croatia.
In previous reports RV-C (and in lesser extent RV-A) have been associated with more severe illness (Bizzintino et al., 2011;Lauinger et al., 2013;Linder et al., 2013;Chen et al., 2015); however, more recent studies failed to report the connections between species and disease severity (Rahamat-Langendoen et al., 2013;Jacobs et al., 2015;van der Linden et al., 2016;Zhao et al., 2018). In this study there were no significant differences observed in clinical symptoms among three species, except more frequent utilization of antibiotic therapy in patients with RV-A species which can be result of subjective clinical assessment of more severe disease and empirical introduction of therapy.
RV circulated throughout the 2-year period covered by this study with peaks in autumn and winter months. Previous studies reported that RV infections occur all year round, with peaks of infection usually in spring and autumn months (Čivljak et al., 2019). However, recent research endeavors also report peaks in autumn and winter months (Zhao et al., 2018), and some of them specifically note that RV-C demonstrate peak in winter months (Linder et al., 2013). Our study also observed no RV-C detection in spring and summer season.
There are several limitations of the study. Due to crosssectional study approach, a control group of asymptomatic patients could not have been included (which would facilitate assessment of RV infection severity). Furthermore, we performed only 5 UTR targeted RT-PCR assay and did not confirm our results with VP4/VP2 or VP1 sequences analysis. Some authors report discordance between proposed phylogeny groups when sequences from the 5'UTR and VP4/VP2 coding regions were analyzed, which can result in imprecise classification (Ratnamohan et al., 2016). More specifically, samples that clustered as RV A using 5 UTR analysis can be revealed as RV-C when VP4/VP2 region is analyzed (Ratnamohan et al., 2016). Also, complicated rhinovirus infections were excluded from the study. Notwithstanding the aforementioned limitations, this is the first study endeavor cataloging a circulation of RV species in Croatia over 2 years.

CONCLUSION
In conclusion, in our study more than half of children infected with RV presented with LRTI, which (together with other newer studies in the field) underlines the need for paradigm shift where RVs will not be merely associated with URTI, but also considered in cases of LRTI. Regardless of the diversity of RV found in this study and the purported heterogeneity of the RV strains infecting the children, the similarity of clinical presentation negates the notion that certain RV species might be more virulent, at least in our case. More clinical and epidemiological studies are warranted to further elucidate this issue, with inevitable use of animal models to study pathogenic specificities of RV infection.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in the GenBank-accession numbers from MN369460 to MN369528.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Ethics Committee of the Dr. Andrija Štampar Teaching Institute of Public Health and conducted as part of the Croatian Science Foundation project entitled "New and neglected respiratory viruses in vulnerable groups of patients" (No. IP-2016-06-7556). Written informed consent to participate in this study was provided by the participants' legal guardian/ next of kin.

AUTHOR CONTRIBUTIONS
SLJ-S, II-J, and JV designed the research. SLJ-S, AS, and DF performed the experiments. MM and TT collected the data. BK, AS, and DF analyzed the data. SLJ-S, TM, MM, TT, II-J, AS, and DF interpreted the data and prepared the draft of the manuscript. JV and TM critically reviewed the draft. SLJ-S and TM wrote the final version of the manuscript. All authors reviewed and approved the final version of the manuscript.

FUNDING
This work has been fully supported by the Croatian Science Foundation under the project No. IP-2016-06-7556 titled "New and neglected respiratory viruses in vulnerable groups of patients" (principal investigator SLJ-S). The funders had no role in the study design, data collection and analysis, decision to publish or preparation of the manuscript.