Genome Analysis of a Novel Tembusu Virus in Taiwan

We identified and isolated a novel Tembusu virus (TMUV) strain TP1906 (TMUV-TP1906) from a Culex annulus mosquito pool collected from the northern part of Taiwan in 2019. The TMUV-TP1906 genome is a 10,990-nucleotide-long, positive-sense, single-stranded RNA, consisting of a single open reading frame (ORF) encoding a polyprotein of 3425 amino acids, with 5′ and 3′ untranslated regions (UTRs) of 94 and 618 nucleotides, respectively. The nucleotide sequence of the TMUV-TP1906 of ORF exhibited 93.71% and 91.27% similarity with Sitiawan virus (STWV) and the TMUV prototype strain MM1775, respectively. The 3′-UTR variable region of TMUV-TP1906 showed nucleotide sequence divergence with other TMUV strains. Phylogenetic analysis of the complete ORF and polyprotein sequences revealed that TMUV-TP1906 is most closely related to STWV which causes encephalitis and retarded growth in chickens. We found that the TMUV-TP1906 caused a cytopathic effect (CPE) in the DF-1 chicken fibroblast cell line, while no apparent CPE was observed in Vero and C6/36 cells. In this study, we first identified and isolated a novel TMUV strain in Taiwan. In addition, to our knowledge, it is the first time that the TMUV strain was isolated from the Cx. annulus mosquitoes. Further study is warranted to investigate the host range and virulence of TMUV-TP1906.


Introduction
Tembusu virus (TMUV) is a single-stranded, positive-sense RNA virus that belongs to the Flaviviridae family, Flavivirus genus and Ntaya virus group. The TMUV prototype strain MM1775 (TMUV-MM1775) was first isolated from Culex tritaeniorhynchus mosquitoes in Kuala Lumpur, Malaysia in 1955. Since 2000, many TMUV strains (TMUVs) have been isolated from birds and mosquitoes, including chicks, ducks, geese, pigeons, sparrows and Culex mosquitoes. TMUV consists of several genetically closely related virus strains with noticeable pathogenicity for poultry, such as Sitiawan virus (STWV), duck TMUV (DTMUV), Perak virus and Baiyangdian virus [1]. STWV was isolated from sick broiler chicks in Sitiawan District of Perak State, Malaysia in 2000 [2]. STWV-infected chicks showed encephalitis, growth retardation and increased blood sugar levels. In 2010, DTMUV, also called duck egg-drop syndrome virus, was found to extensively infect ducks with high morbidity (up to 90%) and mortality (5% to 30%) rates in southeast China [3,4]. DTMUV outbreaks also occurred in layer and broiler duck farms in Pekin ducks in Malaysia in 2012 and in Thailand in 2013 [5,6].
TMUVs have been isolated from several mosquito species in Asia, mainly Cx. tritaeniorhynchus [8,9] and also Culex species, including Cx. vishnui, Cx. gelidus and Cx. pipiens [5,9,10]. A mosquito vector competence study by O'Guinn et al. demonstrated that Cx. vishnui developed high viral titers after feeding on TMUV-infected chicks and could readily transmit TMUV to naive chickens [11]. Although TMUV are considered mosquito-borne viruses [4], in vivo studies suggested that TMUVs could be transmitted among ducks by direct contact and aerosol transmission [12,13]. The role of mosquitoes as vectors in the transmission of TMUV to avian hosts still needs further study.
Taiwan is an island off the southeastern coast of mainland China that is near the epicenter of DTMUV outbreaks. Taiwan straddles the Tropic of Cancer, having a warm tropical-subtropical climate and large variety of mosquito species and is also an important rest point for migrating birds [14]. To monitor arboviruses in Taiwan, we conducted a survey on mosquitoes collected from wetlands, pig farms and parks in Taiwan between April and September 2019. In this study, we identified novel TMUV strains in Cx. annulus and Cx. tritaeniorhynchus mosquitoes collected from northern and central parts of Taiwan, respectively. We isolated and characterized the genetic sequence of this novel virus.

Collections of Mosquitoes
The mosquitoes were collected on pig farms near rice paddy fields in the northern (Wujie Township in Yilan County), central (Wufeng District in Taichung City), southern (Xiaying District in Tainan City), and eastern (Shoufeng Townships in Hualien County) regions; from a wetland habitat for waterbirds in the northern (Beitou District in Taipei City) region; and also from community parks in the northern (Shilin District in Taipei City and Luodong Township in Yilan County), central (Nantun District in Taichung City) and southern (North District in Tainan City) regions of Taiwan between April and September, 2019 ( Figure 1). The mosquitoes were collected using dry ice traps in wetland and parks and sweep nets in pig farm and were transported either alive or on dry ice to the laboratory. Mosquito collections by sweep nets were conducted only on the same day between 18:30 and 20:30 in the pig farms, while dry ice traps were set up overnight from 17:00 to 8:00 the next morning in wetland and parks as previously described [15]. Only female mosquitoes were analyzed in this study. The mosquitoes were frozen at -20 • C and further examined and classified according to the characteristics of the species. The pooled (1 to 50) mosquitoes were used for RNA extraction and virus isolation. The mosquito pools were homogenized in a TissueLyzer (Qiagen GmbH, Hilden, Germany) with two cycles at 4 • C for 90 sec at a frequency of 30 Hz after adding a 3 mm steel ball to each tube and 0.6 mL RPMI 1640 liquid medium. The pools were then clarified by centrifugation. The supernatants were sterilized by filtration and removed for RNA extraction and virus isolation.

RNA Extraction and Real Time RT-PCR
Viral RNA was extracted from the mosquito suspensions or cell culture supernatants using the QIAamp viral RNA mini kit (Qiagen GmbH, Hilden, Germany). The flavivirus-specific primers (1370F: 5 -TGY GTB TAC AAC ATG ATG GG; 1442F: 5 -ATA TGG TAC ATG TGG CTA GGA GC and 1620R: 5 -GTG TCC CAN CCH GCT GTG TCA) were mixed and used for the RT-PCR screening assay. Real-time RT-PCR was performed using a QuantiTect SYBR Green RT-PCR kit (Qiagen GmbH, Hilden, Germany) with the following parameters: 50 • C for 30 min; 95 • C for 15 min and 45 cycles of 95 • C for 15 s, 55 • C for 30 sec, 72 • C for 20 s, and 77 • C for 20 s; and 95 • C for 1 min, melting curve program from 68 • C to 95 • C [16,17]. DNA sequencing of positive RT-PCR products was performed, and RT-PCR positive mosquito pools were subjected to virus isolation.

Virus Isolation and Immunofluorescence Assay (IFA)
Filtered homogenates of mosquito pools were initially cultured with the Aedes albopictus C6/36 cell  line or Vero cell line. C6/36 cells were infected with TMUV-TP1906 and tested by IFA. Briefly, cells  were fixed and incubated with anti-flavivirus monoclonal antibody (D56.3) for 1 h. The cells were subsequently washed and incubated with FITC-conjugated goat anti-mouse IgG antibody (Invitrogen, Carlsbab, CA, USA) for 1 h. After washing, the images of the cells were examined using a fluorescence microscope (Carl Zeiss, Oberkochen, Germany).

Complete Genome Sequencing
RNA extracted from the cell culture supernatants and pellets from the first cell culture passage were used for RT-PCR and DNA sequencing. RT-PCR was performed using the Superscript III One-Step RT-PCR system with Platinum Taq High Fidelity (Invitrogen, Carlsbab, CA, USA). The RT-PCR parameters were as follows: 55 • C for 30 min; 94 • C for 2 min and 40 cycles of 94 • C for 15 s, 50 • C or 54 • C for 30 s, and 68 • C for 1 min; and a prolonged final extension at 68 • C for 5 min. The cDNA of the 5 end of the viral genome was generated by rapid amplification of cDNA ends (RACE) using a 5 /3 RACE kit, 2nd generation (Roche, Mannheim, Germany) with a gene specific primer (GSP-5 end, Table 1). For the 3 end, polyadenylated (polyA) tails were added to the 3 ends of virus genomic ssRNA using T4 RNA ligase. Tailing reactions were performed at 37 • C for 1 h and used for RT-PCR directly. The 5 and 3 end fragments were amplified using the genome-specific primer GSP-5 end and 3 UTR-10 with oligo d(T)-anchor primer (provided by the 5 /3 RACE kit), respectively (listed in Table 1). The PCR products were sequenced directly by using the BigDye Terminator Cycle Sequencing Kit and the ABI 3730xl DNA analyzer (Applied Biosystems) according to the manufacturer's protocols. Each forward and reverse primer was used for sequencing, and the overlapping sequences were combined. The sequence identity of the ORF and protein sequence were analyzed by BLAST (https://blast.ncbi.nlm.nih.gov/Blast.cgi). The sequence identities of the 3 -UTR and 3 -UTR variable region were analyzed by Clustal Omega (https://www.ebi.ac.uk/Tools/msa/clustalo/).

Phylogenetic Analysis
Genome sequences of TMUVs were aligned, edited, and analyzed using Clustal W software [20]. The complete open reading frame sequences of TMUVs including sequences representing the most closely related strains to the TMUV strain TP1906 isolated in this study were obtained using BLAST and retrieved from GenBank for phylogenetic analysis ( Table 2). The phylogenetic analysis was performed using MEGA version 7 (http://www.megasoftware.net/) [21]. To construct the phylogenetic trees, the maximum likelihood method using the Tamura-Nei model (Figures 3 and 5) and the maximum likelihood method using the JTT matrix-based model ( Figure 4) were utilized. The reliability of the analysis was evaluated by a bootstrap test with 1000 replications. Sequences of Bagaza virus (GenBank numbers: KR108245 and MF380425), Japanese encephalitis virus (GenBank number: KF667310) and West Nile virus (GenBank number: M12294) were used as the outgroup.

Identification of Novel TMUV Strains from Culex Mosquitoes in Taiwan
Mosquitoes were collected from nine localities in Taiwan from April to September 2019 ( Figure 1). A total of 15,374 mosquitoes were collected and analyzed. Table 3 shows a summary of the mosquito species, numbers and RT-PCR results. Nineteen mosquito species from nine genera of the Culicidae family were identified. The most commonly collected species was Cx. tritaeniorhychus (70.90%, n = 10,900), followed by Cx. annulus (11.87%, n = 1825) and Mansonia uniformis (5.45%, n = 838). A flavivirus-specific primer set targeting the NS5 gene (amplicon is approximately 200 bp) was initially used for screening the mosquito pools by real-time RT-PCR. A total of 429 mosquito pools were tested, and 55 pools were positive; RT-PCR products were then sequenced and analyzed using the NCBI BLAST platform. Of the 55 positive pools, 53 pools contained genomic Japanese encephalitis virus (JEV) sequences, while a Cx. annulus mosquito pool (TP1906) collected from a wetland habitat for waterbirds in the Beitou District, Taipei City, in northern Taiwan, and a Cx. tritaeniorhynchus mosquito pool (TC1906) collected from a pig farm near rice paddy fields in Wufeng Township, Taichung City in central Taiwan, contained TMUV-like sequences. The partial NS5 gene sequences of RT-PCR products detected in the TP1906 and TC1906 mosquito pools were 100% identical and are closely related to STWV (GenBank number: JX477686) and TMUV-MM1775 (GenBank number: MH414569) sequences with 92.9% and 92.31% nucleotide similarities, respectively. To further study whether the novel flaviviruses detected in Cx. annulus and Cx. tritaeniorhychus mosquito pools were identical, we amplified a partial NS1 gene sequence (964 bp) by using a TMUV NS1-specific primer set (TMUV-NS1_F and TMUV-NS1_R, in Table 1) [19]. The amplified NS1 gene sequences from these two mosquito pools were 99.79% identical with two nucleotide differences. The sequence obtained from the Cx. tritaeniorhynchus mosquito pool (TC1906) had a serine at amino acid position 131 of the NS1 protein (GenBank number: MN958524), while the Cx. annulus mosquito pool (TP1906) had an asparagine at position 131 (GenBank number: MN747003). The results indicated that the viruses from these two mosquito pools were very similar.

Isolation of a Novel TMUV from the Cx. Annulus Mosquito Pool
The C6/36 and Vero cell lines were used for propagation of the novel TMUV from the TP1906 and TC1906 mosquito pools. The novel TMUV strain TP1906 was successfully isolated from the Cx. annulus mosquito pool (TP1906). Virus isolation from the Cx. tritaeniorhynchus mosquito pool (TC1906) was unsuccessful.

Different Cell Lines Were Permissive for TMUV-TP1906 Infection
Immunofluorescence assay using anti-flavivirus antibody indicated that the TMUV-TP1906 was capable of infecting C6/36 cells (Figure 2a). To determine whether the TMUV-TP1906 replication can take place in different cell lines derived from mosquito (C6/36), monkey (Vero), hamster (BHK-21) and chicken (DF-1), we inoculated the TMUV-TP1906 to these cell lines and observed the cytopathic effects and detected viral RNA by using real-time RT-PCR. Although lacking a noticeable cytopathic effect (CPE) in C6/36 and Vero cell cultures, we observed apparent CPE in DF-1 and BHK-21 cell lines 72 h postinfection (Figure 2b). Notably, all the infected DF-1 and BHK-21 cells were detached five days postinfection. The TMUV-TP1906 was capable of replication in DF-1 cells and reached its highest titer at 48 h postinfection (Figure 2c). In contrast, the virus replication was relatively slowly in Vero cells compared to the BHK-21 and C6/36 cells.

Phylogenetic Analysis of TMUV-TP1906 and Other Flaviviruses Based on ORF Sequences
To understand the phylogenetic relationship between TMUV-TP1906 and other mosquito-borne flaviviruses, the complete ORF sequences of 40 representative flaviviruses, namely 36 TMUVs (listed in Table 2), Bagaza virus (GenBank numbers: KR108245 and MF380425), Japanese encephalitis virus (GenBank number: KF667310) and West Nile virus (GenBank number: M12294), were analyzed. Phylogenetic analysis indicated that TMUVs could be divided into two lineages (lineages TMUV and DTMUV). The TMUV lineage contained TMUV-MM1775, STWV, and TMUV-TP1906, and TMUV-TP1906 was most closely related to STWV. The DTMUV lineage could be further grouped into three clusters (clusters 1 to 3) (Figure 3). Cluster 1 contained DTMUV strains from Malaysia and Thailand. Clusters 2.1 and 2.2 contained many avian DTMUV strains found in Thailand and China. DTMUV-SD14 isolated from mallards (Anas platyrhynchos) in China was classified into cluster 3.  (10,278 nt). The evolutionary history was inferred by using the maximum likelihood method based on the Tamura-Nei model [22]. The tree with the highest log likelihood is shown. The percentage of trees in which the associated taxa clustered together is shown next to the branches. Initial tree(s) for the heuristic search were obtained automatically by applying Neighbor-Join and BioNJ algorithms to a matrix of pairwise distances estimated using the maximum composite likelihood (MCL) approach, and then selecting the topology with superior log likelihood value. The tree is drawn to scale, with branch lengths measured in the number of substitutions per site. The analysis involved 40 nucleotide sequences. All positions containing gaps and missing data were eliminated. Evolutionary analyses were conducted in MEGA7 [21]. The reliability of the analysis was calculated using 1000 bootstrap replication. Bootstrap support values >75 are shown. The scale bar indicates nucleotide substitutions per site.

Amino Acid Sequence Comparison between TMUV-TP1906 and Other TMUVs
Genomic sequences of TMUV-TP1906 and other TMUV strains contain a single ORF that encodes a polypeptide of 3425 amino acids, except the ORFs of two DTMUV strains (D1977 and D1921) that encode 3424 amino acids with a deletion in the NS5 protein at amino acid position 86. The ORFs of TMUV-TP1906 and other TMUVs shared 86.52%-93.71% similarity at the nucleotide level and 95.71%-98.92% similarity when the amino acid sequences of the polyproteins were compared ( Table 4). The amino acid sequences of the prM, E, and NS1 proteins of TMUV-TP1906 showed the highest similarity with TMUV-MM1775, while sequences of the C, NS2A, NS3, NS4A, NS4B, and NS5 proteins of TMUV-TP1906 showed the highest similarity with STWV. Figure 4 shows the phylogenetic relationship based on the amino acid sequences of polyproteins of TMUV-TP1906 and TMUVs. The results indicated that TMUV-TP1906 is most closely related to STWV.

Sequence Alignment and Phylogenetic Analysis of 3 -UTR Variable Region of TMUV-TP1906 and Other TMUVs
Comparison of 5 -UTR and 3 -UTR sequences between TMUV-TP1906 and other TMUVs showed 90.43%-94.68% and 89.28%-95.15% sequence similarities, respectively ( Table 4). Analysis of 3 -UTR variable region (VR) sequences (lengths varied from 84 to 94 nucleotides) between TMUV-TP1906 and other TMUVs showed that TMUV-TP1906 contained 1 gap (gap 1) and DTMUVs contained 2 gaps (gaps 1 and 2), while no gaps were found in the VR sequences of TMUV-MM1775 and DTMUV-SD14 (Figure 5a). The novel TMUV-TP1906 had a 3 nucleotide insertion in gap 2 and a nucleotide insertion in gap 1. In addition, there were nine unique nucleotide substitutions in the VR of TMUV-TP1906 compared with those of the other TMUVs. The 3 -UTR VR sequences were more consistent among cluster 2 of the DTMUV lineage, but more variable in the TMUV lineage and clusters 1 and 3 of the DTMUV lineage (Figure 5b). Comparisons of the 3 -UTR VR sequence of TMUV-TP1906 with those of TMUV-MM1775, DTMUV-SD14 (cluster 3), DTMUV-D1921 (cluster 1), and DTMUV-D1977 (cluster 1) showed 74.71%, 71.26%, 65.00%, and 66.25% sequence similarities, respectively, and 56.96%-60.76% sequence similarity with cluster 2 of DTMUVs. Figure 5b shows the phylogenetic tree based on the 3 -UTR VR sequences of TMUVs, and the overall topology was similar to the phylogenetic tree based on analysis of ORF sequences (Figure 3). . The evolutionary history was inferred by using the maximum likelihood method based on the JTT matrix-based model [23]. The tree with the highest log likelihood is shown. Initial tree(s) for the heuristic search were obtained automatically by applying Neighbor-Join and BioNJ algorithms to a matrix of pairwise distances estimated using a JTT model, and then selecting the topology with superior log likelihood value. The tree is drawn to scale, with branch lengths measured in the number of substitutions per site. The analysis involved 40 amino acid sequences. All positions containing gaps and missing data were eliminated. Evolutionary analyses were conducted in MEGA7 [21]. The reliability of the analysis was calculated using 1000 bootstrap replication. Bootstrap support values >75 are shown. The scale bar indicates amino acid substitutions per site. The unique amino acid residues in TP1906, Sitiawan virus, MM1775, SD14, and DTMUV cluster 1 compared to other TMUV strains were colored in blue, brown, gray, green, and yellow, respectively. The unique amino acid residues in DTMUV cluster 1 and 2, TMUV lineage, and TP1906 and STWV compared to other TMUV strains were colored in black, purple, and red, respectively.

Discussion
In this study, we identified two novel TMUV strains, TMUV-TP1906 and TMUV-TC1906, from Cx. annulus and Cx. tritaeniorhynchus mosquitoes collected from a wetland in northern Taiwan and a pig farm in central Taiwan, respectively. TMUV-TP1906 was isolated from the Cx. annulus mosquito pool and grew well in Vero and C6/36 cells without significant cytopathic effects. Through virus isolation and comprehensive genomic and protein sequence analysis, we found that this novel flavivirus is a TMUV-related virus and is most closely related to STWV. STWV, DTMUV, and Bagaza virus, which belong to Ntaya serocomplex virus; and JEV and West Nile virus, which belong to JEV serocomplex virus, have been shown to be primarily Culex mosquito-associated viruses that cause severe diseases in avian species [7]. Whether TMUV-TP1906 is a pathogen in birds or other animals needs further study.
Phylogenetic analysis based on the complete ORF of TMUVs showed that TMUVs can be divided into two lineages, TMUV and DTMUV. Analysis of amino acid sequences showed that there are 22 amino acid substitutions unique to the TMUV lineage and 29 amino acid substitutions unique to cluster 1 and cluster 2 of the DTMUV lineage ( Figure 4). The DTMUV-SD14 strain which belongs to cluster 3 of the DTMUV lineage, could be considered as a transitional group between the TMUV lineage and clusters 1 and 2 of the DTMUV lineage. The DTMUV-SD14 strain contains 48 amino acid mutations, which are different from other TMUVs. In addition, there are five informative amino acid substitutions (NS1-105, NS1-205, 2K-13, NS4B-30 and NS5-562) for the differentiation of the TMUV lineage, DTMUV-SD14 (cluster 3 of the DTMUV lineage) and clusters 1 and 2 of the DTMUV lineage. The TMUV lineage contains three virus strains, the TMUV-MM1775 prototype strain, STWV, and TMUV-TP1906. TMUV-MM1775 has not been documented as a pathogen in birds or other animals; however, STWV, which is most closely related to TMUV-TP1906, causes encephalitis and retards the growth of chicks. Analysis of amino acid sequences showed that TMUV-TP1906 contains 12 unique amino acid substitutions compared with other TMUVs. In addition, comparison of protein sequences of TMUV-TP1906 and STWV with other TMUVs showed 10 amino acid substitutions. These informative amino acid positions may be useful for the classification of TMUVs and may play crucial roles in the evolution and protein function of TMUVs.
The Asn glycosylation site at residue 154 of the E protein (154-Asn) of TMUVs is critical for virus tissue tropism and transmissibility in poultry [24]. TMUV-MM1775 contains 156-Pro, which can disrupt the N-linked glycosylation of 154-Asn of the E protein, and thus may affect viral virulence and transmissibility in avian species [24]. Except for TMUV-MM1775, all the other TMUVs, including TMUV-TP1906, contain 156-Ser which would not disrupt 154-Asn glycosylation and thus may not affect the glycosylation-associated infectivity of TMUVs in avian species. A recent study showed that the Thr to Lys mutation at residue 367 of the E protein plays a predominant role in viral cell-adaptation and virulence attenuation in ducks [25]. TMUV-TP1906 contains 367-Thr of the E protein and thus may retain the virulence in the host. In addition to E protein, NS1 protein contains three N-linked glycosylation sites at residues 130, 175, and 207. A recent study showed that some DTMUV strains had a glycosylation motif mutation at residues 175-178 from NTTD to NITD but the glycosylation site at 130 was conserved [19]. Interestingly, we found NS1 glycosylation motif at residues 175 and 207 were conserved in TMUV-TP1906 and TMUV-TC1906, however, TMUV-TC1906 had a glycosylation motif mutation at residues 130-133 from NNTF to NSTF. Whether these mutations affect protein functions or virulence needs further study.
The 5 and 3 UTRs of flaviviruses are important for virus replication and transmission [26]. The 5 end of the 3 UTR of flavivirus is known to have a VR with different nucleotide lengths [26], and this VR sequence may be associated with the production of subgenomic flavivirus RNA and immune modulation, hence contributing to virus adaptation in vector and non-vector hosts [27,28]. The 3 UTR VRs of TMUVs consist of 84 to 94 nucleotides immediately downstream of the ORF. Phylogenetic analysis of 3 UTR VR and ORF sequences showed similar topological features, indicating that 3 UTR VRs may serve as a potential marker for TMUV evolution. Sequence analysis also showed that the 3 UTR VR of TMUV-TP1906 is very unique among TMUVs ( Figure 5). More research is needed to understand the genetic diversity and function of the 3 UTR VR of TMUVs.
Prominent CPE was observed in BHK-21, C6/36, Vero and DF-1 cell lines after infection with DTMUV [5,10,18,29]. STWV infection induced CPE in the BK3 cells (chicken bursal lymphoma cell line) but not in CPK (porcine kidney cell line), MARC-145 (monkey kidney cell line) and Vero cells [2]. In our study, we found that the TMUV-TP1906, which is closely related to STWV, grew in C6/36 and Vero cell without inducing significant CPE, however, the virus caused drastic CPE in DF-1 and BHK-21 cell lines (Figure 2).
TMUVs are emerging pathogenic flaviviruses causing severe avian diseases in Malaysia, Thailand, and China. Taiwan is located south of East Asia and is close to the epicenter of DTMUV outbreaks in China. In this study, we first identified a TMUV in Taiwan. Interestingly, TMUV-TP1906 is most closely related to STWV and TMUV-MM1775 found in Malaysia and is different from the DTMUVs found in China (Figures 3 and 4 Figure 3; Figure 4). Previous study has divided TMUVs into the Chinese mainland TMUV lineage and Southeast Asian TMUV lineage [8]. According to the phylogeny, TMUV-MM1775, STWV, TMUV-TP1906, and two TMUV strains identified from Cx. tritaeniorhynchus in Yunnan Province of China were grouped as the Southeast Asian TMUV lineage. The phylogeographical analysis suggested that the TMUVs might spread from Southeast Asian countries such as Malaysia and Thailand to the Yunnan Province of China (southern China) and further spread to areas in northern China such as Shandong Province. However, the mechanism of long-distance spread of TMUVs is still unknown. Liu et al. suggested that the spread of TMUVs might occur via migratory birds and ornithophilic Culex spp. mosquitoes, since TMUVs have been identified in wild birds such as pigeons and mallards (Anas platyrhynchos) [7]. The East Asian-Australasian flyway is one of the world's great flyways for migratory birds, and the flyway passes through many countries including Malaysia, Thailand, China, and Taiwan; thus, it is possible that the spread of TMUVs might be through this flyway. In addition, it has been suggested that the spread of avian influenza virus and JEV might also occur through this flyway [30,31].
In this study, we reported the first isolation of a novel TMUV strain from Culex mosquitoes in Taiwan. In addition, it is the first time that the TMUV strain was isolated from Cx. annulus mosquitoes. We performed comprehensive genomic and amino acid analyses between TMUV-TP1906 and other TMUVs, and our results may help to increase understanding of the diversity and evolution of TMUVs. Further study is warranted to investigate the virulence and epidemiology of TMUV-TP1906 in Taiwan.