Molecular epidemiology of Akabane virus in Taiwan.

BACKGROUND: Akabane virus (AKAV) is a teratogenic and neuropathogenic arbovirus that infects livestock and wild animals. AKAVs are endemic arboviruses from dairy farms in Taiwan in 1989, and the first sequence was detected in cattle with nonsuppurative encephalitis in 1992.
OBJECTIVES: This study aims to understand the epidemiological relationships of the akabane viruses between Taiwan and nearby places.
METHODS: In this study, 17 specimens were identified or isolated from vector insects, and ruminant fetuses collected from 1992 to 2015 were sequenced and analysed.
RESULTS: Sequence analyses revealed all Taiwanese AKAVs belonged to genogroup Ia but diverged into two clusters in the phylogenetic trees, implying that at least two invasive events of AKAV may have occurred in Taiwan.
CONCLUSIONS: The two clusters of AKAVs could still be identified in Taiwan in 2015, and a reassortment event was observed, indicating that the two clusters of AKAVs are already endemic in Taiwan.

INTRODUCTION

Akabane virus (AKAV; belonging to the genus Orthobunyavirus, family Peribunyaviridae) is an arthropod‐borne virus that causes premature birth, abortion, stillbirth and congenital abnormalities with arthrogryposis hydranencephaly syndrome in ruminants (Yanase et al., 2020). This virus may also cause bovine encephalitis and encephalomyelitis in calves and cows (Kono et al., 2008; Liao, Lu, Goto & Inaba, 1996; Oem et al., 2012). There are two major endemic zones of AKAV distribution worldwide, where one extends from East Asia to Southeast Asia to Australia and the other extends from the Middle East to South Africa (Al‐Busaidy et al., 1987; Oğuzoğlu et al., 2015; Taylor & Mellor, 1994; Wang et al., 2019; Yanase et al., 2020).

The AKAV genome, similar to that of other orthobunyaviruses, consists of three negative‐stranded RNA segments, that is S, M and L segments named according to size (small, medium and large, respectively) (Elliott, 2014). The encoded proteins in the genome corresponded with those of other orthobunyaviruses, and the genetic analysis exhibited higher conservation of nucleoprotein (N) on the S segment and RNA‐dependent RNA polymerase (RdRp) on the L segment than those on the M segment, which encodes a polyprotein composed of virion glycoproteins (Gc and Gn) and a nonstructural protein (NSm) (Elliott, 2014). Four genogroups (I–IV) and two lineages in genogroup I (Ia and Ib) were identified through phylogenetic analysis of two full‐length open reading frames of S and M segments and partial sequences of the L segment (Kobayashi et al., 2007).

A few AKAV isolates have been sequenced in Taiwan. PT‐17 was first recognised in 1992 and was isolated from a calf with nonsuppurative encephalitis (Liao, Lu et al., 1996b); its sequence was found to be similar to that of the Iriki strain. NT‐14 was identified in swine tonsils in 2000, indicating that AKAV could naturally infect not only ruminants but also swine (Huang et al., 2003). CY‐77 was isolated from bovine erythrocytes in 1993, and the isolate was found to belong to genogroup Ia (Yamakawa et al., 2006).

Culicoide‐borne diseases are important transboundary diseases. Viral carriers can undergo long‐distance dissemination through wind or monsoons, resulting in spread of the virus from hundreds to thousands of kilometres away from the original outbreak site (Aguilar‐Vega et al., 2019; Braverman & Chechik, 1996; Walker & Klement, 2015). To date, nationwide serosurveillance and regional surveillance of vector specimens have shown that AKAV infection is prevalent in Taiwan (Tzeng et al., 2019). However, the origin of Taiwanese AKAVs and their epidemiology remains unclear.

Accordingly, in this study, we reconstructed the epidemiological history of Taiwanese AKAVs by analysing their phylogenetic relationships to improve our understanding of transboundary Culicoid‐borne diseases in Taiwan and elsewhere.

MATERIALS AND METHODS

Sample source

In total, 15 samples of Taiwanese AKAVs were identified from Culicoides biting midges, mosquitoes, bovine tissues and goat tissues from 1992 to 2015, whereas the other two isolates Tainan/17H8 and 93H78 AKAV were provided by the Council of Agriculture National Institute for Animal Health, Taiwan (Table 1). The locations of all the AKAVs evaluated in this study are labelled on the map of Taiwan shown in Figure 1. The S, M and partial L segments were sequenced in all 17 samples, 14 samples and 12 samples, respectively. The GenBank accession numbers of the genome sequences are as follows: MF278855–MF278882, MT676825–MT676826 and MZ501031–MZ501043.

TABLE 1

Taiwanese AKAVs isolated from 1992 to 2015

Isolate/sample	Source	Year	County	Accession no. of genome segment
				S	M	L
PT‐17	Bovine blood	1992	Pingtung	AF034940
CY‐77	Bovine erythrocytes	1993	Chiayi	AB232319	AB297851	AB297886
Tainan/17H8	Bovine	1993	Tainan	MF278879	MF278865	MZ501032
93H78	Bovine	1993	Tainan	MF278877	MF278864	MZ501031
NT‐14	Swine tonsil	2000	Nantou	AF529883
08H65‐1	Bovine	2008	Hualien	MF278881	MF278862	MZ501033
08H65‐2	Bovine	2008	Hualien	MF278880	MF278863	MZ501034
YL‐1/AKA/C/12/TW	Culicoides spp.	2012	Yunlin	MF278866	MF278855	MZ501035
YL‐2/AKA/C/12/TW	Culicoides spp.	2012	Yunlin	MF278867	MF278857	MZ501036
PT‐1/AKA/C/12/TW	Culicoides spp.	2012	Pingtung	MF278870	MF278857	MZ501038
HL‐1/AKA/C/12/TW	Culicoides spp.	2012	Hualien	MF278868	MF278858	MZ501039
YL‐3/AKA/C/12/TW	Culicoides spp.	2012	Yunlin	MF278871	MF278859	MZ501037
HL‐3/AKA/M/12/TW	Anopheles spp.	2012	Hualien	MF278869
12H37	Bovine	2012	Tainan	MF278873	MF278861	MZ501040
13H33	Bovine	2013	Yunlin	MF278872
13H35	Bovine	2013	Taoyuan	MF278875	MF278860	MZ501042
13H22	Goat	2013	Changhua	MF278876
KH‐1/AKA/C/13/TW	Culicoides spp.	2013	Kaohsiung	MF278882	MT676825	MZ501043
YL‐1/AKA/C/15/TW	Culicoides oxystoma	2015	Yunlin	MF278874	MT676826
YL‐2/AKA/M/15/TW	Anopheles sinensis	2015	Yunlin	MF278878	MT676827

FIGURE 1

The geographic distributions of Akabane orthobunyavirus (AKAV) isolates in Taiwan. The circles (●) represent farms at which tissue specimens were collected, and the squares (■) represent the collection of vector specimens. The triangles (▲) show the counties where PT‐17 (Pingtung), NT‐14 (Nantou), CY‐77 (Chiayi), 93H78 (Tainan) and Tainan/17H8 (Tainan) were isolated. The map was provided by CraftMAP (http://www.craftmap.box‐i.net/)

RNA preparation, polymerase chain reaction amplification and cloning

RNA extraction from vector insects was conducted using the standard protocol from the PureLink RNA extraction kit (Thermo Fisher Scientific, Waltham, MA USA), as described in our previous study (Tzeng et al., 2019). RNA from bovine tissues of infected cattle was obtained from the Animal Health Research Institute. cDNA synthesis was conducted using an iScript cDNA synthesis kit (Bio‐Rad Laboratories, Hercules, CA, USA), and viral genes were amplified using an AccurPrime pfx DNA polymerase kit (Thermo Fisher Scientific). The primers for the S segment were obtained from Yamakawa et al. (2006), and the primers for the M segment were obtained from An et al., and Park (2010). The primers for the L segment were designed from the consensus sequences of AKAV‐7/SKR/2010 (JQ308779), OBE‐1 (NC_009894.1) and DHL10M110 (KY284021) by an alignment result. All primers used are listed in Table S1. The amplicons were analysed on 1% agarose gels, and amplicons with the predicted size were extracted with a Quick Gel Extraction Kit (Thermo Fisher Scientific). To insert the amplicons into the desired plasmid vectors, a 3’ A was added to the amplicons extracted with rTaq DNA polymerase (Takara Bio, Shiga, Japan) at 72°C for 10 min, and the products were then directly ligated into the pGEM‐T easy vector with T4 ligase (Promega, Madison, WI, USA). The transformation and colony selection procedures were described previously (Tzeng et al., 2019).

Sequencing and genome assembly

After the amplicon of each AKAV segment was cloned into a pGEM‐T easy vector, at least three vectors were selected for sequencing, and both sites on the pGEM‐T easy vector using T7 and SP6 primers were implemented for each sample for confirmation. Sequencing was then conducted using a sequencing service (Genomic Ltd., Taipei, Taiwan). The sequencing data were trimmed using VectorNTI 10 Advanced (Thermo Fisher Scientific) to remove plasmids and unqualified sequences. Then, the partial genome sequences of the M segment were assembled into genome contigs using ContigExpress in VectorNTI 10 Advanced (Thermo Fisher Scientific).

Sequence analysis for elucidation of genetic characteristics and phylogeny

Phylogenetic analysis was carried out as described by Kobayashi et al. (2007) and Yanase et al. (2018). Each nucleotide sequence and its deduced amino acid sequence were aligned with downloaded sequences from GenBank (Table S2) using ClustalW in the AlignX program in VectorNTI 10 Advanced (Thermo Fisher Scientific). Construction of a neighbour‐joining tree was performed using the alignment results of each nucleotide sequence with Molecular Evolutionary Genetics Analysis version X or MEGA X (Kumar et al., 2018), and the best nucleotide substitution model was determined using the best‐fit model in MEGA X. The reliability of tree construction was examined using 1000 bootstrap test repeats. The sequences of the Sabo virus were employed in the analysis to determine the relative position of the Taiwanese AKAVs in the phylogenetic tree. All Taiwanese AKAV sequences were uploaded to GenBank, and the accession numbers for the obtained sequences are listed in Table 1.

RESULTS

All Taiwanese AKAVs belonging to genogroup Ia were examined by phylogenetic analyses of the S, M and L segments (Figure S1). Determination of clusters in genogroup Ia relied on two criteria: (1) clusters should be monophyletic groups and (2) partial members in clusters should have relationships identified in prior phylogenetic or spatial–temporal analyses. In total, nine clusters were identified (clusters A–J) in the phylogenetic tree of the S segment (Figure 2a), and the same clusters were also identified in the phylogenetic trees of the M and L segments (Figure 2b and c). The Taiwanese AKAVs were distributed in clusters A and J on the phylogenetic trees. The intracluster nucleotide identities were 98.7–100% in cluster A and 98.3–100% in cluster J; intercluster identities for the S, M and partial L (543 bp) segments were 95.6–96.1%, 95.3–97.7% and 93.5–96.4%, respectively. In addition, the intracluster amino acid identities were 99.1–100% in each cluster; intercluster identities for the S, M and partial L segments were 98.3%–100%, 96.1–100% and 97.7–100%.

Three phylogenetic subtrees of Akabane orthobunyavirus (AKAV) in Taiwan were constructed with the sequences of S, M and L segments. The three comprehensive trees are shown in Figure S1 in the Supplementary Information. (a) The phylogenetic tree was constructed with the 699‐bp full‐length nucleoprotein (N) gene on the S segment using the T92+G nucleotide substitution model. (b) The phylogenetic tree was constructed with the 4206‐bp full‐length Gc‐NSm‐Gn polyprotein gene on the M segment using the GTR+G+I nucleotide substitution model. (c) The phylogenetic tree was constructed with the 543‐bp partial RNA‐dependent RNA polymerase gene on the L segment using the T92+G nucleotide substitution model. Clusters A–J were determined using the S segment‐based tree, and the same clustering is labelled on the M segment‐ and L segment‐based subtrees behind the names of the sequences

The AKAVs within cluster A were isolated since 2012 in Taiwan and in 2012 from Okinawa, Japan. Cluster J consisted of genogroup Ia of AKAVs isolated from mainland Japan in the 1980s and Okinawa, Japan in the 1990s. The AKAVs within‐cluster J were highly related to the Taiwanese isolates, including PT‐17, NT‐14 and CY‐77. The same clustering on the S segment tree was also mapped on the phylogenetic tree of the M and L segments, indicating that there was some consistency in the Taiwanese AKAVs in comparison to the phylogenetic tree of the S segment, with the exception of KH‐1/AKA/C/13/TW (Figure 2). A reassortment AKAV, KH‐1/AKA/C/13/TW, was observed between Taiwanese AKAVs in the cluster A and J and was determined by analysing nucleotide identities using pairwise sequence alignment (Table 2).

TABLE 2

Genetic identities of genomic sequences among KH‐1/AKA/C/13/TW and Taiwanese AKAVs

AKAVs	KH‐1/AKA/C/TW/13
	S segment	M segment	L segment*
In cluster A	96.0–96.4%	99.1–100%	99.2–100%
In cluster J	98.4–100%	95.6–97.5%	93.5–96.4%

Only 1408 bp was used for analysis.

The subtrees of the S, M and L segments exhibited relationships among the occurrence of genogroup Ia AKAVs in East Asia in the past. On the subtree of the S segment, major divergence was observed between cluster J and the other clusters, and cluster J included the earliest genogroup Ia AKAV as the Iriki strain from Honshu, Japan, and others from Okinawa, Japan and Taiwan. Cluster C had a single origin in cattle experiencing the same encephalomyelitis epidemic in Kagoshima and Kumamoto, Japan in 2006 (Kono et al., 2008). Cluster D also consisted of the AKAV isolated from cattle with encephalomyelitis and that from South Korea in 2010 and Japan in 2011 (Oem et al., 2012; Yanase et al., 2018). Cluster E was a small cluster composed of AKAVs from Okayama, Kagoshima and Yamaguchi, Japan. Clusters F contained viruses from Honshu and Okinawa, Japan isolated prior to 2001. Cluster G consisted of all isolates obtained from China in 2010, 2013, 2016 and 2018 (Cao et al., 2019; Tang et al., 2017; Tang et al., 2019) as well as isolates from Honshu, Japan obtained in 2013 (Yanase et al., 2018) Cluster H consisted of early isolates in Honshu in 1984 and the coming isolates in Okinawa in 1993 (Figure 2a). In the subtree of the M and L segments, two characteristics were observed when compared with the subtree of the S segment: first, the relative positions of clusters A, B and C were consistent in all subtrees, and second, the members of cluster G were separated into two groups in both subtrees, although on the M segment subtree (Figure 2b and c).

After the publication of the Iriki strain, which was the earliest genogroup Ia strain, the complete open‐reading frame of the S segment has been sequenced in 94 isolates or samples from genogroup Ia. These viruses are all distributed in Taiwan, Japan, South Korea and China (Table S2). In Taiwan, AKAVs, which were first sequenced in 1992, have become endemic viruses, and the new more invasive virus has become endemic since its initial identification in 2012. Moreover, the Taiwanese AKAVs were found to be highly correlated with the Okinawa isolate.

DISCUSSION

The genetic divergence of AKAV clusters (clusters A and J) observed in our analytical results indicated that at least two invasive events established the current AKAV populations in Taiwan. The serological prevalence of AKAV in dairy farms indicated an 89% infection rate in 1989 (Liao et al., 1996a); however, the first S segment of Taiwanese AKAV, that is, PT‐17, from cattle with encephalitis was isolated and sequenced in 1992 (Liao et al., 1996b). From these findings, it was still unclear whether the Taiwanese AKAVs found in 1992 were caused by the resurgence of the endemic AKAV prevalent before 1989 or represented a new invasive strain. Our data combined with records supported that a single invasion occurred before 1989 to establish Taiwanese AKAVs (cluster J). Based on Liao et al. (1996a) that observed a sustained epidemic of AKAVs in Taiwan from 1989 to 1994 that corresponded to PT‐17, 93H78 and CY‐77 were isolated in 1992 and 1993, which were grouped in cluster J. This indicated the isolates from 1989 to 1994 might have come from the same epidemic event, although AKAVs did not be isolated or sequenced in 1989 and 1994. In addition, Taiwanese AKAVs (cluster J) formed a monophyletic group on the phylogenetic trees of the S and L segments, indicating that they may have originated from a single strain and then spread to the whole island of Taiwan, although the members of cluster J did not form a monophyletic group in the phylogenetic tree of the M segment. A possible reason for this may be that members of cluster J included Japanese AKAVs isolated from 1984 to 1985. Our phylogenetic analysis showed that these strains had multiple origins; however, the M and partial L segments were sequenced for only two of these isolates, suggesting that insufficient information was available for analysing the M segment of cluster J. Despite this, the phylogenetic results obtained using sequences of the S and L segments still supported that there was a single invasion of Taiwanese AKAVs (cluster J). Thus, the population of Taiwanese AKAVs (cluster J) may have been established via a single invasion event before 1989. In the case of Taiwanese AKAVs isolated in 2012 (cluster A), the Okinawa isolate ON‐1/P/12 was grouped with Taiwanese AKAVs from 2012 to 2015 in cluster A on the trees of the S, M and L segments, indicating that those of Okinawa and Taiwan had the same origin, although that specific origin has not yet been determined. Overall, our findings supported that the population of Taiwanese AKAVs was established through a single invasion event that may have occurred in 2012 and before 1989; however, the data were still insufficient to conclusively determine the origins of Taiwanese AKAVs.

Our results revealed that Taiwan and the islands in Okinawa Prefecture shared AKAV epidemics. Indeed, the geographic distance between Okinawa and Taiwan is closer than that between Okinawa and mainland Japan (Figure 3). Notably, some AKAVs from Ishigaki Island and Irimote Island in the Yaeyama Islands were found to have closer genetic features to Taiwanese AKAVs, indicating that these isolates may have originated from Taiwan (Figures 2a and 3). However, most AKAVs of Okinawa are genetically similar to those from mainland Japan (Yamakawa et al., 2006), implying that the islands of Okinawa may be an intersection for transboundary transmission of AKAVs from both Taiwan and mainland Japan. The long‐distance dissemination of Culicoides biting midges by wind is determined to involve in the transmission of bluetongue virus, bovine ephemeral fever virus and AKAV (Aguilar‐Vega et al., 2019; Braverman & Chechik, 1996; Walker & Klement, 2015); however, further studies are needed to elucidate the migration patterns of Culicoides biting midges among Taiwan and the islands of Okinawa Prefecture.

FIGURE 3

Proposed route of Taiwanese‐related AKAVs in clusters A and J. The map was generated by QGIS

The geographic distributions of AKAV genogroups exhibited indistinct boundaries. Genogroup Ia is prevalent in southern China and Taiwan but sporadically appears in Japan and South Korea (An et al., 2010; Cao et al., 2019; Kono et al., 2008; Liang et al., 2018; Tang et al., 2017; Tang et al., 2019; Tzeng et al., 2019; Yamakawa et al., 2006; Yanase et al., 2018). Genogroup Ib has been detected sporadically in Israel, Turkey and Indonesia (Oğuzoğlu et al., 2015; Purnomo Edi et al., 2017; Şevik, 2017; Stram et al., 2004; Yamakawa et al., 2006). Furthermore, genogroup II is endemic to Japan and South Korea (An et al., 2010; Kobayashi et al., 2007; Yamakawa et al., 2006), whereas genogroup III only exists in Australia, and genogroup IV is found only in Africa (Wang et al., 2019; Yamakawa et al., 2006). In East Asia, all provinces of China have experienced AKAV epidemics (Wang et al., 2017); however, AKAV sequences have only been obtained from the southern provinces of Yunnan, Guangdong, Guangxi and Hainan. Our results showed that the isolate GD18134, which was from Guangdong, China, was phylogenetically related to isolates from Japan in 2013 (Wu et al., 2020; Yanase et al., 2018) (Figures 2 and 3a). Further analyses of AKAV sequences from China and Southeast Asia are essential to improve our understanding of the temporal and spatial patterns of transboundary transmission of viruses. These studies will also provide insights into the origins of Taiwanese AKAVs. Herein, we provided the first comprehensively phylogenetic data to describe the interaction between AKAVs in Taiwan and nearby islands, also filling up a blank for understanding about molecular epidemiology of AKAV in East Asia.

AUTHOR CONTRIBUTIONS

H. Y. Tzeng and W. C. Tu devised the project, the main conceptual ideas and proof outline. C. L Tsai provided the details of phylogenetic techniques. H. Y Tzeng and K. M. Liao collected all haematophagous insects from sampling farms. L. J. Ting collected the fetus samples and extracted RNAs from the samples. H. Y Tzeng work for cloning and sequencing of target viral sequences. H. Y. Tzeng and W. C. Tu wrote the manuscript.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

ETHICS STATEMENT

The conduction of all experiments followed by the guideline of experimental animal management in Taiwan.

PEER REVIEW

The peer review history for this article is available at https://publons.com/publon/10.1002/vms3.887.

Figure S1. Comprehensive phylogenetic trees of the (a) S, (b) M and (C) L segment

Click here for additional data file.

Table S1. Primer sets for genome amplification of AKAVs in this study

Click here for additional data file.

Table S2. Sequences of AKAVs used for phylogenetic analyses in this study

Click here for additional data file.