Human rhinovirus spatial-temporal epidemiology in rural coastal Kenya, 2015-2016, observed through outpatient surveillance [version 2; peer review: 2 approved]

Background: Human rhinovirus (HRV) is the predominant cause of upper respiratory tract infections, resulting in a significant public health burden. The virus circulates as many different types (168), each generating strong homologous, but weak heterotypic, immunity. The influence of these features on transmission patterns of HRV in the community is understudied. Methods: Nasopharyngeal swabs were collected from patients with symptoms of acute respiratory infection (ARI) at nine out-patient facilities across a Health and Demographic Surveillance System between December 2015 and November 2016. HRV was diagnosed by real-time RT-PCR, and the VP4/VP2 genomic region of the positive samples sequenced. Phylogenetic analysis was used to determine the HRV types. Classification models and G-test statistic were used to investigate HRV type spatial distribution. Demographic characteristics and clinical features of ARI were also compared. Results: Of 5,744 NPS samples collected, HRV was detected in 1057 (18.4%), of which 817 (77.3%) were successfully sequenced. HRV species A, B and C were identified in 360 (44.1%), 67 (8.2%) and 390 (47.7%) samples, respectively. In total, 87 types were determined: 39, 10 and 38 occurred within species A, B and C, respectively. HRV types presented heterogeneous temporal patterns of persistence. Spatially, identical types occurred over a wide distance at similar times, but there was statistically significant evidence for clustering of types between health facilities in close proximity or linked by major road networks. Open Peer Review


Introduction
Human rhinovirus (HRV) is the predominant cause of upper respiratory tract infections (URTIs) referred to as the common cold [1][2][3][4] . The virus has also been associated with lower respiratory tract illnesses, including exacerbation of asthma and chronic obstructive pulmonary disease 5,6 . HRVs have also been identified in mild and asymptomatic cases 7 . HRV is globally ubiquitous and infections usually occur all year-round, although peaking in the early autumn and late spring in many temperate or subtropical countries, and in the rainy season in tropical countries [8][9][10] . Despite most HRV infections being mild, they pose a considerable social and economic burden due to time lost from work or school, medical attendance and reduced performance of regular duties 11 . HRV transmission, infection patterns and diversity have been rarely studied in low-income settings despite bearing the majority burden of acute respiratory illnesses (ARI).
HRVs fall under the genus Enterovirus (HEVs) in the family Picornaviridae 12 . Their genome occurs as a positive-sense, single-stranded RNA molecule of approximately 7.2 kb flanked by a 5' untranslated region (UTR) and 3' poly-A tail. There are three HRV species denoted as HRV-A, HRV-B and HRV-C within which over 168 genetically distinct types have been identified 13 . Identification of these types presently relies on molecular typing methods based on sequence analysis of VP1 or VP4/VP2 proteins encoding regions 12,14 .
Studies investigating the transmission of HRV show that multiple infections frequently occur in individuals over short periods of time of the order of a few weeks [15][16][17] . In household studies, family members experience rapid rates of reinfection (for instance, up to 5 infections in adults or 12 infections in young children per year) where each infection tends to be caused by a different HRV-type and rarely with the same HRV-type 4,18-20 . At a population level of observation, repeated reinfections in individuals appears to reflect the introduction and circulation of new types in the community, rather than persistence of the same HRV types [21][22][23][24] . These observations are consistent with evidence from immunological studies on HRV that indicate strong homologous type responses and little heterotypic immunity to protect individuals against multiple types [25][26][27] . Moreover, studies have shown that individuals develop relatively long-lasting (~ 1 year) type-specific neutralizing antibodies (IgG and IgA) against a specific infecting HRV-type but susceptibility to other different HRV types remains and infection will occur as so long as there is an exposure 28,29 .
Hence it might be hypothesized that at the population level the virus types spread in a manner independent of each other. However, the statistical evidence to support this assertion has yet to be undertaken. This would require spatially structured study designs to capture temporal-spatial HRV transmission. Data of this sort could improve understanding on the nature of spread of these viruses at the population level in relation to homotypic and heterotypic immunity.
A recently published study 30 , described the epidemiology of respiratory viruses in coastal Kenya through surveillance of acute respiratory infection presentations at outpatient facilities spatially structured across the well-defined population of the Kilifi Health and Demographic Surveillance System (KHDSS) 31 . The present study utilized HRV positive samples from this surveillance between December 2015 through November 2016, to explore the temporal and spatial circulation patterns and genetic diversity of HRV types in the wider Kilifi community.

Study area
The study was undertaken within Kilifi County in Coastal Kenya, at nine health facilities namely: Matsangoni, Ngerenya, Mtondia, Sokoke, Mavueni, Jaribuni, Chasimba, Pingilikani and Junju, all located within the KHDSS as previously described 30 . The county has a population of about 1,109,735 and is predominantly rural; the main economic activities include fishing, tourism and subsistence farming. The area has a tropical climate with two rainy seasons (long rains in April to July and short rains in October to December). The KHDSS area in which the study was restricted covers 891km 2 and has a population of 287,014 (mid-point 2016) which is heterogeneously distributed with the highest density in Kilifi township 30 .

Study design
Nasopharyngeal swab (NPS) samples and demographic data were collected from patients presenting with ARI symptoms at the nine selected health facilities within the KHDSS between December 2015 and November 2016 inclusive. The KHDSS area has 21 public health facilities operating under the Kenya Ministry of Health (MOH); the 9 health facilities were purposively selected

Amendments from Version 1
Revisions have been made to the manuscript based on reviewers' comments. The major changes are as follows: • In the Laboratory procedures section, we have added text to the manuscript to clarify on our genotyping assay.
• In the Results section, we have added information to highlight there was no significant difference between species A and C (predominant species) for most clinical signs reported, except for nasal discharge (p<0.004) and difficulty in breathing (p<0.019). In addition, we have also added a column on Supplementary

Sequence alignment and phylogenetic analysis
Raw sequence reads were assembled into contigs using Sequencher software (version 5, Gene Codes Corporation, Ann Arbor, USA). Multiple sequence alignments (MSA) were prepared using MAFFT v7.220 36 . Phylogenetic trees were constructed in MEGA v.6.0 37 with maximum likelihood methods under the GTR model and branch support was assessed using 1000 bootstrap iterations. Types were assigned based on >90% nucleotide similarity to rhinovirus prototype sequences (also referred to as reference sequences) 13,14 and phylogenetic clustering with the reference sequences (with a bootstrap support value above 70%).

Statistical analysis
Initial statistical analysis was conducted in STATA version 13.1 (College Station, Texas). Statistical comparison of demographic characteristics and clinical features in HRV species was conducted using chi square test of association and Fisher's exact test where applicable. Frequency distribution and temporal pattern graphs for HRV were generated in R 3.5.1.

Spatial distribution analysis
Analysis centered on investigating spatial heterogeneity of HRV detections of HRV types, between health facilities. Our null hypothesis was that HRV is introduced and spreads randomly within a population, and that HRV types circulate independent of each other. Hence, we asked: a) Is there evidence that any health facility had significantly more HRV detections in the sample of NPSs than others? b) Are the frequencies of HRV types found per health facility noticeably different from independent random mixing? c) Moreover, was the pattern of types detected at each facility similar to nearby health facilities compared to further away health facilities? To answer these questions we used a combination of the classification algorithms provided by the Orange data mining toolbox 38 and logistic regression implemented by the MATLAB v9.4 fitglme function. The null hypothesis of independent random HRV type detection was tested using the G-test statistic associated with the contingency table of type frequency per health facility. Type pattern similarity between health facilities was analysed using multi-dimensional scaling (MDS) implement in Orange v3.0. MDS represents the similarity between high-dimensional data-points optimally on a low dimensional (in this case two dimensional) Euclidean metric space, thereby revealing patterns and hidden dimensions in the data 39 . The Orange implementation of MDS includes a pairwise similarity network linking all points whose high-dimensional similarity is close compared to the typical similarity between any two data-points; this is a threshold quantity set to the algorithm's maximum therefore including all pairwise connections inferred by the Orange MDS algorithm in our analysis.
Nucleotide sequence data set All newly obtained sequences in this study have been deposited into GenBank under accession numbers MH459421-MH460237.

Results
Between December 2015 and November 2016 inclusive, a total of 5750 participants were recruited and NPS samples collected from 5744 (99.9%) participants as shown in Table 1. Of the 5744 samples collected, 5741 (99.9%) were tested for HRV. HRV was detected in 1057 (18.5%) of the samples tested.
Of the 1057, 817 (77.3%) were successfully amplified and sequenced for the VP4/VP2 genomic region; the remaining 240 (22.7%) samples either totally failed amplification or had short consensus sequences recovered (<300 nucleotides) and were not included in phylogenetic analysis.
Among the HRV positive cases, the median age was 2.0 years, 64.0% were under 5 years of age, and 53.9% were females ( Table 2). The percentage of samples found to be of species A, B and C was 44.1 % ( 360), 8.2% (67) and 47.7% (390), respectively, and did not differ by age or sex (Table 2). Cough (94.5%) and nasal discharge (79.9%) were the most common clinical presentations. Difficulty in breathing, chest in-drawing, crackles, lethargy, nasal flaring, sore throat, wheezing were also recorded ( Table 2). Clinical presentations differed significantly between the typed HRV samples for nasal discharge (p< 0.012) and difficult breathing (p< 0.022) There was no significant (p < 0.05) difference between species A and C (predominant species) for most of the clinical signs reported, except for nasal discharge (p< 0.0042) and difficult breathing (p<0.0194) (Supplementary Table 1, Supplementary File 1).
Phylogenetic analysis revealed 3 major clades representing the HRV-A, B and C species (Supplementary Figure 1). Three samples were not assigned to any HRV type and these were identified as Enterovirus D68 (n=2) and Coxsackievirus B5 (n=1, see offshoot in Supplementary Figure 1). A total of 87 different HRV types were identified: 39 within species HRV-A, 10 within HRV-B and 38 within HRV-C. Within HRV-A, A15 was the predominant type (n=62), followed by A58 (n=36) and A41 (n=20). Other HRV-A types occurred at lower frequencies, ranging from 1 to 19 cases ( Figure 1 panel A). Within HRV-B, B35 presented as the dominant type (n=41) with the other HRV-B types being detected in 1 to 8 cases ( Figure 1 panel B). Within HRV-C, C22 type was the predominant type (n=48), followed closely by C11 (n=45) and C38 (n=35). The other HRV-C types were detected as shown in Figure 1 panel C.
HRV circulated throughout the study period as shown in Figure 2A. All HRV species were detected in all months except in December 2015 when only A and C were detected, although fewer samples were collected in this month. HRV incidence peaks were observed in August, September and October, while troughs were observed in December 2015 and January 2016, which was likely due to fewer samples collected in the first two months of the study. Despite HRV-C being the most prevalent, there was a varying dominance between HRV-A and HRV-C with time ( Figure 2B). HRV species A appeared dominant in December,
Temporal patterns for each of the most frequently occurring types (with cases >16 during the year) are displayed in Figure 3. Most types displayed a unimodal distribution, with peak occurrence including August to October. However, there was considerable variation in the spread of the seasonal occurrence and modal month.
Regression analysis showed that the month of sample collection and age of participant were the most informative predictors for determining whether an individual tested positive for any HRV type. These predictors were identified whether undetermined samples were used in the regression analysis or not. Other factors considered were gender and health facility. First, five standard classification models (logistic regression, neural network using logistic activation functions and one hidden layer, k-nearest neighbor classification, tree and random forest regression) were compared using random partitions of the data into 2/3 training set, 1/3 testing set. Over 100-fold repetitions of training and testing each model, both logistic regression and the neural network classifiers had 77.2% classification accuracy (this increased to 81.0% and 80.9% respectively if undetermined samples were removed from the data set). Other classification models were less successful and because logistic regression is easier to interpret, and likelihood-based, we chose to continue our analysis using only logistic regression. Second, we removed uninformative predictor variables by comparing Akaike information criterion (AIC) for each additive combination of predictor variables and selecting the logistic regression model with lowest AIC. Whether undetermined samples were included in the regression analysis or not, in either case the lowest (i.e. best) AIC model included only month of sample collection and the age of the sample participant as predictors (Supplementary  Table 2, Supplementary File 1). The most informative models implied that the odds ratio of the participants being detected HRV positive declined by 0.1-0.2% per month of life, and that collection in August to October increased the odds ratio of finding HRV positives (Supplementary Table 3, Supplementary  File 1).
Notably, the health facility of sample collection was not included amongst the predictor variables for our classification analysis of HRV positive samples. However, this does not exclude the possibility of correlations in the types collected between different regions. Most sequences typed as belonging to the same HRV type were from different regions (represented by different health facilities) of the KHDSS and some were detected in several months. On the basis of phylogenetic analysis, closely related HRV-A15, C11, C22, C28 types were present in all 9 regions of KHDSS, and A66, B35, C14 and C44 types detected in 8 different regions. In other cases, identical HRV-types were shared in the different geographical regions ranging from 1 to 7 geographical regions. Nonetheless, the contingency table of type collected at each health facility (Supplementary Table 4, Supplementary File 1) reveals highly significant evidence to reject independence of type occurrence at each health facility (G statistic = 849.7, dof = 704, P = 1.25×10 -4 ) in favour of clustering of type occurrence by health facility.
In addition to clustering of individual type occurrence at the health facilities we also found evidence of greater similarity in HRV type distribution between the health facilities located near the road network in central Kilifi area (Chasimba, Jaribuni, Mavueni, Mtondia, Ngerenya, Sokoke) compared to the health facilities off the main road in the south of the KHDSS area (Junju, Pingilikani) and further north up the main highway (Matsangoni) ( Figure 4A). We measured the pairwise similarity of the type distribution at each health facility as the total absolute differences in their type occurrences shifted by the data median and rescaled by the median absolute distance from median (i.e. the rescaled Manhattan metric recommended for high dimensional data sets 40 ). We used the Orange implementation of MDS (see methods) to optimally represent the type distribution similarities on a plane along with a similarity network ( Figure 4B). The similarity network distinguishes between the type distributions of health facilities in the central part of Kilifi area on the road network, in particular the central similarity clique (Mtondia, Ngerenya, Sokoke), and the type distribution observed on the outskirts and harder to travel to parts of the Kilifi area. The MDS representation condenses and quantifies the information derived by direct visual comparison of the type distributions ( Figure 4C).

Discussion
This study investigated the HRV infection, diversity and type distribution in individuals presenting with ARI symptoms for outpatient care at nine spatially structured health facilities in rural coastal Kenya over 12 months period. All three HRV species co-circulated, with HRV-A and HRV-C co-predominant. This observation agrees with a previous study conducted in rural coastal Kenya 41 . Remarkably, a total of 87 HRV types were identified in circulation over the twelve months, representing 51.7% of all known HRV types. Occurrence of HRV throughout the year appears to be sustained by the existence of simultaneous and successive mini-epidemics each caused by a different HRV type introduced into the community independently. We assume each HRV type generates strong and lasting homotypic immunity, that leads to type-specific herd immunity and subsequent fadeout. This observation could be attributed to the natural ability of each type to independently cause an infection with limited cross-immunity with other types 42 . Therefore, as long as there are new introductions and exposure to new HRV-types, the population experiences a series of new HRV infections throughout the year 21-23 .
As expected, there was a decrease in the detection of HRV cases with an increase in the age of the patients. The majority of the HRV positive patients were children below the age of 5 years (63.95%) and a statistically significant decrease in HRV positives amongst older ages was observed. The consistent finding of high proportions of HRV infections in children below the age of 5 years reiterates the need to focus control strategies to this age group since they are the most vulnerable, and presumably make the greatest contribution to community transmission. The low rates of detection of HRV types in adults may reflect the gradual development of type-specific immunity due to increased exposure throughout life.
We identified two Enterovirus D68s and one Coxsackievirus B5 in the nasopharyngeal samples initially classified as HRV positive using a real-time RT-PCR method. This is not unusual as the real time RT-PCR target region is genetically close between HRV and other members of the genus Enterovirus 35 . A previous study conducted in Kilifi using the same molecular diagnostic assay reported a similar detection of non-HRV enteroviruses 41 . Further, a recent study in Tanzania observed a relatively high prevalence of non-HRV enteroviruses in NPS sample including poliovirus type 1, enterovirus-D68, A71, echovirus-6, 7, 9 , 11 and a variety of coxsackievirus serotypes 43 .
Phylogenetic analysis revealed close genetic association between sequences from different health facilities, with high intra-type genetic identities between the sequences from different health facilities (87 -100%). In some cases, the identical HRV-types were circulating simultaneously in different geographical regions separated by a distance as far as 30 km apart. Supported by the low genetic variation in the VP4/VP2 coding sequences and the concurrent distribution between these identical HRV-types it is probable that there is a rapid spread of HRV-types within the KHDSS once introduced and or many introductions of the same type into different areas of the KHDSS.
Despite phylogenetic evidence for rapid spread across the KHDSS, there is highly significant evidence for variation in the distribution of HRV-types between health facilities. An obvious mechanism that accounts for the variation in HRV-type distribution between health facilities is that transmission occurs more frequently between people attending the same health facility compared to those who attend other health facilities. In this case, sharing a health facility is a proxy for being more likely to live nearby and share the same social amenities and gatherings which are hotspot for transmission of respiratory viruses. Additionally, multidimensional scaling of the type distributions reveals greater type distribution similarity between the health facilities servicing parts of central Kilifi area that are more easily travelled between by road compared to the health facilities located significantly off the highway, or much further north along the highway. Spatial differentiation of type distribution, along with greater similarity in type distribution between the areas where we expect more human co-transit, is consistent with the expectations of spatial metapopulation models of infectious disease transmission 44,45 . However, it is comparatively rare to be able to demonstrate direct evidence of metapopulation effects in transmission, as we have done (see the discussion in Grenfell et al. 46 ). We ascribed (above) the observation of multiple mini-epidemics in this community, each caused by a different HRV type, to the generation of homotypic herd immunity. If there was stronger heterotypic immunity, we would expect fewer HRV cases, and even periodic HRV fadeout. With strong homotypic and weak heterotypic immunity, we therefore expect newly introduced HRV types to spread unconstrained by pre-existing heterotypic immunity in the population. The differences in circulation periods observed among HRV-types could be as a result of stochastic differences in frequency of introduction and onward transmission of the HRV types, or variation in type-specific immunity. Similar observation has been reported in 47 where simultaneous and successive epidemics caused by different HRV types contributed to HRV high incidence and enabled HRV to remain in the local population for extended periods. In some cases, type-specific epidemics were served with different variants of the same HRV type probably as a result of separate introductions into the community at different times over the year. Moreover, comparing the HRV incidence with other respiratory viruses, as shown in the preceding study 30 , HRV tend to peak in the second half of the year as other respiratory viruses report low incidence during this period. An explanation for this observation could be a case of virus interactions 48,49 that whatever triggers (either environment or biological) high incidence in HRV in the 3 rd quarter of the year, causes low incidence in the other viruses. Further studies are needed to evaluate whether virus interaction affects incidence and transmission of individual viruses in the population.
Although the study has strength in its structured design and its implementation, it faced a number of limitations. First, there was low patient recruitment in the initial study months (December 2015 and January 2016) due to clinic closure on public holidays or people migrating to other areas of the country during festive season. This may have contributed to the lower observed HRV prevalence compared to the other months of study. Second, we failed to sequence the VP4/VP2 coding region for 22.2% of HRV positives samples, possibly due to low virus titers in the nasopharyngeal samples inferred from high Ct values (Ct value>33) in these samples. Third, since a maximum sample number per clinic was fixed, changes in prevalence of one virus, for example during a seasonal peak of coronavirus, would lead to fewer samples testing positive for other viruses. Hence there would be lack of independence of numbers of each virus type, or even between HRV types, over time, affecting temporal patterns of absolute numbers. However, the prevalence of each virus or virus type at any timepoint will reflect that circulating in the community.
In summary, we observed co-circulation of the three HRV species in all nine health facilities scattered in the KHDSS area of coastal Kenya, and combined we documented the circulation of a majority of all known HRV types to date (87/165) within a single year period in this small geographical area. Some of the HRV-types circulated in the KHDSS population close to all months of observation (10/12) suggesting marked local persistence of some types while others appeared and faded from circulation quite rapidly possibly due to low herd homologous immunity for the former and stronger herd homologous immunity in the latter. HRV transmission in the community is enhanced in people living close to one another and between areas linked by road network. Our study reports a substantial HRV burden among patients seeking outpatient care in this low-income setting of tropical Africa and a differential prevalence of the HRV species and types with significant differences in their local spatial-temporal distribution.

Data availability
The replication data and analysis scripts for this manuscript are available from the Harvard Dataverse: https://doi.org/10.7910/ DVN/DUQNDX 50 .
As the dataset contains potentially identifying information on participants, it is stored under restricted access. Details on the eligibility for access and request form are available from http://kemri-wellcome.org/about-us/#ChildVerticalTab_15 for consideration by our Data Governance Committee (dgc@kemriwellcome.org)

Grant information
This work was supported by the Wellcome Trust [Grant No. 102975] and the DELTAS Africa Initiative . The DELTAS Africa Initiative is an independent funding scheme of the African Academy of Sciences (AAS)'s Alliance for Accelerating Excellence in Science in Africa (AESA) and supported by the New Partnership for Africa's Development Planning and Coordinating Agency (NEPAD Agency) with funding from the Wellcome Trust [107769/Z/10/Z] and the UK government. The views expressed in this publication are those of the author(s) and not necessarily those of AAS, NEPAD Agency, Wellcome Trust or the UK government.
The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Open Peer Review
showed a high prevalence and diversity of HRV in out-patient with symptoms of acute respiratory infection.
Globally, the paper is very interesting because it allows not only to know the different serotypes of rhinovirus circulating in Kenya but also their spatial-temporal circulation. The study is also a contribution to the scarcity of data observed in low-income settings regarding HRV circulation and diversity.
Some remarks and comments about the study are noted: Study design: The case definition used should be more specific and explain why children under 7 are excluded from the study. A more appropriate reference for acute respiratory infections should be used (WHO Guideline describing standards for influenza surveillance including case definitions for ILI and SARI 1 ).
Authors should specify if the sample is a nasopharyngeal or nasal specimen swab. The cited technique appears to be a nasal swab and not a nasopharyngeal swab.

Laboratory procedures
Why are samples with Ct> 35 considered negative? This was not explained in reference 27, which the author refers to us.
The multiplex real-time PCR assay system used for the detection of rhinoviruses does not differentiate between rhinoviruses and enteroviruses. Other enteroviruses such as Cox B5 and EV-D68 have been found, and this may lead to overestimation of the prevalence of rhinoviruses

Discussion
Remarkably, a total of 87 HRV types were identified in circulation over the twelve months, representing 52.7% of all known HRV types: Could the authors provide the total number of HRV types to better understand how the 52.7% were obtained Other members of the genus Enterovirus was identified: the authors could discuss the use of more specific detection methods for rhinovirus than to mention similar studies. This part must be more discussed.
The case definition used should be more specific and explain why children under 7 are excluded from the study. A more appropriate reference for acute respiratory infections should be used (WHO Guideline describing standards for influenza surveillance including case definitions for ILI and SARI1).

Response
All details of eligibility and sample collection are described in the baseline report on the study (Nyiro et al. 2018). Individuals of any age were eligible excepting those aged less than 7 days. This exclusion was made since it was believed few children less than 7 days of age would present with respiratory virus infections and considered that subjecting the child to a deep NPS was not justified. We did not use ILI or SARI definitions, instead we used cough, sneezing, nasal congestion, difficulty in breathing, or fast breathing for age (the latter as defined by WHO guidelines on patient management of common childhood illnesses (2013).
Authors should specify if the sample is a nasopharyngeal or nasal specimen swab. The cited technique appears to be a nasal swab and not a nasopharyngeal swab.

Response
A deep nasopharyngeal swab was collected from participants, as described in Nyiro et al 2018. For clarity we have added the text "to a distance where the tip located the deep nasopharynx," in the last paragraph of Study design.
Why are samples with Ct> 35 considered negative? This was not explained in reference 27, which the author refers to us.

Response
The reviewer is correct that the reason for considering as negative samples with Ct>=35.0 was not explained in ref 27. We apologize for this oversight. The cut-off value of Ct 35.0 was given in the paper describing the protocol that this study followed (Hammitt et al. 2011), but the rationale for setting the cut-off at 35.0 was not given. The reference for this paper has now been added to the manuscript on Laboratory procedures section. In addition, experience has shown us that VP4/VP2 RT-PCR often (95%) fails for samples with Ct > 35.0.

Response
There is potential for over-estimation of rhinovirus prevalence for the reason stated.HRV and other enteroviruses are closely related at the genetic level. However, we sequenced in the VP2/VP4 region 78% of the rhinovirus positives (817/1057), identifying only 3 nonrhinovirus enteroviruses. This suggests that the over-estimation of prevalence is small.
Remarkably, a total of 87 HRV types were identified in circulation over the twelve months, representing 52.7% of all known HRV types: Could the authors provide the total number of HRV types to better understand how the 52.7% were obtained Response Current HRV classification method puts the total number of HRV types at 168, www.picornaviridae.com, (accessed on 13 February 2019). Our study identified and typed 87 HRV types: these represent 51.7% of the known HRVs. This has been discussed in the first paragraph of Study design.
Other members of the genus Enterovirus was identified: the authors could discuss the use of more specific detection methods for rhinovirus than to mention similar studies. This part must be more discussed.

Response
We are not entirely sure what the reviewer is requesting. What we can say is that the study used a sensitive real-time PCR for the detection of HRV. Further that the assay targeted the 5' UTR which is genetically similar for both rhinoviruses and other enteroviruses. As a result, non-rhinovirus enteroviruses were detected. However, these were uncommon. Using primers that anneal to only the highly conserved region on the 5 'UTR of all the HRV species can provide the most sensitive assay for the detection of HRV strains (Bochkov et al. 2014) The authors mention a predominance of rhinovirus C followed by rhinovirus A: authors can discuss this predominance by analyzing the clinical signs of the study population. It is known that rhino C is more prevalent in asthmatics.

Response
Based on test of proportions assuming equality, there was no significant (p < 0.05) difference between species A and C (predominant species) for the clinical signs reported except for nasal discharge (p< 0.0042) and difficult breathing (p<0.0194), which differed significantly. We have added this information on the second paragraph of Results section and also inserted a column on the same in supplementary Table 1 When discussing occurrence of multiple typical mini-epidemics caused but specific types in the community, the authors should refer to other studies elsewhere and add these references to the discussion.

Response
We are aware of one other study showing multiple typical mini-epidemics (Sansone et al. 2013).The study demonstrates that simultaneous and successive epidemics caused by different HRV types contributed to high HRV incidence and enabled HRV to remain in the local population for extended periods. We have added text in the Discussion section (page14) to this effect.
Partly epidemiology of HRV species among the populations have revealed substantial diversity of HRV types in this region (L'Huillier et al. 2015) (Milanoi et al. 2016). The present study not only identifies the HRV types in circulation but describes their spatial-temporal patterns in the community. These data provide added perspectives on the spread patterns of HRV types in a community.
Could the authors add some information about why they excluded infants younger than 7 days and people with ARI symptoms for longer than 30 days?

Response
We provided an explanation for the 7day exclusion in response to question 3 by the first reviewer. We felt that the inclusion of this age group would not significantly increase cases (of all respiratory viruses) and hence considered the collection of a deep NPS could not be justified. Viral ARI are self-limiting, often lasting between 7-14 days (Pappas et al. 2008) for immune-competent individuals. For persistent ARI (>14 days), the individual could be having other underlying chronic infections e.g. TB or immune compromised due HIV-AIDS or malignancy (Chu et al. 2016) (Piralla et al. 2015).
with Ct values >35.0 (i.e. for samples with low viral load).
Did the evidence 300nt was detrimental to calling genotypes? For example, they use truncated longer amplicons and the outcome by BLAST to determine a real difference?

Response
Indeed, shorter sequences were detrimental to calling genotypes, we had to set 300nt of the VP4/VP2 region as the minimum nucleotide size to call HRV types because, with reduction of nucleotide size a typed sequence tends to shift to other branches of genotype/type or fails to cluster with a reference sequence thus falling as untypable.
The authors might consider the term "species" to describe the varying dominance between HRV species discussed around Figure 2.

Response
We are not entirely clear what the reviewer is suggesting. However, we have now added "species" in the text describing the varying dominance between HRV species-A and C in the first paragraph on page 10 of the manuscript.

Response
The circulation of HRV throughout the year appears to be maintained by the existence of simultaneous and successive mini-epidemics each caused by a different HRV type introduced into the community independently. We assume that each HRV type generates strong and lasting homotypic immunity, that leads to type-specific herd immunity and fade out. We have added text in the first paragraph of Discussion section.
Could the author add some referenced to pre-existing knowledge where discussion comments have about (a) of heterotypic immunity to RVs and (b) virus interactions?

Responses
Heterotypic immunity against HRV infections is generally weak and short-term (lasting between 5 weeks and 16 weeks) (Fleet et al. 1965(Fleet et al. , 1968Stott et al. 1969) compared to strong homotypic immunity (~1 year) (Barclay et al. 1989). Although heterotypic immunity can be established it is not significant to protect an individual from other multiple HRV types (Stott et al. 1969). The duration of homotypic immunity may be enhanced by reinfection with other types due to non-specific boosting of titers of types to which individuals were previously exposed (Fleet et al. 1965). Our own recent study of household HRV infection revealed rapid heterotypic reinfection in younger individuals (e.g. the median time to reinfection in 1-10 year olds was 1month), with a decreasing rate in older individuals ) who presumably had past experience of more types and hence reduced susceptibility.
We have also added references (Mackay 2008;Brunstein et al. 2008) in the Discussion section (page 14), which comments on virus interactions. Mackay's review of human rhinoviruses (Mackay 2008) indicates that virus interactions could be the cause of seasonal variation in the prevalence of any virus, whereby the high prevalence of one respiratory virus moderates or prevents processes that allow one or more other viruses to establish themselves in the host at the same time. The review also suggests that HRV partnership with host immunity can be mutual, incidentally giving the host an advantage by protecting against more cytopathic viral pathogens in the respiratory system while the host provides a vessel for HRV replication and transmission.
In their summary, the authors imply that the patterns in genotypes and species exchange detected in the nine health facilities are the same as those in the community. How sure can they be that there is not a mild and asymptomatic nonpresenting infection?

Response
The reviewer is right, we might have missed the mild and asymptomatic non-presenting infections since we only sampled from patients who presented to health centers with symptoms of ARI. However, we believe our surveillance is representative of circulation of all HRV infections because (a) we have a representative spatial distribution of health facilities and (b) our study of the household where we collected samples irrespective of symptoms showed similar type distributions for symptomatic and asymptomatic individuals  .
Could the authors expand on their 22.2% VP4/VP2 sequencing failure in the context of the comment about changed methods above

Response
We are not certain of the reasons for the 22.7% failure (240 samples) in sequencing. Viral load was slightly lower for failed samples: the mean Ct value for sequenced versus failed samples was 29.7 (standard deviation 3.1, 817 samples) and 30.7 (standard deviation 3.2, 240 samples) respectively. Furthermore, we had 154/240 samples failed at the PCR stage, which could have resulted from mismatches in the primer binding sites or assay sensitivity given that we were using a single step rather than nested PCR.