Two Rhabdoviruses, One Novel, Isolated from Armigeres subalbatus in China

The family Rhabdoviridae contain important human and mammalian pathogens that are vectored by different arthropod species. The ground supernatants of mosquitoes were used to inoculate in BHK-21 and C6/36 cells for virus isolation. Then, the viral complete genome sequence was obtained and used for phylogenetic analysis. In this study, we observed a cytopathic effect (CPE) in mosquito cells (C6/36) and rod-like virion after inoculating a pool of Armigeres subalbatus samples collected in Shanxi Province, China, in 2019 (SX1916). Meta-transcriptomics sequencing revealed the presence of two distinctive rhabdoviruses with similar abundance levels, namely, Shanxi Armigeres subalbatus rhabdovirus (SXARV) and Shanxi Arboretum virus (SXABTV). Despite the fact that the SXARV genome (9590 nt) was much shorter than that of SXABTV (11,480 nt), both belonged to the Almendravirus group within Rhabdoviridae whose genomes encoded five proteins (N, P, M, G, and L) and a small hydrophobin (U1) and the difference in lengths is mainly caused by a substantially shorter N protein encoded by SXARV. On the phylogenetic tree, SXABTV was closely related (90.7% amino acid identity at L protein) with the Arboretum virus isolated from Psorophora albigenu mosquitoes in Peru in 2014, whereas SXARV was distantly related to Rio Chico virus (63.3% amino acid identity), a genetic distance large enough to be defined as a new species within Rhabdoviridae. Collectively, we report a simultaneous isolation of two related rhabdoviruses from Armigeres subalbatus that marked the circulation of almendraviruses in Shanxi, China.


Introduction
Rhabdovirus is a single-stranded negative-sense RNA virus that is rod or bullet shaped, with a particle length of 100-430 nm and diameter of 45-100 nm. According to the latest virological classification report from the International Committee on Taxonomy of Viruses (ICTV), the family Rhabdoviridae comprises 45 genera and 275 species [1]. Rhabdoviridae viruses have a wide host range, including arthropods, such as mosquitoes, ticks, sandflies, birds, fish, and mammals [2][3][4].
Almendravirus (Rhabdoviridae) virus had a full-length genome of approximately 11,000 nt, which mainly encode five proteins, namely, nucleoprotein (N), phosphoprotein (P), and matrix protein (M), glycoprotein (G), RNA-dependent RNA polymerase protein (L), as well as a small hydrophobin (U1) between the G and L genes, following the order 3 N-P-M-G-U1-L5 [4]. The U1 protein in the viral genome has unique characteristics that are used to identify Almendravirus viruses [5], which have been isolated in North America, South America, and Asia [5][6][7].
To date, the ICTV has established seven viral species within the genus Almendravirus: Arboretum virus and Puerto Almendras viruse isolated from Peru [7], Balsa virus isolated from Colombia [5], Coot Bay virus isolated from USA [5], Menghai rhabdovirus isolated from China [6], and Japan Rhabdovirus [8] and Rio Chico virus isolated from Panama [5]. Among these, only Menghai rhabdovirus was extensively described in China, including a strain isolated from Aedes albopictus collected in 2010 [6] and five strains isolated from midges collected in 2018. In addition, Almendravirus has not been isolated from bloodsucking insects collected in China.
From June to August 2019, we conducted comprehensive sampling of different species of vectors, including sandfly and mosquito in Wuxiang County, Shanxi Province, Central China. A number of viruses have been identified from these vectors, including Wuxiang virus isolated from sandfly samples [9]. In this study, we report two virus strains isolated from Armigeres subalbatus, both belonging to the Almendravirus genus (Rhabdoviridae).

Virus Isolation and Preliminary Identification
The collected mosquito specimens were pooled into 16 samples and homogenized. Each pool contains 30~50 mosquitoes and there are a total of seven pools for Anopheles sinensis, seven for Armigeres subalbatus, one for Aedimorphus vexans, and one for Culex pipiens pallens. The resulting supernatant was inoculated onto BHK-21 or C6/36 cells and the cells were cultured continuously and observed daily for CPE. For the pooled Armigeres subalbatus sample SX1916, the first passage in C6/36 cells showed no CPE. CPE appeared in the second generation on day 4, characterized by the aggregation of large numbers of cells and cell space enlargement ( Figure 1B). The cell supernatant showing CPE was re-inoculated onto C6/36 cells for passages 3 and 4, and CPE again appeared on day 4, suggesting that CPE caused by SX1916 can be stably passaged. On the other hand, the supernatant from SX1916 isolate was inoculated onto BHK-21 cells for three consecutive generations and no CPE was observed. To identify the infectious agent that causes CPE, we first used multiple arbovirus group/genera-specific gene amplification primers (flavivirus, alphavirus, and bunyavirus) and virus species-specific gene amplification primers (Japanese encephalitis virus, Banna virus, Tibet orbivirus, Densovirus, and Totivirus) for gene amplification testing via reversetranscription (RT)-PCR of the SX1916 virus isolate (Table S2). None of the primers used resulted in positive results.
To identify the viruses from the SX1916 isolate, we extracted the total RNA of the supernatant of C6/36 cells infected with the SX1916 isolate for next-generation sequencing (NGS). The sequencing resulted in 9,144,290 reads, which were assembled into 4969 contigs. Taxonomy annotation of all contigs revealed two rhabdovirus genome sequences in the SX1916 isolate, with lengths of 11,480 and 9590 nt, respectively. These viruses were tentatively named as Shanxi Arboretum virus (SXABTV, strain SXABTV1916-1) and Shanxi Armigeres subalbatus rhabdovirus (SXARV, strain SXARV1916-2). Importantly, both rhabdoviruses had very high mean coverage, namely, 33,038-fold for SXABTV and 76,101-fold for SXARV (Figure 2), which account for 28.5% and 5.6% of total RNA sequenced, suggesting the two viruses are highly abundant in the cell supernatant. In addition, we searched for other potential virus contigs from our annotation results (Table S3). However, other virus sequences contain only Guadeloupe mosquito quaranja-like virus 3 (423 bp) and Guato virus (474 bp) and they are too short for analysis.

Virus Plaque Purification
Because the NGS results identified two rhabdoviruses with large genomic differences in the SX1916 isolate, we applied virus plaque purification methods to purify each rhabdovirus. The SX1916 isolate was used for four consecutive passages of the fourth-passage C6/36 cells for virus plaque purification. In the first and second plaque purification, the virus-infected C6/36 cells were cultured for 7 days (a short-diameter plaque was observed on day 6) to obtain a single virus plaque. We obtained nine and three single virus plaques in the first and second purification tests, respectively ( Figure S1). All 12 plaques had similar diameters (mean, 2.4 ± 0.5 mm; n = 12); they were thawed and diluted with cell culture medium to extract RNA and subjected to gene detection. The results of parallel detection of SXARV and SXABTV genes showed that all 12 plaques were SXARV gene positive and SXABTV gene negative. Plaques that tested positive for SXARV were diluted with cell culture medium and subjected to additional plaque purification. After three rounds of plaque purification, an SXARV isolate (SXARV1916-2) was obtained. All subsequent studies involving SXARV came from the SXARV1916-2 strain.
To confirm that sample used before the plaque purification test contained both viruses, we performed parallel genetic testing of both viruses in fourth-generation C6/36 cells inoculated with the SX1916 isolate to verify whether both viruses were present in the cell supernatant used in the plaque purification test. The results showed simultaneous amplification of the viral gene amplification products of both SXARV and SXABTV in the fourth-generation supernatant of the C6/36 cells ( Figure S2).
The plaque-purified SXARV1916-2 isolate was inoculated into C6/36 cells to observe its cytopathic characteristics. In second-generation C6/36 cells inoculated with plaquepurified SXARV1916-2, CPE was clearly observed after 4 days of culture, manifesting as large numbers of fused cells and cell shedding ( Figure 1C). However, in fourth-generation C6/36 cells inoculated with the SX1916 virus isolate, CPE on day 4 was characterized by the aggregation of large numbers of cells, but not cell shedding ( Figure 1B). Thus, CPE differed significantly between C6/36 cells inoculated with purified SXARV1916-2 and those inoculated with the SX1916 virus isolate containing two virus strains.

Electron Microscopy of Virus Morphology
Plaque-purified SXARV was inoculated into C6/36 cells and ultra-thin section detection was performed when the cells showed clear CPE. Rod-shaped virus particles (length, 300 nm; diameter, 70 nm) were observed via electron microscopy ( Figure 1D and Figure S3).

Nucleotide Sequence of the Viral Genome Coding Region
According to the virus genome sequence information obtained via NGS sequencing, we designed amplification primers to cover the full-length genome nucleotide sequence of SXARV and SXABTV, respectively, using Primer-BLAST (https://www.ncbi.nlm.nih.gov/ tools/primer-blast (accessed on 5 April 2021)), and the primers used are listed in Table S4. The cDNA library of the fourth-generation supernatant of C6/36 cells inoculated with the SX1916 isolate was used as a template. PCR amplification and nucleotide sequence determination of the whole genome of each rhabdovirus were conducted using the designed primers. The non-infected control C6/36 cell lines were negative for both viruses discovered in this study. The obtained whole-genome nucleotide sequences of SXABTV and SXARV were uploaded to GenBank under accession nos. MW890015 and MW890016, respectively.
Sequence analyses of SXABTV and SXARV genomes showed that both viruses encode six proteins: N proteins of 1305 nt/434 aa and 471 nt/156 aa, P proteins of 840 nt/279 aa and 636 nt/211 aa, M proteins of 480 nt/159 aa and 480 nt/159 aa, G proteins of 1431 nt/476 aa and 1350 nt/449 aa, and L proteins of 6189 nt/2062 aa and 6189 nt/2062 aa. The viral genomes of SXABTV and SXARV also contained genes U1 (222 nt/73 aa and 162 nt/53 aa) between the G and L proteins of the small hydrophobin ( Figure 3). The difference in length between SXABTV and SXARV is mainly caused by a substantially shorter N protein encoded by SXARV.

Nucleotide and Amino Acid Sequence Identity in Viral Gene Coding Regions
We examined the identity of nucleotides and amino acids in the coding regions of SXABTV and SXARV and other Almendravirus viruses, including 13 isolate strains in GenBank, and found that SXABTV and Arboretum virus (ABTV) had the highest identity. We identified nucleotide (amino acid) identity of the coding region of the N gene at 83.9% (97.5%), the P gene at 70.8% (68.5%), the M gene at 79.3% (88.5%), the G gene at 75.6% (84.1%), the U1 gene at 69.4% (52.9%), and the L gene at 78.9% (90.7%).
Nucleotide and amino acid identity of the coding regions of the SXARV isolate had the highest identity with the Rio Chico virus. We identified nucleotide (amino acid) identity of the coding region of the N gene at 59% (52.8%), the P gene at 55% (54.8%), the M gene at 59.5% (48.1%), the G gene at 53.4% (41.7%), the U1 gene at 57.5% (49.4%), and the L gene at 64.6% (63.3%) ( Table 1).
The virus nucleotide and amino acid identity data (

Evolutionary History of SXABTV and SXARV
To analyze the phylogenetic relationship between the SXABTV and SXARV and other Almendravirus viruses, we prepared an analysis dataset of the L gene sequence of Rhabdoviridae viruses, including all Almendravirus viruses and representative strains of 29 other genera of Rhabdoviridae (Table S5). Phylogenetic analysis results obtained using the maximum likelihood method showed that both SXABTV and SXARV were in the Almendravirus genus, and that SXABTV was closely related to an ABTV isolate obtained from mosquitoes collected in Peru, whereas SXARV was closely related to a Rio Chico virus isolate obtained from mosquitoes collected in Panama (Figure 4).

Discussion
In this study, we obtained virus isolate SX1916 from Armigeres subalbatus collected in Shanxi Province in Central China in 2019. SX1916 was found to cause CPE in C6/36 cells and was stably passaged ( Figure S5). Virology and molecular biology analyses showed that SX1916 contained two rhabdoviruses, SXABTV and SXARV.
Plaque purification can be used to clone virus isolates or variants with different morphological characteristics, including plaque diameter [10]. We applied this technique to clone separate isolates of SXARV and SXABTV from the SX1916 isolate. However, all of the plaques obtained in two plaque purification assays belonged to SXARV, and no SXABTV plaques were obtained, which suggested that SXABTV may not be capable of forming plaques. A study on the plaque-forming ability of 212 arboviruses (including serotypes) in Vero and LLC-MK2 cells showed that 169 arboviruses can form plaques in these two cell lines, whereas 17 and 18 arboviruses could not form plaques in VERO and LLC-MK2 cells, respectively, and the remaining seven could not form plaques in either cell line [11]. The results of that study suggested that some important mosquito-borne viruses, such as Zika virus, can cause CPE, but not form plaques, in Vero cells [12]. Thus, not all arboviruses have the ability to form plaques. ABTV can cause CPE in C6/36 cells, but whether it can form plaques in this cell line has not been reported [7]. The SXABTV isolate failed to form plaques in C6/36 cells; however, further research is required to conclude that it has no plaque-forming properties. It is possible that the presence of SXARV inhibited plaque forming by SXABTV in the present study. Interference between the two rhabdoviruses may also have affected their cytopathic properties. After the SX1916 isolate was inoculated with C6/36 cells, the cells showed significant aggregation, but not cell shedding. However, after plaque-purified SXARV was inoculated onto C6/36 cells, these cells showed substantial cell fusion and shedding. These differences may be caused by mutual interference between the co-infecting SXABTV and SXARV, which was then eliminated in plaque-purified SXARV. A study of snakehead retrovirus (Sn RV) and grouper nervous necrosis virus (GNNV) co-infection of GF-1 cells showed that Sn RV enhanced the infectivity and CPE of GNNV [13], suggesting a complicated co-infection mechanism driven by interaction between the two viruses.
Viruses in the Almendravirus genus encode the N, P, M, G, and L proteins. Between the G and L proteins, they also encodes a small hydrophobin (U1) [7]. The conservation of amino acids in the coding region of the rhabdovirus L gene is higher than that of the other four proteins (N, P, M, and G) [6]; RdRp in the L gene coding region is a core protein that determines viral genome replication [7]. Recent studies have reported that amino acid identity of the coding region of the L protein of <90% is sufficient for identifying new rhabdovirus species [14].
As of 10 March 2021, GenBank contained the genome sequences of 13 Almendravirus species (Table S5), including the 2 rhabdoviruses isolated in this study; a total of 15 Almendravirus viruses have been sequenced. Our amino acid identity analysis of the L protein coding region of these 15 virus isolates showed 95.9-100% identity among the seven Menghai rhabdoviruses isolated in Yunnan and Guangdong Provinces in China, 99.6% identity between the two Balsa viruses isolated from Peru, and 90.7% identity between SXABTV isolated in this study and the ABTV isolated from Peru. Thus, these 11 virus isolates belonged to three Almendravirus species (Figure 4), among which 7 isolates belonged to Menghai rhabdovirus virus and 2 to Balsa virus. The nucleotide and amino acid lengths of the coding regions of the six genes from SXABTV isolated in this study were identical to those of ABTV; SXABTV had only one more amino acid T than ABTV, at the 188-190 position of the P gene. The amino acids of the L protein coding region had 90.7% identity between SXABTV and ABTV, indicating that SXABTV was an isolated strain of ABTV, which is found in different regions and different mosquito species.
SXARV isolated in this study had a full-length genome of 9590 nt, whereas that of SXABTV was 11,480 nt, and the nucleotides and amino acids of the coding regions of the six protein genes differed significantly between the two virus isolates (Table 1), with only 45.5% amino acid identity in the L gene coding region between the two viruses. The SXARV and SXABTV isolates obtained in this study clearly belong to two completely different virus species. Interestingly, although the results of our phylogenetic analysis showed that SXARV and the Rio Chico virus isolated from the Panamanian mosquito were of the same evolutionary branch, the total length of the SXARV isolate was different from known Almendravirus viruses, including Rio Chico virus. All Almendravirus viruses except SXARV had a genome length between 10,735 nt (Rio Chico virus) and 11,876 nt (Puerto Almendras virus) ( Table 1), whereas the full length of the SXARV genome was only 9590 nt, representing the shortest genome among known Almendravirus viruses. The difference in length between SXARV and other Almendravirus viruses is mainly caused by a substantially shorter N protein encoded by SXARV. The highest identity of the L gene of SXARV was 63.3%, with a Rio Chico virus isolate. Since this value is <90%, SXARV can be considered a new rhabdovirus species [14]. In summary, SXARV is a novel virus species in the Almendravirus genus (Rhabdoviridae).

Specimen Collection
From June to August in 2019, we collected blood-sucking insects overnight (18:00 p.m. to 6:00 a.m.) at animal pens including chicken pens, dog houses, and mule pens in Wuxiang County, Shanxi Province, China (36 • 39 -37 • 8 N, 112 • 26 -113 • 22 E) using a blood-sucking insect collection tool (model no.: MM200BL; Guangzhou Changsheng Chemical Technology Service Co., Ltd., Guangzhou, China, https://guangzhou.11467.com/info/8049209.htm (accessed on 5 June 2019)) comprising a 12-V, 12-Ah battery, 0.1-A bulb, and 0.14-A fan. All mosquito specimens were placed in a −40 • C low-temperature refrigerator for 30 min and then removed for processing in an ice bath. Each specimen was labeled according to morphology, collection time, collection environment, and then stored in liquid nitrogen until laboratory testing.

Virus Isolation
To isolate viruses from mosquitoes, we pooled groups of 30-50 individuals according to collection site and species, put them into a glass grinder, and washed them twice with grinding solution consisting of 93% Eagle medium, 5% penicillin and streptomycin (100 U/mL), 1% glutamine (30 g/L), and 1% NaHCO 3 . We then added 1.5 mL grinding solution and ground the pooled sample in an ice bath. After centrifugation at 4 • C and 12,000 rpm for 30 min, 70 µL of ground supernatant was inoculated into 80% monolayer BHK-21 and C6/36 cells in a 24-well plate (Corning Inc., Corning, NY, USA). BHK-21 and C6/36 cells were continuously cultured in an incubator at 37 • C with 5% CO 2 and at 28 • C, respectively. Cytopathic effect (CPE) was monitored under a microscope every 12 h. When CPE appeared, the virus solution was collected and stored in the refrigerator at −80 • C until identification. First-generation cell culture supernatant that did not show CPE was blind passaged for three generations and discarded if CPE did not appear after that [15][16][17].

Virus Plaque
C6/36 cells were transferred into 6-well culture plates (Corning Inc., NY, USA) to grow into monolayer cells with a coverage rate of 80%. Virus dilutions (10 -1 -10 -6 ) were sequentially added to the 6-well culture plate (0.1 mL/well). The wells were placed in a 28 • C incubator for 1 h, and 1% agarose-MEM medium containing 2% FBS was added to each well to cover the cells (3 mL/well). After culturing for 5 days, a second layer of 1% agarose-MEM medium containing 7% neutral red and 2% FBS was spread over the cells (3 mL/well). The time and number of instances of plaque appearance were observed, and the plaque diameter was measured as previously described [18].

Electron Microscopy
When CPE reached 75%, virus-infected C6/36 cells were collected by centrifugation, and the sample was fixed with 2% formaldehyde and 2.5% glutaraldehyde solution. Following stepwise dehydration, soaking, embedding, polymerization, and trimming, 80 nm sections were prepared using an ultra-thin slicer and dried at room temperature for staining. The sections were stained with 1% uranyl acetate for 10 min and cleaned with ultrapure water. Next, the slices were stained with 0.02% lead citrate for 5 min and then cleaned and dried at room temperature for inspection. Observations were conducted using a transmission electron microscope (model no.: TF20; FEI Co., Hillsboro, OR, USA) [19].

Viral RNA Extraction and cDNA Library Preparation
Total RNA was extracted from the supernatant of virus-infected cells using the QI-Aamp Viral RNA Mini Kit (Qiagen, Valencia, CA, USA) according to the manufacturer's instructions. The extracted RNA was placed in a 65 • C water bath for 10 min and then immediately placed in an ice bath for 2 min. Next, 32 µL RNA was added to the first-strand reaction tube of the Ready-To-Go kit (GE Healthcare, Little Chalfont, Buckinghamshire, UK), which was maintained at room temperature for 1 min before the addition of 1 µL pd(N)6 random primer (TaKaRa, Kyoto, Japan). This mixture was centrifuged immediately and incubated in a water bath at 37 • C for 1 h. The total volume of the cDNA library was 33 µL; it was stored at −40 • C until later use [19,20].

Virus Identification Using Total RNA Sequencing
To examine the viral genome of the SX1916 isolate, we extracted the total RNA of the supernatant of C6/36 cells infected with the SX1916 isolate for next-generation sequencing (NGS) using the QIAamp Viral RNA Mini Kit (QIAamp, Qiagen, Valencia, CA, USA). We used the TruSeq total RNA library kit to construct the cDNA library, and sequencing was performed on the HiSeq4000 platform (Illumina, San Diego, CA, USA). After low-quality reads were removed from the sequencing data, we performed de-novo assembly. We extracted virus-related contigs based on the taxonomic information of the blast hits. To avoid mis-assembly, reads were mapped back to the virus genome with Bowtie2 (version 2.3.5.1, Ben Langmead) and inspected with Geneious Program (version 9.1.5, Auckland, New Zealand) [21]. The viral genome sequences obtained by NGS was used as template to design overlapping PCR primers for sequence confirmation.

Viral Gene Amplification and Nucleotide Sequencing
We applied polymerase chain reaction (PCR) for viral gene amplification in a 25 µL system comprising 2 µL cDNA as template, 12.5 µL 2× GoTaq Green Master Mix, 8.5 µL nuclease-free water, and 1 µL each of 10 pmol/µL upstream and downstream primers. After the PCR was complete, 5 µL of gene amplification products was detected by 1% agarose gel electrophoresis [19,20]. The primers used in this study are listed in Tables S2 and S3. The positive PCR products were sequenced following the four-color fluorescent labeling dideoxy termination method (Tianyi Huiyuan Biotechnology Co., Ltd., Beijing, China).

Nucleotide Sequence Analysis
The nucleotide sequence was compared to non-redundant protein (nr) database available at the National Center for Biotechnology Information (NCBI) website (https: //blast.ncbi.nlm.nih.gov/Blast.cgi (accessed on 5 May 2021)) using blastn program. We used the Seqman software (DNAStar, Madison, WI, USA) for nucleotide sequence splicing and quality analysis and the BioEdit software (Version 7.0, Tom Hall) for multiple sequence alignment. We used the MEGA software (Version 7.0, Kumar, Stecher, and Tamura 2015) to complete the systematic evolution analysis based on the maximum likelihood method, with a bootstrap value of 1000; the General Time Reversible model is used in analysis. We used the MegAlign (Madison, WI, USA) tool for identity analysis of nucleotide and amino acid sequences [16,19]. Strain information used in virus molecular genetic evolution analysis is provided in Table S5.

Conclusions
In the study, two rhabdoviruses, SXARV and SXABTV, were found in a pool of Armigeres subalbatus samples collected in Shanxi Province, China. Amino acid identity analysis shows that SXARV is a novel virus species in the Almendravirus genus (Rhabdoviridae). Despite the isolation of a variety of Almendravirus viruses in Asia, North America, and South America, it remains unclear whether these viruses could infect mammalian hosts. Therefore, further investigations of the infection status of Almendravirus in domestic and wild animals will be of great importance.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/ 10.3390/pathogens11060624/s1, Table S1: Summary of mosquito specimens collected from June to August 2019 in Wuxiang County, Shanxi Province; Table S2: Viral gene amplification primers used for initial virus identification in this study; Table S3: Taxonomy annotation of all virus contigs. Table  S4: Shanxi Armigeres subalbatus rhabdovirus (SXARV) and Shanxi Arboretum virus (SXABTV) genome nucleotide sequence amplification primers used in this study; Table S5: Virus strains used in virus molecular genetic evolution analysis in this study; Figure S1: Plaque purification of the SX1916 virus isolate; Figure S2: SXARV and SXABTV gene amplification in fourth-generation C6/36 cells inoculated with the SX1916 virus isolate; Figure S3. Electron microscopic morphology of SXARV. Figure S4: Alignment of Amino acid sequence of the Shanxi Arboretum virus (SXABTV) and Shanxi Armigeres subalbatus rhabdovirus (SXARV). Figure S5. SX1916 virus (SXABTV and SXARV) and SXARV gene amplification.