A five-year retrospective study on the epidemiology of hand, foot and mouth disease in Sabah, Malaysia

Hand, foot, and mouth disease (HFMD) is endemic in Malaysia, with the number of cases increasing. Sabah has experienced several HFMD outbreaks, but information on the epidemiology and molecular characteristics of responsible viruses is scarce. In this study, data of 17,574 reports of HFMD cases in Sabah from 2015 to 2019 were extracted from a public health disease surveillance system and analyzed. Twenty-one swab samples from 13 children were collected from Beaufort, Sabah, during an outbreak in August 2018 for detection and serotyping of causative viruses by semi-nested reverse transcription-polymerase chain reaction (snRT-PCR) of the VP4–VP2 region and consensus degenerate hybrid oligonucleotide primer PCR of the VP1 region, respectively. Nucleotide sequencing and phylogenetic analysis were conducted by the neighbor-joining method. The average annual incidence of HFMD was 94.3 per 100,000 people, with the greatest yearly increase between 2017 and 2018. Swabs from six children were tested positive for enterovirus, of which five were positive for CVA16 and one for EV71. All CVA16 strains belonged to sub-genotype B1a, and the EV71 strain belonged to sub-genotype B5. Phylogenetic analyses indicate that enterovirus genotype shift might be responsible for the increasing trend of HFMD in Sabah, however, further study is needed.

www.nature.com/scientificreports/ Several researchers have characterized EV71 and CVA16 strains from HFMD cases in Peninsular Malaysia and Sarawak [10][11][12][13] . Although a number of HFMD outbreaks have been reported in Sabah, information on the genetic characteristics of these strains in Sabah is limited. So far, only two studies have been conducted on EV71 strains isolated from Sabah in 1999 11 and 2008 12 , but none on CVA16. Since then, no studies have been conducted on HFMD in Sabah for almost a decade. The limited research could be due to the assumption of low burden of HFMD in Sabah and a lack of standardized sample collection protocols that would maximize the yield of viruses for molecular characterization. More surveillance is necessary to understand and monitor the emergence and spread of new strains. Besides, understanding the epidemiological pattern of HFMD is important for the implementation of appropriate intervention strategies. Therefore, this five-year retrospective study explored the epidemiological pattern of HFMD in Sabah using the data of 17,574 reports of HFMD cases recorded in a Malaysian public health disease surveillance system, eNotifikasi, and investigated the genotype characteristics of seven enterovirus-positive samples that caused HFMD in the region.

Results
Temporal and demographic trends. From the eNotifikasi system (http:// enoti fikasi. moh. gov. my), a total of 17,574 reports of HFMD cases were recorded in Sabah during the five-year study period, with no deaths reported. While cases were notified year-round, seasonal analysis showed that reported cases peaked annually between January and March, except in 2018, where HFMD peaked in July (Fig. 1). The number of cases gradually increased from 2015 to 2019, but there was a slight decrease in 2017 from the previous year (Table 1). The annual incidence rate from 2015 to 2019 varied from 39.9 per 100,000 people to 166.1 per 100,000 people, with a   (Table 2). On the contrary, Tongod district had the lowest number of cases, with just 19 cases reported throughout the five years. Of the 25 districts, a peak in the number of cases was observed in 17 (68%) districts in 2019, five (20%) districts in 2018, two (8%) districts in 2016, and one (4%) district in 2015. The annual spatial trend dynamics are shown in Fig. 2, indicating that there was a yearly variation across districts. Throughout the five-year study period, the district with the largest variation in annual incidence rate was Beaufort (58 per 100,000 people in 2016 to 657 per 100,000 people in 2018) while Kinabatangan had the smallest variation in annual incidence rate, ranging from 8 per 100,000 people in 2015 to 18 per 100,000 in 2018. A prominent spatial change was observed in 2018, when incidence increased markedly in the western regions of Sabah from the previous year. High incidence rates were not only found in densely populated districts, such as Kota Kinabalu, Putatan and Penampang, but also in low population density districts, including Sipitang, Beaufort, Kuala Penyu and Kota Belud. The highest incidence district was Kuala Penyu with a five-year average annual incidence of 257 per 100,000 people, while the lowest were Tongod and Kunak districts, both with 9 per 100,000 people.
Phylogenetic analysis of enteroviruses. A total of 21 swab samples were collected from 13 children during an outbreak reported in Beaufort, Sabah, in August 2018. Swabs from six children comprising of four females and two males, ranging from one to seven years old were tested positive for enteroviruses using snRT-PCR. Sequence analysis of the VP1 gene of the positive snRT-PCR samples by BLAST identified five (83.3%) CVA16 and one (16.7%) EV71 serotypes. The partial VP1 sequences of the CVA16 and EV71 viruses were 254 bp and 359 bp in length, respectively. These partial VP1 sequences were used to construct the phylogenetic trees of the two serotypes. This is the first study that shows the genetic characteristics of CVA16 strains from Sabah. According to the phylogenetic analysis, all five CVA16 strains detected in this study were classified as genotype B, specifically subgenotype B1a. Three CVA16 strains clustered with Thai strains from 2017, while the other three clustered with The only EV71 strain detected in this study was classified as a genotype B, precisely sub-genotype B5 (Fig. 4). The strain was clustered closely together with five EV71 strains from Thailand detected in 2012, 2014, and 2017, and a strain from France detected in 2013. However, it did not cluster with any previously detected Malaysian strains. Strains isolated from Sarawak and Sabah in 2003 and 2008, respectively, also belonged to the same subgenotype but of different lineage. Sub-genotype B4 was found circulating in Sabah and Sarawak in 1999 and 2000, respectively, while B3 was the predominant sub-genotype during the 1997 outbreak in Sarawak. From this phylogenetic analysis, genotypic changes of EV71 have been observed over time even in the same region, for instance Sarawak.

Discussion
There has been a considerable increase in reports of HFMD cases in Sabah throughout the five-year study period. HFMD cases peaked in the early months, between January and March every year, except in 2018, where the highest number of reported cases was in July. The exact causes of the increment of cases and the notable seasonal characteristic of HFMD as well as the different peak period in 2018 are unclear. However, there are several factors that could possibly affect HFMD incidence.
A study by Zhao and Hu (2019) found that school terms and a major public holiday in China, the Chinese Spring Festival, affected the HFMD transmission seasonality in mainland China, contributing to two major waves in the region 14 . The study revealed that the Chinese Spring Festival reflected the seasonal contact rate in the population, which dominantly caused the larger peak in March, while the school terms reflected the seasonal contact rate in children, resulting in the smaller peak in autumn 14 . The school year in Malaysia starts in early January and it is divided into two school terms: from January to May for the first term and from June to November for the second term. Besides, Chinese New Year, which is the same as the Chinese Spring Festival, is one of Malaysia's biggest holidays celebrated in the beginning of the year, between late January and mid-February. The first day of school and the Chinese New Year celebration in Malaysia when combined, could be a driving factor for the seasonal peaks (between January and March every year, except for 2018) revealed in this study. However, further in-depth research is necessary to determine the effects of school terms and holidays on HFMD transmission in Malaysia.
In 2018, a large HFMD outbreak was reported in Malaysia, which affected over 76,000 children nationwide 15 . Selangor recorded the highest number of HFMD cases on August 4, 2018 with 12,868 cases, followed by Kuala Lumpur with 4996 cases and Sarawak with 4988 cases 16 . Although the number of HFMD cases in Sabah was not as high as these states, the outbreak undoubtedly had a significant impact on the local community as the region in general, experienced a sudden spike in HFMD cases in 2018. As shown in this study, cases in 2018 continued to soar upwards until July of the same year, unlike the previous years and 2019, where cases started to decrease after March, which could be due to exhaustion of susceptible hosts 14 or effect of control measures. A possible cause for the unusual pattern in 2018 might be due to the movement of people traveling to their respective locations for polling during Malaysia's 14th General Election on May 9, 2018. This was later followed by two major celebrations, the Harvest Festival, which is an important regional event and Eid-ul-Fitr (Festival of Breaking the Fast) celebrated in late May and mid-June, respectively.
Meteorological factors, such as air temperature, rainfall, relative humidity, air pressure and wind speed, have been shown to impact HFMD differently depending on the climatic zones and latitudes 17,18 . A meta-analysis by Duan et al. (2019), revealed that HFMD incidences in subtropical and temperate regions were significantly associated with meteorological factors but not in tropical regions 17 . On the contrary, other studies of tropical www.nature.com/scientificreports/ climate found significant correlation between HFMD and weather variables, including air temperature, rainfall, pressure, and wind speed 19,20 . These show that meteorological factors have inconsistent and varying effects even in regions with the same climate, However, the effects of meteorological factors on HFMD incidence in this study were not examined and thus, the question on whether these factors contribute to the seasonal peaks in Sabah remains an enigma. www.nature.com/scientificreports/ In this study, CVA16 and EV71 were detected from the swab samples collected during the HFMD outbreak in Beaufort in August 2018, where the former was the predominant serotype. Other recent studies of HFMD in Peninsular Malaysia also reported the presence of CVA16 and EV71 in their samples, besides CVA6, which was not detected in this study. Ling et al. (2014) identified CVA6, CVA16 and EV71 from specimens collected from patients in Seri Kembangan, Selangor, where CVA6 was the most common among the three types of enterovirus 10 . Another study by Lee et al. (2021) reported that only CVA6 and CVA16 were found in clinically diagnosed HFMD patients in Kuala Lumpur and Selangor, where CVA6-infected cases were more common in children under 12 months old 15 .
EV71 and CVA16 are regarded as the two major etiological agents of HFMD in Malaysia 9 but only in recent years has CVA6 become one of the predominant serotypes in the country, as reported in the other studies mentioned previously. In fact, the emergence of CVA6 as the pathogen responsible for HFMD outbreaks has been observed long ago in other countries, dating back to 2008 21 , but not in Malaysia. This might be due to the limited number of studies on the characterization and identification of enterovirus types that caused HFMD in Malaysia, especially in Sabah.
Four sub-genotypes, B3, B4, C1 and C2 of the neurovirulent EV71 were found co-circulating and caused the 1997 HFMD outbreak in Peninsular Malaysia, while sub-genotypes B3 and C1 co-circulated in Sarawak in the same year 9 . Sub-genotype B3 was the predominant type that was associated with fatal cases in 1997 7 . Later in 2000, sub-genotypes B4, B5 and C1 co-circulated during the outbreaks in Peninsular Malaysia and Sarawak, where B4 emerged as the predominant sub-genotype 7,9 . From 2003 onwards, EV71 in Malaysia was mainly of subgenotype B5, with minor contribution from sub-genotype C1 until 2005 [7][8][9] . The predominance of sub-genotype B5 in Malaysia remains even after a decade, as shown in Ling et al. (2014) 10 . The EV71 strain in the present study was also classified as sub-genotype B5 but with the limited number of samples in this study, further investigations with a larger sample size are necessary to determine the prevalence of sub-genotype B5 in Malaysia.
Genotype  11 . According to the sequence analysis of VP1 gene in this study, all CVA16 strains were classified as sub-genotype B1a, which has never been reported from Malaysia and is distinct from those isolated from Sarawak and Peninsular Malaysia. Phylogenetically, our CVA16 strains were close to the strains from Thailand, suggesting that they share a common ancestor. The first characterization of CVA16 and the identification of sub-genotype B1a in Sabah can shed some light on the genomics of HFMD-causing enterovirus in the region.
As in our study, EV71 and CVA16 have been identified circulating in several countries, including China, Japan, and Thailand 6 . Further, co-circulation of both CVA16 and EV71 serotypes during outbreaks has been reported in various countries, including Malaysia, Taiwan, and China 22 , which is a significant public health concern. It has been suggested that co-circulation of multiple serotypes of enterovirus may cause genetic recombination and co-infection, leading to increased risk of severe disease 23 .
In support of this notion, the present study found that HFMD mainly affected children under-five years of age, especially one-year-old children. This may be attributed to a decrease in maternal antibodies when they start weaning and more interaction with other children. A similar trend has been reported in other countries showing that HFMD cases were mostly found in children under the school-age 24 . However, infants less than one-yearold had a lower infection rate, presumably due to the lack of contact with other children or the protective effect of maternal antibodies from breast milk 25 . In agreement with other studies 24,26 , this study also revealed male predominance in HFMD cases, possibly attributed to the different degree of attention and treatment given by caregivers between boys and girls, especially in Asian countries. Caregivers are more likely to seek medical care for boys than girls 24 . On the contrary, in Singapore and Taiwan, gender did not significantly affect the infection rate 27,28 .
As shown in this study, other countries in Asia, such as China 29 , Singapore, and Japan 7 , also observed an increment of more than two-fold in the number of cases throughout the last decade. Besides, this study also found that areas with high population density and more developed economies, such as Kota Kinabalu and Penampang had high number of HFMD cases, which is consistent with previous studies 30 . This may be attributed to the economic development and urbanization in these areas leading to more frequent communication and close contact with individuals, resulting in a higher risk of HFMD transmission 31 . However, high incidence rates were also found in low populated rural areas, such as Sipitang and Kuala Penyu, which may be due to inadequate sanitary conditions, lack of clean water supply, poor hygiene practices and limited knowledge on the disease 32 .
In conclusion, HFMD in Sabah is often overlooked, and the disease warrants attention since the number of cases has increased remarkably over the last five years. Knowledge of the epidemiological and molecular characteristics of HFMD in Sabah is essential to assist in developing appropriate prevention and effective treatment strategies. This will contribute to the overall success of reducing the burden of communicable diseases in Malaysia. This study, for the first time, described the epidemiological features of HFMD in Sabah and investigated the enterovirus serotypes that might be responsible for the disease in the area. Due to the limited number of samples in this study, further studies are necessary to include more HFMD samples from different areas of Sabah for molecular analysis to gain a better understanding of the etiologic causes of HFMD in the region.  Enterovirus detection using semi-nested reverse-transcription polymerase chain reaction (snRT-PCR). Enteroviral RNA was detected in samples by snRT-PCR of the VP4-VP2 region. The first primer pair, MD91-OL68-1 (Table 3), was used to amplify the partial 5' non-coding region, VP2, and complete VP4. 3 µL of RNA was used in the PCR using the AccessQuick RT-PCR System (Promega, USA). The PCR condition consisted of 1 cycle of initial denaturation at 95 °C for 2 min, 40 cycles of denaturation at 95 °C for 30 s, annealing at 55 °C for 30 s, and extension at 72 °C for 1 min, followed by 1 cycle of final extension at 72 °C for 7 min. snRT-PCR with the second primer pair, EVP4-OL68-1 (Table 3), was performed by using 2 µL of the PCR product under the following PCR condition: 1 cycle of initial denaturation at 95 °C for 6 min, 40 cycles of denaturation at 95 °C for 30 s, annealing at 55 °C for 30 s, and extension at 72 °C for 1 min, followed by 1 cycle of final extension at 72 °C for 7 min.
Enterovirus serotyping using consensus degenerate hybrid oligonucleotide primer (CODE-HOP). The enterovirus VP1 gene sequence was amplified using a CODEHOP protocol described by Nix, Oberste, and Pallansch (2006) 35 . Briefly, the cDNA generated by the four different primers (AN32-AN35) was used for the first round of PCR with primers 224 and 222 (Table 3). The condition for the first PCR consisted of 1 cycle of initial denaturation at 95 °C for 6 min, 40 cycles of denaturation at 95 °C for 30 s, annealing at 42 °C for 30 s, and extension at 60 °C for 1 min. The first PCR generated a product of 992 bp, which was used for the second round of PCR with primers AN89 and AN88 (Table 3) for nested amplification. The condition for the second PCR consisted of 1 cycle of initial denaturation at 95 °C for 6 min, 40 cycles of denaturation at 95 °C for 30 s, annealing at 60 °C for 30 s, and extension at 72 °C for 30 s, followed by 1 cycle of final extension at 72 °C for 5 min. The second PCR generated a product of 375 bp. Sequence analysis. Partial nucleotide sequences of the VP1 gene of CVA16 (254 bp) and EV71 (359 bp) were used for phylogenetic analyses. The amplified DNA of each positive sample was sequenced using the Big-Dye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems, Foster City, CA) according to the manufacturer's instructions, and the product was run on an ABI Prism 3100 Genetic Analyzer (Applied Biosystems). Nucleotide sequences of VP1 were checked against the NCBI database by BLAST to find the enterovirus serotype with the highest identity. The nucleotide sequences of other CVA16 and EV71 strains were extracted from Gen-Bank. Multiple sequence alignment was performed using Clustal W, and phylogenetic trees were constructed using MEGA 6.0, applying the neighbor-joining method based on the Tamura-Nei substitution model 36 . Bootstrap analysis of 1000 replicates was conducted to determine the significance of the branching of the constructed tree. Nucleotide sequences analyzed in this study have been submitted to the DNA Data Bank of Japan (DDBJ). Data analysis. Data analysis was performed using Microsoft Excel for Mac 2017 (Version 15.30). QGIS v3.14.1 (Pi) was used to generate spatial maps of HFMD cases in Sabah.  www.nature.com/scientificreports/