An Investigation of Culicoides (Diptera: Ceratopogonidae) as Potential Vectors of Medically and Veterinary Important Arboviruses in South Africa

Culicoides-borne viruses such as bluetongue, African horse sickness, and Schmallenberg virus cause major economic burdens due to animal outbreaks in Africa and their emergence in Europe and Asia. However, little is known about the role of Culicoides as vectors for zoonotic arboviruses. In this study, we identify both veterinary and zoonotic arboviruses in pools of Culicoides biting midges in South Africa, during 2012–2017. Midges were collected at six surveillance sites in three provinces and screened for Alphavirs, Flavivirus, Orthobunyavirus, and Phlebovirus genera; equine encephalosis virus (EEV); and Rhaboviridae, by reverse transcription polymerase chain reaction. In total, 66/331 (minimum infection rate (MIR) = 0.4) pools tested positive for one or more arbovirus. Orthobunyaviruses, including Shuni virus (MIR = 0.1) and EEV (MIR = 0.2) were more readily detected, while only 2/66 (MIR = 0.1) Middelburg virus and 4/66 unknown Rhabdoviridae viruses (MIR = 0.0) were detected. This study suggests Culicoides as potential vectors of both veterinary and zoonotic arboviruses detected in disease outbreaks in Africa, which may contribute to the emergence of these viruses to new regions.


Introduction
Culicoides midges have been implicated in the transmission of various pathogens of medical and veterinary importance [1]. The unprecedented spread of bluetongue virus (BTV) (Reoviridae: Orbivirus) [2] followed by the emergence of Schmallenberg virus (SBV) (Peribunyaviridae: Orthobunyavirus) into northern Europe [3] and the high incidence of African horse sickness in South Africa [4] highlights the importance of Culicoides midges (Diptera: Ceratopogonidae) as arbovirus vectors in establishing epidemics in animals. To date, over 50 arboviruses have been isolated from Culicoides, 45% of which have not been detected in other arthropods [1]. A recent study demonstrated human feeding to be as high as 11% in field collected Culicoides species, and as yet, little information exists on the role of Culicoides as vectors for zoonotic pathogens [5].
Biting midges have been implicated as the primary vectors of SBV and Oropouche virus (OROV) (Peribunyaviridae: Orthobunyavirus) which are responsible for congenital deformities in sheep and cattle and from mild febrile illness to severe disease in humans [3,6]. Arboviruses with zoonotic potential such as West Nile virus (WNV) (Flaviviridae: Flavivirus) [7], SHUV (Peribunyaviridae: Orthobunyavirus) [8], Dugbe virus (Nairoviridae: Nairovirus) [8], and Rift Valley fever virus (RVFV) (Phenuiviridae: Phlebovirus) [9] have also been isolated or detected in Culicoides pools, although the role of midges in the epidemiology of these viruses is unknown. Culicoides midges were separated from other insects using a light microscope and pooled (100 ≤ n ≤ 500), according to site and month, and homogenized using a TissueL-yser™ (QIAGEN, Valencia, CA, USA). Following maceration, RNA was extracted using the QIAamp ® Viral RNA mini kit (QIAGEN), according to the manufacturer's instructions. During the summer months, when insect numbers were high, testing was limited to only two randomly selected pools per site per month. Whereas during winter months, insect numbers declined, therefore, only one pool was available and tested and, in some cases, no insects were collected, and no pools tested (n < 100). Pool sizes were, therefore, determined based on the number of midges collected per site, per month, with a maximum of 500 and a minimum of 100 midges per pool.
Various published reverse transcription polymerase chain reactions (RT-PCRs) were used to screen pools for arboviruses in the genera Alphavirus [18], Flavivirus, [19] Orthovirus [20], and Phlebovirus [20], as well as equine encephalosis virus (EEV) [14] and Rhabdoviridae [21]. All orthobunyavirus PCR positives were followed up with a One-Step Bunyaviridae RT-PCR. The One Step Bunyaviridae conventional RTPCR was used to obtain a > 500 nt gene regions for known orthobunyavirus positives. SuperScript ® III/Platinum ® Taq Mix (ThermoFisher Scientific, MA, USA) was used in combination with 20 pmol each of BUN 1 and BUN 2 primers [22] in a 50 µL reaction. These primers amplify a 550 base pair fragment of the N gene of the S segment of orthobunyaviruses. Cycling conditions were as follow: 60 • C for 1 min, 50 • C for 30 min, 94 • C for 15 min, followed by 35 cycles of 94 • C for 45 s, 50 • C for 45 s, and 72 • C for 45 s, and a final extension of 10 min at 72 • C [12]. All SHUV positives were followed up with a SHUV M segment PCR (Appendix A). Good laboratory practice with separated pre and post PCR areas, as well as the inclusion of negative and positive controls, were used in all assays. PCR positive results were confirmed by sequencing at Inqaba Biotec (Pretoria, South Africa) and sequence data were analysed and trimmed using CLC Main Workbench version 8.0.1 (https://www.qiagenbioinformatics.com (accessed on 19 September 2019)) and MEGA 6.06 software (https://www.megasoftware.net (accessed on 20 September 2019)). Multiple sequence alignments were compiled using the online version of MAFFT version 7 software (http://mafft.cbrc.jp (accessed on 19 September 2019)). All assembled sequences were compared to the National Center for Biotechnology Information (NCBI) using the Basic Local Alignment Search Tool (BLAST) (http://www.ncbi.nlm.nih.gov/BLAST/ (accessed on 19 September 2019)). Sequences were also compared to positive controls to exclude false positives. The sequences >300 nucleotides (nt) were submitted to GenBank (Accession numbers MN270991-MN271022, Table S2). Maximum likelihood analyses were conducted on all the datasets in RaxML [23], starting with a random tree selection and estimated models. The bootstrap support values were calculated using the autoMRE bootstopping criterion in RaxML [24].

Results
In total, more than 2,223,296 specimens of Culicoides were collected at six sites in 532 collections (mean (x) = 4179) ( Table 1). A total of 25 different Culicoides species were morphologically identified from 418 individuals (Table S1). From this, Culicoides imicola was the predominant species and present at all six sites. Seasonal trends could be observed with an increase in Culicoides numbers during summer and autumn months (January to May) with peak activity in February and a decrease in numbers during winter months (June to August), although never zero. In total, 66/331 pools tested positive for one or more arbovirus, resulting in a MIR of 0.4 (95% confidence interval (CI) 0.0-0.5) ( Table 2). The highest MIR was recorded in Lapalala Wilderness (MIR = 0.8, 95% CI 0.5-1.2) and Marakele National Park (MIR = 0.5, 95% CI 0.2-0.6), both in Limpopo province ( Table 2).
Orthobunyaviruses ( Table 2).   The phylogenetic analysis of a 503 nt sequence of the S segment ( Figure 2) as well as the BLAST results against reference sequences confirmed the results obtained from the 152 nt sequence ( Figure 1) for two samples as SHUV (KYA299_16, and MN055_16) and one sample as SBV (MAR538_16) (Figure 2). The phylogenetic analysis and BLAST of a 503 nt sequence of KYA233_16 identified the sample as SHUV, with this sample lacking a 152 nt sequence. For two samples, KYA358_17 and LAP066_13, conflicting phylogenetic-and p-distance results were seen. KAY358_17 has a 100% nucleotide identity to MAR538_16 (152 nt sequence of the S segment) that clustered within clade B, although it collapsed, probably due to incomplete sampling or lack of reference sequences. In contrast, based on the phylogenetic analysis ( Figure 2) and BLAST results (Table S2) of the 503 nt sequence of the S segment, KAY358_17 groups with SHUV, and therefore was confidently identified as such. Sample LAP066_13 was recovered as sister to SHUV, PEAV and AINV, and MAR278_14 (sister to LAP066_13 + (SHUV + PEAV + AINV)), and therefore could not confidently be identified based on the 152 nt sequence of the S segment ( Figure 1). LAP066_13 also clustered separately (bs = 100) to SHUV, PEAV, and AINV on the 503 nt sequence of the S segment ( Figure 2) re-affirming aforementioned grouping (Figure 1), and therefore remains unknown.  The phylogenetic analysis of a 503 nt sequence of the S segment ( Figure 2) as well as the BLAST results against reference sequences confirmed the results obtained from the 152 nt sequence ( Figure 1) for two samples as SHUV (KYA299_16, and MN055_16) and one sample as SBV (MAR538_16) (Figure 2). The phylogenetic analysis and BLAST of a 503 nt sequence of KYA233_16 identified the sample as SHUV, with this sample lacking a 152 nt sequence. For two samples, KYA358_17 and LAP066_13, conflicting phylogeneticand p-distance results were seen. KAY358_17 has a 100% nucleotide identity to MAR538_16 (152 nt sequence of the S segment) that clustered within clade B, although it collapsed, probably due to incomplete sampling or lack of reference sequences. In contrast, based on the phylogenetic analysis ( Figure 2) and BLAST results (Table S2) of the 503 nt sequence of the S segment, KAY358_17 groups with SHUV, and therefore was confidently identified as such. Sample LAP066_13 was recovered as sister to SHUV, PEAV and AINV, and MAR278_14 (sister to LAP066_13 + (SHUV + PEAV + AINV)), and therefore could not confidently be identified based on the 152 nt sequence of the S segment ( Figure 1). LAP066_13 also clustered separately (bs = 100) to SHUV, PEAV, and AINV on the 503 nt sequence of the S segment ( Figure 2       the PCR results of MN009_15, MN23_15, and MN47_16, although the Rhabdoviridae positives remain unclassified. Samples MN009_15, MN047_16, and MN023_15 formed a well-supported sister group (bs = 100) to Riverside-(RISV) and Tongilchon virus (TCHV) (bs = 99) which, in turn, was well supported as a yet unclassified monophyletic group (bs = 99). The pdistance analysis (not shown) reveals the highest nucleotide identities (86.6%) of samples MN009_15, MN23_15, and MN47_16 to RISV (KU248086, KU248087, and KU248085, respectively).

Discussion
This study used field surveillance and molecular screening of Culicoides to detect viral RNA previously implicated in outbreaks in animals in South Africa and establish a possible association of Culicoides as vectors for the identified viruses. Midge pools, collected from 2012 to 2017, tested positive for viruses from the Alphavirus and Orthobunyavirus genera and family Rhabdoviridae as well as EEV. No Flavivirus or Phlebovirus positives were reported. The highest positivity of viruses was detected in February followed by January and March. This is typical seasonality for arboviruses in South Africa as the arbovirus season is late summer to autumn (February to May). Pools collected during winter months were small or not collected at all, resulting in no positives during July. Overwintering of viruses has, however, been demonstrated for BTV [26] and AHSV [27] in South Africa, even at low insect activity.
Orthobunyaviruses, specifically, SHUV had the highest prevalence across the sites with the highest minimum infection rate in Marakele National Park, while no orthbunyaviruses were detected in KNP in our collection. Phylogenetic analysis based on a 152 nt sequence of the S segment clustered the newly sequenced SHUV strains with known SHUV sequences from horses [28] possibly indicating transmission from midges. In addition, the strains that clearly grouped with SHUV, 15 other orthobunyviruses were

Discussion
This study used field surveillance and molecular screening of Culicoides to detect viral RNA previously implicated in outbreaks in animals in South Africa and establish a possible association of Culicoides as vectors for the identified viruses. Midge pools, collected from 2012 to 2017, tested positive for viruses from the Alphavirus and Orthobunyavirus genera and family Rhabdoviridae as well as EEV. No Flavivirus or Phlebovirus positives were reported. The highest positivity of viruses was detected in February followed by January and March. This is typical seasonality for arboviruses in South Africa as the arbovirus season is late summer to autumn (February to May). Pools collected during winter months were small or not collected at all, resulting in no positives during July. Overwintering of viruses has, however, been demonstrated for BTV [26] and AHSV [27] in South Africa, even at low insect activity.
Orthobunyaviruses, specifically, SHUV had the highest prevalence across the sites with the highest minimum infection rate in Marakele National Park, while no orthbunyaviruses were detected in KNP in our collection. Phylogenetic analysis based on a 152 nt sequence of the S segment clustered the newly sequenced SHUV strains with known SHUV sequences from horses [28] possibly indicating transmission from midges. In addition, the strains that clearly grouped with SHUV, 15 other orthobunyviruses were identified as SHAV, SABV, SANV, SBV, SATV, and INGV, although the phylogenetic clustering of the Simbu serogroup could not confidently be resolved to the species level due to short-sequenced regions. The Shamonda, Sabo, and Sango viruses have previously been identified in parts of Africa, the Middle East, and Asia as a livestock disease responsible for fetal congenital deformities [29,30] and little is known on the zoonotic potential of these viruses. The Schmallenberg virus has been responsible for large outbreaks of abortion storms and fetal malformations in cattle resulting in large economic losses in Europe [3]. The Schmallenberg virus RNA was also recently detected in Culicoides pools as well as in malformed livestock in Israel [31]. Antibodies to SBV have been reported in ruminants of neighboring countries of Namibia [32] and Mozambique [33], indicating the virus may be circulating in southern Africa; however, the close relationship with other members of this group, such as is apparent from our study, may also suggest possible cross reactivity of antibodies in the Simbu serogroup [34]. The Ingwavuma virus was the only virus not from clade B in the Simbu serogroup. The Ingwavuma virus has been isolated in South Africa from a bird and two pools of Culex mosquitoes [35]. Since then, antibodies to INGV have been detected in water buffaloes, goats, sheep, and pigs in part of Indonesia and Thailand where pigs are natural hosts of INGV [36,37].
Additional PCR assays were successful in obtaining longer sequences for six orthobunyavirus positive samples and could confidently identify five viruses based on the S segment. One orthobunyavirus positive remained unidentified (LAP066_13). This may suggest a new or uncharacterized virus and requires further investigation including additional reference sequences. The conflicting data for sample KYA358_17 could be because the S segment is highly conserved, where the smaller sequence is very similar to various members of the broader Simbu complex with the larger, 503 nt sequence revealing higher identity to SHUV. The variability of the S segment between strains could also be due to reassortment as the multi-segmented nature of orthobunyaviruses allows for genetic exchange to increase the fitness, spread, and pathogenicity of the virus or serogroup [38]. The inability to generate larger sequences for all samples or sequence to the M and L gene segment could be due to the low level of RNA present in the vector pool or due to genetic differences in the primer area that was designed based on SHUV sequences isolated from horses.
During the six-year surveillance, MIDV was only detected twice in midge pools from KNP, despite detecting animal cases [13]. The low positivity of MIDV in Culicoides pools during the surveillance period suggests that the virus is present in the area and that the midges were able to ingest the virus but that they might play a less important role as a vector. Phylogenetic analyses cluster the MIDV detected in the midges with MIDV strains found in horses [18]. This could indicate circulation of various genotypes in that area or structural changes in the gene caused by the midge ingestion [39]. No flaviviruses, such as WNV, were detected in this study. Previous studies have reported on WNV detection in Culicoides [7], although their role in the epidemiology of WNV is still unclear.
The high prevalence of EEV is not surprising, as Culicoides midges, specifically C. imicola and Culicoides bolitinos, are the primary vectors for EEV. Cases of febrile, neurological, and respiratory disease in horses infected with EEV were also frequently detected across South Africa during the same period [14]. Equine encephalosis virus detected in midges formed two distinct clusters with the majority grouping with the Kyalami and Northrand serotypes (clade B) and not the Bryanston and Potchefstroom serotypes, as previously reported [40]. No EEV was detected in midge pools from Mnisi or KNP with the highest EEV MIR reported in Lapalala Wilderness. Previous studies have indicated antibodies to EEV in zebras (Zebra burchelli) from KNP [41], although less is published on the epidemiology of EEV and the circulation of this predominantly horse pathogen within wildlife areas. Further investigations into the role of wildlife reservoirs are needed.
The four viruses from the Rhabdoviridae family were accidentally discovered after the flavivirus genus PCR produced bands of the right size on gel electrophoreses. Sanger sequence and p-distance analyses revealed uncharacterized viruses from the Rhabdoviridae family. The phylogenetic analysis of the L segment clusters these viruses with unclassified viruses (RISV and TCHV), although further sequence data is needed to characterize the virus species. The Riverside virus was first isolated in Hungary from an unknown Aedes species, while TCHV was isolated in South Korea from Culex bitaeniorhynchus [42,43]. Grouping of MN009_15, MN23_15, and MN47_16 with mosquito-associated viruses rather than with Culicoides-associated viruses such as Sweetwater Branch, Bivens Arm, and Tibrogargan virus indicates possible incidental ingestion. Another reason could be that both RISV and TCHV were only recently discovered and could be associated with Culicoides and warrant further investigation. The Culicoides-transmitted Rhabdoviridae viruses has been implicated as disease-causing agents in Australia, the Americas, Asia, and Africa [44,45], as well as bovine ephemeral fever virus in South Africa [39].
Unfortunately, virus isolation was unsuccessful on mammalian cells (Vero and BHK) and will, in the future, be attempted on insect cells such as C6/36 or KC cell. Further attempts at next generation sequencing and virus isolation will be made to resolve the unplaced and unidentified viruses. A limiting factor could be large pool sizes resulting in decreasing detection rate, while future studies could benefit by decreasing pool sizes to <100.
The six vector surveillance sites were selected based on previous outbreaks of neurological disease in animals due to arbovirus infections as well as the strong humans-animal interface at these sites [10][11][12][13][14][15]. The highest overall MIR was reported in Lapalala Wilderness. Studies have indicated a MIR ≥1 result in outbreaks of arboviruses such as WNV and eastern equine encephalitis virus in the Americas [46,47]. Although the Orthobunyavirus MIR in Lapalala Wilderness was <1, indicating a low probability of an outbreak, a SHUV outbreak in wildlife occurred from the end of 2017 to mid-2018 [12]; only a small portion of the Culicoides population was sampled and more extensive surveillance could have accurately predicted the outbreak. Marakele National Park had the highest virus diversity with various orthobunyaviruses and EEV detected. In Mpumalanga, there is a strong human-livestock-wildlife interface at the collection sites. Spillovers of disease from wildlife to domestic animals have been reported in these areas [48]. In Boschkop and Kyalami, numerous outbreaks of WNV [10], SHUV [27], MIDV [18], and EEV [14] in horses have been reported, which consist of agricultural small holdings rather than wildlife areas. Morphological identification of Culicoides species from subsamples collected at the various surveillance sites demonstrates the high abundance of C. imicola at all sites, supporting previous findings [49,50]. The abundance of C. imicola and Culicoides schultzei in areas with high infection rates could provide valuable insight into these species as vectors for emerging arboviruses in South Africa.
Taking everything into consideration, detection of viral RNA by PCR does not prove vector status for Culicoides midges, because oral susceptibility and transmission as well as the association of the infected arthropod with the vertebrate population in which the infection is occurring need to be demonstrated [1]. Screening of field collected Culicoides, however, represent a non-invasive method to determine virus presence in an area without using invasive methods such as to immobilize and bleed animals or involved human patients. Previous studies have also implicated field collected Culicoides in the spread of SHUV [8] and successful virus dissemination has been demonstrated for SHUV in the European Culicoides nubeculosus and the North American Culicoides sonorensis [51].
This study identified viruses of medical and veterinary importance from RNA isolated from Culicoides collected in wildlife and urban areas in South Africa. Even though vector competence was not shown, the data suggest that Culicoides may play a role in disease outbreaks warranting further investigations into their potential as vectors for the identified viruses and their capacity to spread these pathogens to new regions.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/ 10.3390/v13101978/s1: Table S1, Culicoides species morphologically identified from subsamples collected at six vector surveillance sites in South Africa from 2012-2017 using 220 V Onderstpoort Light traps. Adapted from [3], Table S2, Collection information pertaining to the Culicoides pools screened in this study as well as phylogenetic and BLAST results of the viruses (family and species) detected in midge pools. GenBank accession number are indicated where applicable.
Author Contributions: J.S., conceptualization, methodology, validation, formal analysis, investigation, writing-original draft; G.J.V., investigation, data curation, writing-review and editing; M.V., conceptualization, methodology, resources, writing-review and editing, supervision, funding acquisition. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement:
The data that support the findings of this study are openly available from GenBank, accession numbers MN270991-MN271022.

Acknowledgments:
We would like to thank all current and previous staff and students from the CVZ, especially T Johnson and T Nelufule, that participated in fieldwork and collections of midges and made this study possible over the years, and L Snyman for assisting with phylogenetic analyses and M Ridden and K Labuschagne for assisting with midge identification.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
For one sample, MN055-16, a partial (395 nt) M segment, was obtained (accession number in Table S2). The SHUV M segment primers were designed based on the genomic position of reference SHUV sequences from GenBank. The conserved region between amino acid positions 2306 and 3266 of the M segment was used to identify first round (forward 23065 -AGCATGYGCYAGGGATACAC-3 2326 and reverse 32665 -CTGRCTGATTAATGGCTGC-3 3247) and nested (forward 23415 -TTGAACCAGCAACA CGAGGA-3 2361 and reverse 27595 -ATGTGGGTGCAATRGARGGC-3 2739) primers using the NCBI Primer designing tool (National Center for Biotechnology Information, www.ncbi.nlm.nih.gov (accessed on 30 March 2018)) resulting in a nested PCR product of 418 bp. The cDNA was synthesised using a RevertAid First Strand cDNA Synthesis Kit (Thermo Fisher Scientific, Waltham, MA, USA) and random priming (Sigma-Aldrich, St. Louis, Mo, USA) according to the manufacturer's instructions, incubated for 5 min at 25 • C, then 60 min at 45 • C, and finally for 5 min at 70 • C. The resulting cDNA was used in a first round PCR using Platinum™ Taq DNA polymerase (Thermo Fisher Scientific). Each 25 µL contained 10 pmol forward and reverse primers, Platinum™ Taq DNA polymerase buffer, 50 mM MgCl 2 (Thermo Fisher Scientific), 10 mM dNTPs (Thermo Fisher Scientific), and Platinum™ Taq DNA Polymerase enzyme. Cycling conditions were as follow: 94 • C for 2 min, followed by 35 cycles of 94 • C for 30 s, 55 • C for 1 min, and 72 • C for 1 min, and a final extension of 72 • C for 7 min. A nested PCR was, then, also performed using Platinum™ Taq DNA polymerase (Thermo Fisher Scientific). Each 25 µL contained 10 pmol nested forward and nested reverse primers, Platinum™ Taq DNA polymerase buffer, 50 mM MgCl 2 (Thermo Fisher Scientific), 10 mM dNTPs (Thermo Fisher Scientific), and Platinum™ Taq DNA polymerase enzyme. The cycling conditions were the same as the first round PCR.