Ecology

Characteristics of the single agrometeorological disaster scenariosSAD-f occurred in 49 out of 57 years at different spatial scales, with a maximum IOC of 0.519 in 2013; SAD-d occurred in 33 years with a maximum IOC value of 0.808 in 1995; SAD-s occurred in 5 years, with a maximum IOC value of 0.115 in 1977 (Fig. 4). SAD-d showed a declining trend over the past 57 years, but the SAD-d frequency was higher than SAD-f and SAD-d. Since the mid-1980s, the frequency of SAD-f has increased, while the frequency and scale of SAD-s were relatively small.Figure 4IOC change curve for single agrometeorological disasters (SAD) in different scenarios.Full size imageA large-scale grade SAD-f event occurred in 2013 in Liaoning Province, and regional SAD-f occurred in 17 years. Five years showed a large-scale SAD-d and 11 years had regional SAD-d. There were no large-scale and regional years for SAD-s (Table 5).Table 5 Occurrence years of large-scale and regional single agrometeorological disaster (SAD) in Liaoning Province.Full size tableThe occurrence of SAD-f was recorded at 15 sites with a frequency greater than 20%, 13 sites with frequency in the 10%–20% range, and 24 sites with a low frequency (P ≤ 10%) in Liaoning Province from 1961 to 2017 (Fig. 5). SAD-d occurred at 12 sites with a frequency higher than 20%, 22 sites with frequency between 10% and 20%, and 18 sites with a low frequency (P ≤ 10%). In the three scenarios, the occurrence frequency and distribution of SAD-f was the highest and SAD-s was the lowest.Figure 5Frequency of single agrometeorological disasters (SAD) in different scenarios (a SAD-f, b SAD-d, c SAD-s). Maps generated in ArcGIS 9.3.Full size imageComparison of the characteristics of single agrometeorological disasters and combined agrometeorological disastersThe maximum IOC of SAD was 0.808 in 1995 and the mean value was 0.294 for the 57 years of the study; the maximum IOC of CAD was 0.654 in 1987, and the mean value was 0.180 over the past 57 years; SAD and CAD occurred in all 57 years (Fig. 6). Both SAD and CAD showed declining trends from 1961 to 2017. The IOC was lower for CAD than for SAD for 42 years and higher than SAD for 14 years.Figure 6Change in the IOC for agrometeorological disasters in rice crops.Full size imageThis paper analysed the mean IOC of SAD and CAD over six decades and found that the interdecadal mean value of the IOC in CAD was lower than that of SAD over five of the periods, but the IOC of SAD was lower than that of CAD in 1971–1980 (Fig. 7). The IOC of SAD showed a decreasing trend from the 1970s to the 2010s but showed an increasing trend after 2011. The IOC of CAD showed a decreasing trend from the 1970s to the 2000s, but showed an increasing trend after 2001 (Fig. 7).Figure 7Interdecadal mean value of IOC for agrometeorological disasters in rice crops.Full size imageThere was one site (Fushun) with a SAD frequency of more than 50% in Liaoning Province from 1961 to 2017, 42 sites with a frequency between 20% and 50%, and nine sites with a frequency lower than 20%. There were four sites (Xinfeng, Jianping, Xinbin and Caohekou) with a CAD frequency higher than 50%, 13 sites with a frequency in the range 20%–50%, and 35 sites with a frequency lower than 20% (Fig. 8). The frequency and range of CAD were less than those of SAD.Figure 8Frequency of agrometeorological disaster in rice crops (a combined agrometeorological disaster (CAD), b single agrometeorological disaster (SAD)). Maps generated in ArcGIS 9.3.Full size imageThere has been little research into the temporal or spatial distribution of CAD for rice and its occurrence characteristics: most research has been on SAD. For example, studies have examined the characteristics of SCD, DCD, FD for rice in northeast China26, 27, 29, and the risk of multiple disasters for rice in northeast China30, 31. Han et al.31 analysed the risk of disaster using the reduction rate of rice yield in Liaoning Province from 1980 to 2011, and found that the high-risk areas were distributed in the west and northeast of Liaoning Province; higher rates of yield reduction in lean years were mainly found in western Liaoning and its surrounding areas. In this study, a higher frequency of CAD was mainly distributed in the northwest of Liaoning Province, while that of SAD occurred in the northeast of Liaoning Province. The median frequency of CAD occurred in the northwest and northeast of Liaoning Province, while that of SAD covered most areas in Liaoning Province. The range of medium and higher frequency occurrence in CAD was consistent with the distribution of high-risk and high yield reduction areas in the study of Han et al.31. Therefore, it can be speculated that the CAD scenarios might magnify the effect of each single disaster, and, therefore, CAD would more easily lead to a higher reduction in the rice yield.Comparison of the occurrence of single agrometeorological disasters and combined agrometeorological disastersDuring the rice growing season in Liaoning Province, there were three scenarios of SAD and six scenarios of CAD. Compared with SAD, CAD had more scenarios and more complex processes, and its effect on rice was more difficult to evaluate. In SAD, the occurrence frequency and distribution of SAD-f and SAD-d were both high, when FD and DCD occurred alone in only one rice growth stage. In CAD, the occurrence frequency and distribution of TD-1, when FD and DCD occurred simultaneously, was the highest in the six scenarios. A single or combined occurrence of FD and DCD was most common disaster for rice in Liaoning Province. The occurrence frequency and distribution of OD-1 were both smaller than that of SAD-f, indicating that the occurrence was lower when FD happened at both the seedling and milk stages. SAD-s and OD-2 had the lowest frequency and range in all scenarios, indicating that DSD rarely appeared in SAD and CAD. The occurrence of SCD was not major disaster in the growth and development of rice in Liaoning Province, but the occurrence of DCD or FD, or both, was.In this study, the occurrence frequency and range of SAD and CAD for rice showed declining trends in most sites over the past 57 years, which was consistent with the results of other studies. Studies on rice DCD and SCD concluded that cold damage events of rice in most areas of northeast China showed decreasing trends26, 27. Because of events such as climate warming, earlier warming in spring, delaying first frost dates and fewer low temperature days in summer, the trend of disasters was lower in rice planting areas30. However, although rice disasters showed a decreasing trend, local disasters may increase because of the frequent occurrence of climate anomalies. SAD-f and OD-1 scenarios in this study showed no significant decreasing trend, and even a partial increasing trend. Jiang et al.29 believed that the possibility of frequent SCD in north-east China was still high. According to Xi et al.32, cold periods would still occur in the growing season of rice in northeast China. Hu et al.33 concluded that the increase of SCD in northeast China was mainly because of the increase of climate variability, and most of the sites with increases were located in areas with decreasing temperature or no obvious trend of temperature increase.Rice is a higher temperature-loving crop, which is mainly restricted by temperature conditions during its growing season. Liaoning Province is in the south of the rice planting area of the colder regions in China. Because of the relatively low latitude, heat conditions during the rice growing season were better than those in Jilin and Heilongjiang to the north of Liaoning Province. The climatic risk of cold damage in the rice growing season was lower than other regions in northeast China34. The occurrence of CAD was generally caused by low temperatures, which were the dominant factor. When two or more disasters occur together, there is a coupling or amplifying effect on rice growth compared with a single disaster.A comparison of the rice yield reduction rates in the years when CAD or SAD occurred in more than 50% of stations in Liaoning Province revealed that the former happened in 5 years, 1969, 1974, 1976, 1980 and 1987, whereas the latter happened in 7 years, 1972, 1973, 1985, 1986, 1990, 1995 and 2013. When CAD was the major occurrence, the average yield reduction rate in the five years was 10.6%. The yield reduction rate in 1969 was 34.6%, which was the highest in the past 57 years. When SAD was the major occurrence, the average yield reduction rate in the seven years was 9.8%. The average yield reduction rate in the years when CAD dominated was greater than in the years when SAD dominated. Therefore, it can be speculated that CAD has a greater effect on rice growing than any single disaster within CAD. However, it is difficult to quantify the effect on rice yield of CAD, and further controlled field experiments should be conducted to verify these. It is difficult to control field experiments that are limited by conditions and facilities.Comparison of the occurrence of agrometeorological disasters in years having rice yield reductionsOn the basis of the rice yield reduction rate in calculations Liaoning Province from 1961 to 2017, a total of 10 years (Table 6) were screened. Six years had large-scale disasters (including SAD and CAD) and four years had regional disasters. In 1969, which showed the highest yield reduction rate (34.6%), 30 sites had TD-1 disasters and the other 22 sites had SAD-f disasters. In 1972, the second highest reduction year (29.1%), 11 sites had MD-1 disasters, i.e. three kinds of disasters occurred, seven sites had TD-1 disasters, one site had a TD-2 disaster, 31 sites had SAD-d disasters, one site had a SAD-f disaster, and only one station had no disaster. The TD-1 disaster, i.e. delayed cold damage and frost injury, was the most frequent CAD over the years, and SAD-d, i.e., delayed cold damage, was the most frequent SAD. The occurrence of single and combined agrometeorological disasters in different regions strongly affected the rice yield. Generally, the larger the disaster range, the higher the yield reduction. However, some years were not completely consistent with this conclusion. The yield reduction rate was also related to the type, severity, occurrence period and geographical location of the disasters.Table 6 Comparison of agrometeorological disasters in years having greater than 10% rice yield reduction rates in Liaoning Province.Full size tableIn every year from 1961 to 2017, CAD or SAD occurred in Liaoning Province, and the rice yields declined in 23 of the 57 years owing to meteorological disasters (Fig. 9). Although meteorological disasters occurred in the other 34 years, there was no reduction in rice production, which may be related to the gradient of the disaster or the spatial distribution of the rice planting areas. The rice yield reduction rates in 1969 and 1976 were 34.6% and 15.6%, respectively. In these two years, CAD occurred at 30 stations and SAD occurred at 22 stations, and TD-1 was the main type of CAD, whereas SAD-d was the main type of SAD. Using statistical data, on the rice planting area of each city in Liaoning Province, the provincial area can be divided into four regions. The first region was Shenyang City, which has the largest rice planting area, accounting for 20%–25% of the total rice planting area; the second region was Panjin City, accounting for 15%–20% of the total rice area; the third region encompassed Tieling and other six cities, accounting for nearly 50% of the total rice area, with each city representing 5%–10%; and the fourth region encompassed Jinzhou and five other cities, accounting for 10%–15% of the total, with each city representing 0–5%. As shown in Fig. 10a,b, TD-1 occurred in the first region in both 1969 and 1976 and in the second region in 1969. SAD-d occurred in the second region in 1976. In the third region, TD-1 occurred at more stations of 1969 than in 1976. The rice area in the first three regions accounted for nearly 80% of the total rice area, and CAD occurred more often than SAD in these regions. Thus, there was a greater yield reduction rate in 1969 than in 1976.Figure 9The IOC change curve of all agrometeorological disasters and the rice yield reduction rate from 1961 to 2017.Full size imageFigure 10Distributions of the types of agrometeorological disasters and the percentages of rice planting areas in different regions of Liaoning Province in 1969 and 1976 (a: 1969; b: 1976). Maps generated in ArcGIS 9.3.Full size imageThe occurrence characteristics of single disasters or the risk of yield reduction were analysed in previous studies, but the quantitative effect on rice production was rarely evaluated. Ji et al.26 reported that the delayed cold damage in 1961, 1962, 1969, 1972, 1976, 1989 and 1995 was so severe that there was a large reduction in rice production. In our paper, we examined the occurrence of not just one disaster, i.e. delayed cold damage, over time, but also other types of disasters including SAD and CAD. For example, in 1972 and 1976, the disaster scenario affecting the largest number of stations was TD-1, i.e., both delayed cold damage and frost damage occurred in the growing season of rice. In 1961, the most widespread damage came from a single disaster—frost damage. According to the records35, Liaoning Province experienced frost damage in 1961, 1962, 1969, 1972, 1976 and 1995, and the rice yield was seriously reduced. Most regions of Liaoning Province experienced both delayed cold damage and frost damage in 1976 and 1995. There was a low temperature during the critical period of rice growth (mid-July to mid-August) in 1995. In 1985, the growing season in most areas was characterized by unusually persistent low temperature and little sunshine. These statistics were basically consistent with the conclusion of this study. In the process of rice production, a variety of disasters occurred caused by low temperature, such as delayed cold damage, frost damage and sterile cold damage. More

Fish and experimental protocolThree-hundred-fifty market-size (755–3865 g) common carp Cyprinus carpio were obtained from six sources at various times of year to for effects of variation in rearing conditions. The weight of collected carp corresponded to the weight of carp normally delivered to the market. Fish were obtained from the Faculty of Fisheries and Protection of Waters of the University of South Bohemia in Ceske Budejovice (FFPW USB), Vodnany and the fisheries Chlumec nad Cidlinou, Blatna, Hodonin, Klatovy, Lnare, and Tabor. Ten fish were collected from each fishery at the spring (March/April), summer (June/July), and autumn harvests (October/November) in 2018 and 2019. Carp were transported live to the laboratory of the FFPW, killed by a blow to the head, weighed, measured, and filleted. Two fillets, one with skin removed, from each fish were individually vacuum packed, immediately frozen, and stored at − 32 °C until chemical analysis.Ethics approvalAll the methods used in the present study followed relevant guidelines and regulations. Also, the competent authority (Ethical Committee for the Protection of Animals in Research of the University of South Bohemia, FFPW Vodnany) approved the fish sampling and protocols of the present study and reporting herein follows the recommendations in the ARRIVE guidelines.Chemical analysisSeven-hundred carp fillets were analysed for basic nutritional composition, dry matter, protein, fat, and ash. All samples were homogenized by grinding before analysis.The determination of dry matter followed ISO 1442:1997 Meat and meat products—Determination of moisture content (Reference method)26. The homogenized samples were dried with sand to constant weight at 103 ± 2 °C in a laboratory oven (Memmert UE 500, Memmert GmbH + Co. KG, Germany).The determination of ash was based on the standard ISO 936:1998 Meat and meat products—Determination of total ash27. The homogenized samples were burned in a muffle furnace (Nabertherm A11/HR, Nabertherm GmbH, Germany) at 550 ± 25 °C to a grey-white colour.The determination of total fat was based on the standard ISO 1443:1973 Meat and meat products—Determination of total fat content28. The homogenized samples were hydrolysed by hydrochloric acid, and fat was extracted by light petroleum in SOXTEC 2050 (FOSS Headquarters, Denmark).The determination of nitrogen used the Kjeldahl method based on the standard method ISO 937:1978 Meat and meat products—Determination of nitrogen content (Reference method)29. The homogenized samples were digested by sulphuric acid and a catalyser in a KjelROC Digestor 20 (OPSIS AB, Sweden) digestion unit at 420 ± 10 °C. Organically bound nitrogen was measured on the KJELTEC 8400 with KJELTEC sampler 8420 (FOSS Headquarters, Denmark). Calculation of protein content from nitrogen used the conversion factor for meat of 6.25.All analysis of dry matter, ash, and total fat were performed in duplicate and analysis of nitrogen (protein) was performed in triplicate for each sample.Calculation of fat-free nitrogen (Nff) in g/100 g used the formula24:$$ N_{ff} = frac{{100 times N { }}}{{100 – F { }}}. $$This formula was applied to nitrogen (N) and fat (F) content for all samples, providing a fat-free nitrogen value for each sample.Fish meat content calculated based on nitrogen factor Nf (total fillet) in g/100 g used the formula9:$$ Fish ;content_{Nf} = frac{N times 100}{{N_{f} }}. $$Fish meat content calculated based on fat-free nitrogen factor (Nff) and DCC (defatted carp content) in g/100 g used formulas11:$$ Fishc; content_{Nff} = DCC + F, $$$$ DCC = frac{N times 100}{{N_{ff} }}. $$Statistical analysisKolmogorov–Smirnov and Bartlett’s tests were applied to assess normal distribution data and the homoscedasticity of variance, respectively. A two-way ANOVA and Tukey’s test was conducted to analyse effects of season, weight, fishery, and difference between fillets with and without skin. The significance level was set at p  More

Uganda has been exceptionally successful in scaling-up coverage of LLINs. Following the mass distribution campaigns to deliver free LLINs in 2013–14 and 2017–18, 90 and 83% of households, respectively reported ownership of at least one LLIN7,14. However, despite this success, the burden of malaria remains high in much of the country. Uganda had the third highest number of malaria cases reported in 2019, with reported case incidence increasing since 20142. If Uganda is to achieve the goals established by the World Health Organization’s Global Technical Strategy for malaria including reducing malaria case incidence by at least 90% by 2030 as compared with 201515, additional malaria control measures will be needed. This report highlights the critical role of IRS in substantially reducing the burden of malaria in areas where transmission remains high despite deployment of LLINs. Withdrawing IRS after 5 years of sustained use in three districts in northern Uganda was associated with a more than fivefold increase in malaria cases within 10 months. Restarting IRS with a single round in nine districts in Northern Uganda ~3 years after IRS had been stopped was associated with a transient but important (more than a fivefold) decrease in malaria cases within 8–12 months, returning to pre-IRS levels after 23 months. Initiating and sustaining IRS in five districts in Eastern Uganda was associated with a gradual reduction in malaria cases reaching almost a sevenfold reduction after 4–5 years.Robust evidence supports the widespread use of LLINs for malaria control. In a systematic review of clinical trials conducted between 1987 and 2001, insecticide treated nets reduced all cause child mortality by 17% and the incidence of uncomplicated P. falciparum malaria by almost half16. However, there is concern that the effectiveness of LLINs may be diminishing due to widespread resistance to pyrethroids which until recently were the only class of insecticides approved for LLINs. Similar to many other African countries, high-level resistance to pyrethroids among the principle Anopheles vectors has been reported recently throughout Uganda17,18,19. In addition, behavioral changes in vector biting activity following the introduction of LLINs have been reported which could present new challenges for malaria control20,21,22. Finally, the effectiveness of LLINs may be further compromised by poor adherence and waning coverage in the setting of free distribution campaigns done intermittently. In Uganda, less than 18% of households reported adequate coverage (defined as at least one LLIN per two residents) 3 years after the 2013–14 distribution campaign23 and adequate coverage decreased from 71% to 51% between 6 and 18 months following the 2017–18 distribution campaign24. Although the World Health Organization recommends mass distribution campaigns every 3 years, mounting evidence suggests that LLINs should be distributed more frequently to sustain high coverage25,26,27,28,29,30,31.Given concerns about the current effectiveness of pyrethroid-based LLINs and the persistently high burden of malaria despite aggressive scale-up of LLINs in countries like Uganda, additional malaria control measures are needed. IRS is an attractive option. Historically, IRS programs were used to dramatically reduce and even eliminate malaria in many parts of the world. Thus, while there is some evidence for the impact of IRS in the absence of LLINs32, it is surprising that the evidence base from contemporary controlled trials on the impact of adding IRS to LLINs for vector control is limited. A recent systematic review of cluster randomized controlled trials conducted in sub-Saharan Africa since 2008, reported that adding IRS using a “pyrethroid-like” insecticide to LLINs did not provide any benefits, while adding IRS with a “non-pyrethroid-like” insecticide produced mixed results5. Among the four trials comparing IRS plus LLINs with LLINs alone, three evaluated IRS with a carbamate (bendiocarb) and one evaluated a long-lasting organophosphate, pirimiphos-methyl (Actellic 300CS®)33,34,35,36. Only two trials (both using bendiocarb) assessed malaria incidence; one from Sudan found a 35% reduction when adding IRS to LLINs34, while another from Benin found no benefit of adding IRS33. All four trials assessed parasite prevalence, with an overall non-significant trend towards a lower prevalence when adding IRS to LLINs (RR = 0.67, 95% CI 0.35–1.28)5. However, when the analyses were restricted to include only the two studies with LLIN usage over 50%, adding IRS reduced parasite prevalence by over 50% (RR = 0.47, 95% CI 0.33–0.67)5. Of note, none of the trials that evaluated the impact of adding IRS with a “non-pyrethroid-like” insecticide assessed outcomes beyond 2 years. More recently, a number of observational studies have reported benefits of using IRS with pirimiphos-methyl (Actellic 300CS®). In the Mopti Region of Mali, delivery of a single round of IRS with Actellic 300CS® was associated with a 42% decrease in the peak incidence of laboratory-confirmed malaria cases reported at public health facilities37. In the Koulikoro Region of Mali, villages that received a single round of IRS with Actellic 300CS® combined with LLINs observed a greater than 50% decrease in the incidence of malaria compared to villages that only received LLINs38. In the Northern Region of Ghana, districts that received IRS with Actellic 300CS® reported 26–58% fewer cases of laboratory-confirmed malaria cases reported at public health facilities over a 2-year period, compared to districts that did not receive IRS39. In Northern Zambia, implementation of IRS with Actellic 300CS® targeting only high burden areas over a 3 year period was associated with a 25% decline in parasite prevalence during the rainy season, but no decline during the dry season40. In Western Kenya, the introduction of a single round of IRS with Actellic 300CS® was associated with a 44–65% decrease in district level malaria case counts over a 10 month period compared to pre-IRS levels41. In addition, several recent reports have documented dramatic resurgences of malaria following the withdrawal of IRS with bendiocarb in Benin42, and the withdrawal of IRS with Actellic 300CS® in Mali and Ghana37,39.The results from this study provides additional support for the critical role IRS can play in reducing the burden of malaria in African countries with high LLINs coverage. A strength of the study was its use of a large, rigorously collected dataset. Data were collected over nearly 7 years through an enhanced health facility-based surveillance system covering 14 districts in Uganda where IRS was being withdrawn, re-started, and initiated. This enhanced surveillance system facilitated laboratory testing and provided prospectively collected, individual-level data, allowing for analyses of quantitative changes in laboratory-confirmed cases of malaria over time, controlling for temporal changes in rainfall, seasonal effects, diagnostic practices, and health seeking behavior. Previous work by our group documented a marked decrease in malaria TPRs after 4 years of sustained IRS with bendiocarb in one district of Northern Uganda followed by a rapid resurgence over an 18-month period after IRS was withdrawn11. In this study we expand on these findings by including data from three districts and covering a 31-month period following the withdrawal of IRS. We were able to quantify more than a fivefold increase in malaria cases which was sustained over the 10–31 months following the withdrawal of IRS. This marked resurgence occurred despite the fact the first universal LLIN distribution campaign was timed to occur right after IRS was withdrawn. Given the dramatic nature of the resurgence, the Ugandan government was able to procure funding for a single round of IRS with Actellic 300CS® ~3 years after IRS was withdrawn in 10 districts of Northern Uganda. In this study, we assessed the impact of this single round in nine of these districts. This single round was associated with over a fivefold decrease in malaria cases after 8–12 months, with malaria cases returning to pre-IRS levels after almost 2 years. These data suggest that IRS with longer-acting formulations such as Actellic 300CS® administered every 2 years could be considered as a strategy for mitigating the risk of resurgence following sustained IRS and/or enabling countries to expand coverage when resources are limited, but formal assessment and a cost-effectiveness analyses are needed. This study also evaluated the impact of 5 years of sustained IRS in five districts of Eastern Uganda, starting first with bendiocarb and then switching to Actellic 300CS® after 18 months. Rounds of IRS were initially associated with marked decreases in malaria cases followed by peaks before subsequent rounds until the fourth and fifth years after IRS was initiated when there was a sustained decrease of almost sevenfold compared to pre-IRS level. Given the before-and-after nature of our study design, it is not clear whether the maximum sustained benefits of IRS seen after 4–5 years were due to the cumulative effect of multiple rounds of IRS, the switch from bendiocarb to Actellic 300CS®, improvements in implementation (although campaigns occurred regularly and coverage was universally high across rounds, see Supplementary Table 4), the second universal LLIN distribution campaign which occurred in this area in 2017, and/or other factors.This study had several limitations. First, we used an observational study design, with measures of impact based on comparisons made before-and-after key changes in IRS policy. Although cluster randomized controlled trials are the gold standard study design for estimating the impact of IRS, it could be argued that withholding IRS would be unethical, given what is known about its impact in Uganda. Second, our estimates of impact could have been confounded by secular trends in factors not accounted for in our analyses. However, we feel that our overall conclusions are robust given the large amount of data available from multiple sites over an extended period with multiple complementary objectives providing consistent findings. Third, we could not assess the impact of IRS independent of LLIN use and did not have access to measures of IRS or LLIN coverage from our study populations. It is possible that some of the impacts we observed were from LLIN distributions in combination with IRS campaigns. However, we were able to provide a “real world” assessment of IRS in a setting where LLIN use is strongly supported by repeated universal distribution campaigns that are becoming increasingly common in sub-Saharan Africa. Similarly, we cannot draw conclusions on the impact of different IRS compounds given all sites received the same formulations consecutively. The results from Objective 3 indicate that malaria incidence dropped substantially in the years that districts stopped receiving bendiocarb and began receiving Actellic 300CS®. However, we cannot conclude whether this reduction was a result of this change or rather the cumulative impact of sustained IRS campaigns, as it has been suggested that in very high transmission settings, several years of IRS may be needed to maximize impact on measures of morbidity.43,44 Finally, our study outcome was limited to case counts of laboratory-confirmed malaria captured at health facilities. Thus, we were unable to measure the impact of IRS on other important indicators such as measures of vector distribution, parasite prevalence, or mortality.There is a growing body of evidence that combining LLINs with IRS using “non-pyrethroid-like” insecticides, especially the long acting organophosphate Actellic 300CS®, is highly effective at reducing the burden of malaria in Uganda, and elsewhere in Africa. Despite these encouraging findings, IRS coverage in Africa has been moving in the wrong direction. The proportion of those at risk protected by IRS in Africa peaked at just over 10% in 2010. However, the spread of pyrethroid resistance has led many control programs to switch to more expensive formulations resulting in a 53% decrease in the number of houses sprayed between years of peak coverage and 2015 across 18 countries supported by the US President’s Malaria Initiative45 and an overall reduction in the proportion protected by IRS in Africa to less than 2% in 20192. Given the lack of recent progress in reducing the global burden of malaria coupled with challenges in funding, renewed commitments are needed to address the “high burden to high impact” approach now being advocated by the World Health Organization2. IRS is a widely available tool that could be scaled up, however demands currently exceed the availability of resources. Additional work is needed to optimize the use of IRS, prevent further spread of insecticide resistance, and better evaluate the cost-effectiveness of IRS in the context of other control interventions. More

1.Piertney, S. B. In Current Trends in Wildlife Research (eds Mateo, R. et al.) 201–223 (Springer International Publishing, 2016).2.Jabin, G., Sahajpal, V., Chandra, K. & Thakur, M. In Forensic DNA Typing: Principles, Applications and Advancements (eds Shrivastava, P., et al.) 399–403 (Springer, Singapore, 2020).3.Drechsler, A., Helling, T. & Steinfartz, S. Genetic fingerprinting proves cross-correlated automatic photo-identification of individuals as highly efficient in large capture–mark–recapture studies. Ecol. Evol. 5, 141–151. https://doi.org/10.1002/ece3.1340 (2015).Article 
PubMed 

Google Scholar 
4.Gupta, S. K. In DNA Fingerprinting: Advancements and Future Endeavors (eds Ranjan Dash, H., et al.) 77–87 (Springer, Singapore, 2018).5.Dale, T. D. et al. Enhancement of wildlife disease surveillance using multiplex quantitative PCR: Development of qPCR assays for major pathogens in UK squirrel populations. Eur. J. Wildl. Res. 62, 589–599. https://doi.org/10.1007/s10344-016-1031-z (2016).Article 
PubMed 
PubMed Central 

Google Scholar 
6.Latch, E. K., Heffelfinger, J. R., Fike, J. A. & Rhodes, O. E. Jr. Species-wide phylogeography of North American mule deer (Odocoileus hemionus): Cryptic glacial refugia and postglacial recolonization. Mol. Ecol. 18, 1730–1745. https://doi.org/10.1111/j.1365-294X.2009.04153.x (2009).Article 
PubMed 

Google Scholar 
7.Moscarella, R. A., Aguilera, M. & Escalante, A. A. Phylogeography, population structure, and implications for conservation of white-tailed deer (Odocoileus virginianus) in Venezuela. J. Mammal. 84, 1300–1315. https://doi.org/10.1644/brb-028 (2003).Article 

Google Scholar 
8.Lang, K. R. & Blanchong, J. A. Population genetic structure of white-tailed deer: Understanding risk of chronic wasting disease spread. J. Wildl. Manag. 76, 832–840. https://doi.org/10.1002/jwmg.292 (2012).Article 

Google Scholar 
9.Locher, A., Scribner, K. T., Moore, J. A., Murphy, B. & Kanefsky, J. Influence of landscape features on spatial genetic structure of white-tailed deer in human-altered landscapes. J. Wildl. Manag. 79, 180–194 (2015).Article 

Google Scholar 
10.Green, M. L., Manjerovic, M. B., Mateus-Pinilla, N. & Novakofski, J. Genetic assignment tests reveal dispersal of white-tailed deer: Implications for chronic wasting disease. J. Mammal. 95, 646–654. https://doi.org/10.1644/13-MAMM-A-167 (2014).Article 

Google Scholar 
11.Zachos, F. E. et al. Population viability analysis and genetic diversity of the endangered red deer Cervus elaphus population from Mesola, Italy. Wildl. Biol. 15, 175–186 (2009).Article 

Google Scholar 
12.Villanova, V. L., Hughes, P. T. & Hoffman, E. A. Combining genetic structure and demographic analyses to estimate persistence in endangered key deer (Odocoileus virginianus clavium). Conserv. Genet. 18, 1061–1076. https://doi.org/10.1007/s10592-017-0958-2 (2017).Article 

Google Scholar 
13.Lorenzini, R. DNA forensics and the poaching of wildlife in Italy: A case study. Forensic Sci. Int. 153, 218–221. https://doi.org/10.1016/j.forsciint.2005.04.032 (2005).CAS 
Article 
PubMed 

Google Scholar 
14.Vikas, K. Wildlife DNA evidence: Recognition, collection and preservation. Res. J. Forensic Sci. 3(7), 8–15 (2015).
Google Scholar 
15.Waits, L. P. & Paetkau, D. Noninvasive genetic sampling tools for wildlife biologists: A review of applications and recommendations for accurate data collection. J. Wildl. Manag. 69, 1419–1433 (2005).Article 

Google Scholar 
16.Oyler-McCance, S. J. & Leberg, P. L. In The Wildlife Techniques Manual: Research (ed 7th) 526–546 (John Hopkins University Press, 2012).17.Begley-Miller, D. R., Hipp, A. L., Brown, B. H., Hahn, M. & Rooney, T. P. White-tailed deer are a biotic filter during community assembly, reducing species and phylogenetic diversity. AoB Plants https://doi.org/10.1093/aobpla/plu030 (2014).Article 
PubMed 
PubMed Central 

Google Scholar 
18.Saunders, S. E., Bartelt-Hunt, S. L. & Bartz, J. C. Occurrence, transmission, and zoonotic potential of chronic wasting disease. Emerg. Infect. Dis. 18, 369–376 (2012).Article 

Google Scholar 
19.Miller, W. L., Edson, J., Pietrandrea, P., Miller-Butterworth, C. & Walter, W. D. Identification and evaluation of a core microsatellite panel for use in white-tailed deer (Odocoileus virginianus). BMC Genet. 20, 49. https://doi.org/10.1186/s12863-019-0750-z (2019).Article 
PubMed 
PubMed Central 

Google Scholar 
20.Budd, K., Berkman, L. K., Anderson, M., Koppelman, J. & Eggert, L. S. Genetic structure and recovery of white-tailed deer in Missouri. J. Wildl. Manag. 82, 1598–1607. https://doi.org/10.1002/jwmg.21546 (2018).Article 

Google Scholar 
21.DeYoung, R. W., Demarais, S., Gonzales, R. A., Honeycutt, R. L. & Gee, K. L. Multiple paternity in white-tailed deer (Odocoileus Virginianus) revealed by DNA microsatellites. J. Mammal. 83, 884–892. https://doi.org/10.1644/1545-1542(2002)083%3c0884:mpiwtd%3e2.0.co;2 (2002).Article 

Google Scholar 
22.Poutanen, J., Pusenius, J., Wikström, M. & Brommer, J. E. Estimating population density of the white-tailed deer in Finland using non-invasive genetic sampling and spatial capture-recapture. Ann. Zool. Fenn. 56, 1–16 (2019).Article 

Google Scholar 
23.Brinkman, T. J. et al. Individual identification of Sitka black-tailed deer (Odocoileus hemionus sitkensis) using DNA from fecal pellets. Conserv. Genet. Resour. 2, 115–118. https://doi.org/10.1007/s12686-010-9176-7 (2010).Article 

Google Scholar 
24.de Vargas Wolfgramm, E. et al. Simplified buccal DNA extraction with FTA® Elute Cards. Forensic Sci. Int. Genet. 3, 125–127. https://doi.org/10.1016/j.fsigen.2008.11.008 (2009).CAS 
Article 

Google Scholar 
25.Bunting, S., Burnett, E., Hunter, R. B., Field, R. & Hunter, K. L. Incorporating molecular genetics into remote expedition fieldwork. Trop. Conserv. Sci. 7, 260–271. https://doi.org/10.1177/194008291400700207 (2014).Article 

Google Scholar 
26.McClure, M. C., McKay, S. D., Schnabel, R. D. & Taylor, J. F. Assessment of DNA extracted from FTA cards for use on the Illumina iSelect BeadChip. BMC Res Notes 2, 107. https://doi.org/10.1186/1756-0500-2-107 (2009).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
27.Milne, E. et al. Buccal DNA collection: Comparison of buccal swabs with FTA cards. Cancer Epidemiol. Biomark. Prev. 15, 816–819. https://doi.org/10.1158/1055-9965.epi-05-0753 (2006).CAS 
Article 

Google Scholar 
28.Smith, L. M. & Burgoyne, L. A. Collecting, archiving and processing DNA from wildlife samples using FTA databasing paper. BMC Ecol 4, 4–4. https://doi.org/10.1186/1472-6785-4-4 (2004).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
29.Fryxell, R. T. T. et al. Survey of Borreliae in ticks, canines, and white-tailed deer from Arkansas, U.S.A. Parasit. Vectors 5, 139. https://doi.org/10.1186/1756-3305-5-139 (2012).Article 
PubMed 
PubMed Central 

Google Scholar 
30.Picard-Meyer, E., Barrat, J. & Cliquet, F. Use of filter paper (FTA®) technology for sampling, recovery and molecular characterisation of rabies viruses. J. Virol. Methods 140, 174–182. https://doi.org/10.1016/j.jviromet.2006.11.011 (2007).CAS 
Article 
PubMed 

Google Scholar 
31.Haley, N. J. et al. Antemortem detection of chronic wasting disease prions in nasal brush collections and rectal biopsy specimens from white-tailed deer by real-time quaking-induced conversion. J. Clin. Microbiol. 54, 1108–1116. https://doi.org/10.1128/jcm.02699-15 (2016).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
32.Mas, S., Crescenti, A., Gassó, P., Vidal-Taboada, J. M. & Lafuente, A. DNA cards: Determinants of DNA yield and quality in collecting genetic samples for pharmacogenetic studies. Basic Clin. Pharmacol. Toxicol. 101, 132–137. https://doi.org/10.1111/j.1742-7843.2007.00089.x (2007).CAS 
Article 
PubMed 

Google Scholar 
33.Drew, R. E. et al. Conservation genetics of the fisher (Martes pennanti) based on mitochondrial DNA sequencing. Mol. Ecol. 12, 51–62. https://doi.org/10.1046/j.1365-294X.2003.01715.x (2003).CAS 
Article 
PubMed 

Google Scholar 
34.Soulsbury, C. D., Iossa, G., Edwards, K. J., Baker, P. J. & Harris, S. Allelic dropout from a high-quality DNA source. Conserv. Genet. 8, 733–738. https://doi.org/10.1007/s10592-006-9194-x (2007).CAS 
Article 

Google Scholar 
35.Soanes, K. et al. Evaluating the success of wildlife crossing structures using genetic approaches and an experimental design: Lessons from a gliding mammal. J. Appl. Ecol. 55, 129–138. https://doi.org/10.1111/1365-2664.12966 (2018).Article 

Google Scholar 
36.Cheng, E., Hodges, K. E., Sollmann, R. & Mills, L. S. Genetic sampling for estimating density of common species. Ecol. Evol. 7, 6210–6219. https://doi.org/10.1002/ece3.3137 (2017).Article 
PubMed 
PubMed Central 

Google Scholar 
37.DeYoung, R. W. et al. Evaluation of DNA microsatellite panel useful for genetic exclusion studies in white-tailed deer. Wildl. Soc. Bull. 31, 220–232 (2003).
Google Scholar 
38.Vedicherla, S. & Buckley, C. T. Rapid chondrocyte isolation for tissue engineering applications: The effect of enzyme concentration and temporal exposure on the matrix forming capacity of nasal derived chondrocytes. Biomed. Res. Int. 2395138, 12 (2017).
Google Scholar 
39.Maličev, E., Kregar-Velikonja, N., Barlič, A., Alibegović, A. & Drobnič, M. Comparison of articular and auricular cartilage as a cell source for the autologous chondrocyte implantation. J. Orthop. Res. 27, 943–948. https://doi.org/10.1002/jor.20833 (2009).CAS 
Article 
PubMed 

Google Scholar 
40.Goel, M., Picciani, R. G., Lee, R. K. & Bhattacharya, S. K. Aqueous humor dynamics: A review. Open Ophthalmol J 4, 52–59. https://doi.org/10.2174/1874364101004010052 (2010).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
41.Park, S. J., Kim, J. Y., Yang, Y. G. & Lee, S. H. Direct STR amplification from whole blood and blood- or saliva-spotted FTA® without DNA purification*. J. Forensic Sci. 53, 335–341. https://doi.org/10.1111/j.1556-4029.2008.00666.x (2008).CAS 
Article 
PubMed 

Google Scholar 
42.Forgacs, D., Wallen, R., Boedeker, A. & Derr, J. Evaluation of fecal samples as a valid source of DNA by comparing paired blood and fecal samples from American bison (Bison bison). BMC Genet https://doi.org/10.1186/s12863-019-0722-3 (2019).Article 
PubMed 
PubMed Central 

Google Scholar 
43.Pfeiler, S. S. et al. Costs and precision of fecal DNA mark–recapture versus traditional mark–resight. Wildl. Soc. Bull. 44, 531–542. https://doi.org/10.1002/wsb.1119 (2020).Article 

Google Scholar 
44.Henry, P., Henry, A. & Russello, M. A noninvasive hair sampling technique to obtain high quality DNA from elusive small mammals. J. Vis. Exp. JoVE. https://doi.org/10.3791/2791 (2011).Article 
PubMed 

Google Scholar 
45.Wirsing, A. J., Quinn, T. P., Adams, J. R. & Waits, L. P. Optimizing selection of brown bear hair for noninvasive genetic analysis. Wildl. Soc. Bull. 44, 94–100. https://doi.org/10.1002/wsb.1057 (2020).Article 

Google Scholar  More

Study sites and examined tree speciesThe study was conducted at three sites in the Andes of southern Ecuador along an elevation gradient at ca. 1000 m (Bombuscaro, Podocarpus NP), ca. 2000 m (San Francisco Reserve) and ca. 3000 m elevation (Cajanuma, Podocarpus NP) in the Provinces of Loja and Zamora-Chinchipe. All sites are located in protected forest areas. At each elevation three permanent 1-ha plots were established in 2018, choosing representative portions of old-growth forest without visible signs of human disturbance (Appendix A1).The forest types at the three sites differ in floristic composition, species richness and structural characteristics49: The premontane rain forest (below 1300 m) at the lowermost site reaches 40 m in height with common tree families being Fabaceae, Moraceae, Myristicaceae, Rubiaceae, and Sapotaceae. It is replaced at 1300–2100 m by smaller-statured lower montane rain forest with Euphorbiaceae, Lauraceae, Melastomataceae, and Rubiaceae as characteristic tree families, and above 2100 m by upper montane rain forest with a canopy height that rarely exceeds 8–10 m. Dominant tree families of the latter forest type are Aquifoliaceae, Clusiaceae, Cunoniaceae, and Melastomataceae. Tree species turnover is complete between premontane and upper montane forest, while a few tree species are shared between lower montane and premontane or upper montane forest types.The climate is tropical humid with a precipitation peak from June to August and a less humid period from September to December. Mean annual temperature decreases with elevation from 20 °C at 1000 m to 9.5 °C at 3000 m, while annual precipitation increases from around 2000 mm at the two lowermost sites to 4500 mm at 3000 m. Typically, there are no arid months with More

In total, 13 culicoid species were found in the present study, with 45.5% of the collected specimens belonging to the Obsoletus Complex while species only occasionally present in previous collections in Germany, accounted for approximately 25% of the sampled individuals. Thus, the species composition is only partly in accordance to earlier studies on the German Culicoides fauna according to which 70 to over 90% of the specimens belonged to the Obsoletus Complex and up to 20% represented members of the Pulicaris Complex, while other culicoid species were present in negligible numbers only12,13. However, previous studies were based on UV-light trap catches12,13,14,15 and targeted active culicoid specimens16. The results obtained in this study are very specific as they represent the species compositions associated with the respective breeding substrates.The gender ratio differed strongly between species, revealing no pattern applicable to all species. The dominance of female Culicoides emerging from breeding sites corresponds to earlier results17,18, even though the sex ratio in the present study showed a much higher proportion of females with 70.7% or a female:male ratio of 2.4:1 than the above studies with 55.6%17 or a female:male ratio of 1.06:118.The evaluation of the diversity of each biotope (excluding the ungrazed meadow where no Culicoides were found) revealed clear differences between the agriculturally used habitats and the more natural biotopes. The Shannon–Weaver index depicted very low diversity for all three studied meadows where biting midges were found. The two meadows (with cattle and sheep) of region 2 reached the lowest possible diversity. This seems plausible as only one species was sampled within each biotope. The meadow with cattle of region 1 revealed at least two species. The Evenness factor of 0.24 depicts the dominance of one of them. The low number of species and unbalanced number of specimens within the biotope result in a low Shannon–Weaver index of 0.24, which describes the poor level of biodiversity.The Simpson index measures the probability that two individuals, randomly selected from a sample, belong to the same species. As only one species was sampled on each meadow from region 2, the probability to choose two specimens which belong to one species is 100% (displayed by the value of D = 1.0). The meadow with cattle of region 2 revealed at least two culicoid species, but the dominance of one species leads to a high Simpson index of 0.92 as well.Opposite to the very low biodiversity of all meadows, the four more natural biotopes of region 3 show an overall high level of biodiversity: according to the Shannon–Weaver index, the level of biodiversity is highest within the AFS (H = 2.96). Compared to the other biotopes of region 3, the AFS revealed by far the highest numbers of culicoid species and specimens. This and the relatively high Evenness factor (E = 0.89) lead to the high H value. The Shannon–Weaver indices for CW and MA are 1.91 and 1.92, respectively. Based on the low numbers of species and specimens in both biotopes, the relatively high H value is mainly caused by its high Evenness values of 0.95 (CW) and 0.96 (MA), respectively. Therefore, the almost equal numbers of all present species leads to the relatively high biodiversity, rather than a high number of species.The Shannon–Weaver index of the DW is the lowest of the four biotopes of region 3 with H = 1.42 and rates this biotope as the one with the lowest diversity of region 3. Though the number of species equal the one of the CW and MA, the higher number of specimens and especially the much lower Evenness factor of 0.71 reduces the H value.Other than the Shannon–Weaver index, the Simpson index rates both, the AFS and the MA, as the two most diverse biotopes. With values of D = 0.13, the probability to randomly select two species of the same species is rather low in both biotopes. As the AFS revealed more than double as many species than the MA, the lower number of caught specimens of the MA must have led to the same biodiversity rate.Study 1—Influence of domestic animals on meadows: up to date, dung-breeding Culicoides have been investigated more thoroughly18,19,20 than most other culicoid species. Most studies have focused on examining selectively either dungheaps or cowpats, rather than conducting a direct comparison between grazed and ungrazed meadows under field conditions. In the present study, we were able to show that the ungrazed meadow seems to be an unsuitable breeding habitat for Culicoides. Therefore, it seems plausible that the suitability of meadows as culicoid breeding sites can be largely, if not completely, attributed to the influence of livestock pasturing.The strong dominance of Obsoletus Complex specimens sampled on grazed meadows is not surprising as this species complex is known to contain typical dung-breeders19,20. The high potential of manure as a breeding substrate has been demonstrated before21,22 and explains the high quantity of Culicoides developing on meadows used by cattle in the present study. While 0.83 midges/sample were found on the meadow with cattle in region 1, only 0.21 midges/sample were collected on the meadow with cattle in region 2. The quantitative differences between these two study sites might be caused by the differing time periods of sampling (April to July for region 1 and August to October for region 2). Previous studies observed population peaks of Obsoletus Complex midges in October, though23, giving reason to expect even higher numbers of midges for region 2 than for region 1, particularly so, as region 2 is an agriculturally dominated area with a higher abundance of potential blood hosts and more suitable breeding habitats than region 1.Compared to the much higher total number of midges emerging from cowpats, sheep dung produced only two specimens. The very low number of midges originating from sheep faeces might be due to the very quick decomposition and desiccation of the rather small droppings, which likely reduces the quality of these remains as culicoid breeding sites. Therefore, it can be assumed that, contrary to pastures with cattle dung, sheep-runs might not play an essential role in promoting the distribution of Culicoides. For modeling approaches, it should be considered, though, that this might only apply to single scattered pieces of faeces as the longer persistence of higher volumes of sheep dung, i.e. on muckheaps, might very likely raise its quality as potential breeding sites as observed by21.All grazed meadows revealed very few culicoid species. Besides members of the Obsoletus Complex, only one individual of C. comosioculatus was found. The present investigation represents a case study though as merely one habitat of each type was sampled. More research to confirm the present results is therefore strongly recommended, even more, as ceratopogonid communities of terrestrial ecosystems have been barely investigated24, with the consequence that breeding sites of Culicoides spp. are still poorly known25.Study 2—Quality of forest-dominated biotopes as culicoid breeding sites: In the present study, the AFS turned out to be very productive as a culicoid breeding site in regards to the number of caught specimens and species diversity. Ten of the 13 collected species were found in the AFS. This is 2.5 times as many species as in the three other biotopes of region 3, which contained four species each in different compositions. Therefore, species-specific requirements for larval development seem to be met for more culicoid species in the AFS than in any of the other study sites.The measured pH values are in accordance to soil analyses conducted in German forests26. As the top layers usually are the most acidic ones, the chosen depth of soil sampling in the present study (upper 0–5 cm) persistently produced low pH values. Additionally, the used solvent (CaCl2) is less sensitive to fast changing weather conditions, but also lowers the measured pH value significantly compared to distilled water26—a solvent often used in earlier studies analyzing the distribution of Ceratopogonidae.The wide variances of the soil factors, especially moisture and organic content, were mainly caused by unequal soil conditions within each biotope rather than changes over time (unpublished data). Nevertheless, the statistical analysis revealed that all four biotopes of region 3 were significantly different from each other regarding the three soil factors. Comparing the means of each soil factor revealed that the AFS contained a higher level of soil moisture, a less acidic pH value and a higher organic content than the other three biotopes of region 3. We could show that significantly more midges (0.4 Culicoides/sample) developed in the AFS compared to the three other biotopes of region 3 with 0.12 (DW), 0.07 (CW) and 0.06 (MA) Culicoides per sample.Previous studies have assumed that the level of moisture be a crucial factor for ceratopogonid development17,20. Also, some studies determined the organic content as pivotal17,27. Our statistical analysis revealed that each soil factor has an impact on the probability of Culicoides to occur. Due to high correlations between the various measured soil factors, it could not be clarified, though, whether they influence the number of specimens, too. But as many culicoid species are known to lay their eggs in batches and previous egg-laying encourages females to oviposit at the same site28, an increase in the probability of biting midge presence should indirectly result in a higher number of specimens, too.The aggregation of larvae in terrestrial habitats29 typically results in a high number of samples completely devoid of midges and an overall low number of specimens sampled by emergence traps30. Thus, the obtained low numbers of collected specimens are not surprising. Nevertheless, emergence traps are still considered to be the best tool for the investigation of breeding site productivity, as it offers a safe assignment of species to their specific developmental sites24,29,31.The Culicoides collected in this study are discussed on species level in regards to existing literature.Culicoides achrayi was found in the AFS. A swamp as a breeding site32 and soil located in stagnant water22 have previously been described for this species. We confirm June as the time of emergence32 and add that C. achrayi co-exists with C. pulicaris.Culicoides albicans was collected in the AFS and DW. Specimens hatched from late April to mid-June, representing one generation per year. We confirm co-habitation with C. pictipennis and C. kibunensis11,33 and the preference for very humid substrates which has been described for the wettest parts of boglands5,34 and for artificially waterlogged soil11. Our results show, that C. albicans larvae can tolerate medium moisture levels, too. The mean organic content of their developmental sites reached from moderate to high, and the pH values lay between strong and ultra-acidic.Culicoides comosioculatus was found on the meadow with cattle dung in mid-June. As only one individual (a gravid female with the presumed intention to oviposit) was collected and no literature regarding breeding sites of this species could be found, our finding only indicates that this species might possibly develop in animal dung although in extremely low numbers.Culicoides grisescens was found within the AFS, the CW and the DW from late May until mid-July. Kremer35 listed soils of swamps and boggy grasslands as developmental sites. We collected C. grisescens in three different biotopes with wide variances of the mean moisture level, mean organic content and mean pH value, which reveals the wide tolerance range of this species towards these three soil factors.Culicoides impunctatus was collected in the AFS and the CW from late May to mid-July, representing one generation per year. This finding differs from earlier observations of two generations per year in Scotland36. Previous studies described breeding sites as acidic, oligotrophic grasslands, swamps, boglands or marshes, often of a peaty consistence5,10,33,34,37 and with soil pH values of 5.0–6.5 (dissolved in distilled water)37. This matches the pH values of the AFS in the present study (lower, but dissolved in CaCl2), but excludes the much lower pH values of the CW. The range considered suitable for C. impunctatus larvae should therefore be extended downwards to as low as pH 2.9–3.9 (CaCl2). We found C. impunctatus in two biotopes comprising a wide variance regarding soil moisture and organic content, which illustrates the wide tolerance range of this species. Individuals of C. impunctatus co-exist with Obsoletus Complex specimens as both were collected within the same sample in the AFS.Culicoides kibunensis was collected in the AFS and MA, which matches earlier observations depicting swamps of eutrophic fresh water bodies17,34, soil of stagnant water bodies22 and acidic grasslands in considerable distances to swamps33 as breeding sites. The AFS and MA revealed pH values between 3.4 and 5.4. Soil moisture and organic content displayed wide variances. All specimens hatched from late May to mid-June. Culicoides kibunensis was found to co-exist with C. albicans as observed by Kettle33. Earlier observations of co-habitations with C. obsoletus s.s. and C. pallidicornis5,34 could not be confirmed.Obsoletus Complex members were present in all study sites except for the ungrazed meadow. In the grazed meadows, Obsoletus Complex midges emerged almost throughout the entire sampling period except for the month of September. Two peaks were observed, one in June/July and a smaller one in October. As in the grazed meadows, the biotopes of region 3 also revealed two generations, but emerging at a slightly earlier time of the year with one peak in May/June and the other one in September/October.Members of the Obsoletus Complex are known to be generalists regarding their choice of breeding sites. Only the identified member species, C. chiopterus and C. obsoletus s.s., are considered here.Culicoides chiopterus was exclusively found on meadows grazed by cattle, which is in accordance to several earlier studies as this species is described as a dung-breeding species developing in cowpats and horse droppings5,34,35,38.Culicoides obsoletus s.s. was mostly sampled in the AFS. Only one individual was collected on a meadow grazed by cattle. Previous descriptions of breeding sites differed widely. Acidic grasslands in considerable distance to bogs/swamps33 and leaf litter compost5,35 could not be confirmed in the present study, although the MA and AFS were of a comparable character. While Uslu and Dik17 could not find any C. obsoletus s.s. in wet organic matter-rich soil, we collected most specimens of this species in the AFS and can therefore confirm previous findings11,29,32,39. The time of C. obsoletus s.s. activity in Germany (April–October) as described by Havelka32 agrees with our observations.Culicoides pallidicornis was found in the MA in late June. This species revealed the smallest variances of all sampled biting midge species regarding the three soil factors, using soil with pH values of 3.6–5.0 (CaCl2) and a relatively low level of moisture. This contradicts earlier observations where C. pallidicornis developed in the mud of eutrophic fresh-water swamps5. While C. pallidicornis larvae are known to co-exist with C. kibunensis5, we can add C. subfagineus to share the same developmental site.Culicoides pictipennis was collected in the DW and, to a minor part, in the AFS. The preferred physicochemical breeding conditions were ultra to extremely acidic with a medium moisture level and a moderate to slightly increased organic content. This differs from previous studies, which have found this species to develop only at the margin of stillwater bodies like pools and ponds, and the littoral of lakes or in artificially waterlogged soil11,32,34. Havelka32 observed C. pictipennis between May and June, while in our investigation the first specimen emerged as early as mid-April. We can confirm the co-existence of C. pictipennis and C. albicans as previously observed by Harrup11.Culicoides pulicaris was sampled in the AFS from late June until September, which agrees with observations denoting May to September as the activity time of this species32. Culicoides pulicaris seems to prefer breeding substrates with a high moisture level and a high organic content, as previously described17,32,34. We can add that C. pulicaris breeds in soil showing pH values at least between 4.0 and 5.4. We collected C. pulicaris together with C. achrayi and found it to simultaneously emerge from one biotope with C. obsoletus s.s. Additionally, we can confirm the co-existence of C. pulicaris with C. punctatus5,40, since both species have similar breeding habitat preferences11.Culicoides punctatus was sampled in the AFS and, to a minor part, in the CW. Time of emergence was from mid-June to late September, which is in accordance with earlier observations listing April-August and October as times of activity32. In the present study, a strong preference for swampy conditions with soil of high moisture, high organic content and a strong to very strong acidity was found. This is in agreement to previous findings11,32,41. The co-existence of C. punctatus with C. pulicaris is well known5,40 and can be confirmed once more. Additionally, we found C. punctatus to co-occur with C. subfasciipennis.Culicoides subfagineus was caught in the MA in late June. The soil was oligotrophic and contained a relatively low moisture level with pH values between 3.6 and 5.0. The first record of this species in Germany was in 2014, when C. subfagineus was observed to attack cattle42.Culicoides subfasciipennis was sampled in mid-June in the AFS. The time and choice of breeding site are in accordance to previous findings17,32. Breeding conditions for the only individual collected revealed a medium soil moisture factor, a pH value of 5.2 and a medium organic content. The species was found to co-develop with C. punctatus. More

Mosquito lines and maintenanceTo introduce Wolbachia into Brazilian Ae. aegypti, an Australian line infected with the wMel strain21 was backcrossed for 8 generations to a natural mosquito population of Rio de Janeiro, Brazil24. Following the genetic background introgression, additional crosses and knockdown resistance (kdr) screening were undertaken to replicate natural insecticide resistance profiling and generate the line wMelRio. To assure a minimal variation in this profiling overtime, and sustain a homogeneous genetic background, wMelRio colony was refreshed with 10% wild males once in every five generations40.To maintain wMelRio, immatures (i.e. larval stages L1 to L4) were reared in dechlorinated water, at 28 °C, and fed Tetramin Flakes (Tetra GmbH, Herrenteich, Germany) until pupal formation. Following adult emersion, groups of 1000 females and 800 males were sorted and kept in BugDorm cages (MegaView Science Co Ltd., Taiwan) at 25 °C, with 10% sucrose solution ad libitum. Every three days, females were fed human blood (from blood donation centers; see details under ‘ethical considerations’), through Hemotek artificial feeders (Hemotek Ltd, UK). Note that, to avoid arboviral contamination of our colony, all blood samples were formerly tested negative for DENV, ZIKV, CHIKV, MAYV and YFV by multiplex qPCR assays36,68. Egg-laying was induced by placing dampened strips of filter paper (i.e. partially immersed in water-containing cups) inside the cages for 2–3 days, after which they were gradually dried at room temperature. Strips loaded with eggs (i.e. ovistrips) were kept at room temperature until further use, either for colony maintenance or field release. Eggs older than 40 days were discarded due to a decay in overall quality60.Egg releasesMass-reared wMel-infected Brazilian Ae. aegypti, wMelRio, were released as eggs in Jurujuba (22°56′ 00″ S, 43°07′ 00″ W), a lower-middle-class community in the city of Niterói (state of Rio de Janeiro, Brazil). Located by the shores of Guanabara bay, this community has grown from a typical fisherman settlement, with informal occupancy, to a total population of 2797 residents in 890 houses. Jurujuba encompasses a total area of 2.53 km2, divided into seven smaller sectors (i.e. sub-areas or localities within the neighborhood): Ponto Final, Várzea, Brasília, Cascarejo, Praia de Adão e Eva, Peixe-Galo and Salinas.wMelRio eggs were released in the field through special deployment devices, referred to as mosquito release containers (MRCs), which consisted of small white plastic buckets (19 cm height × 18 cm top diameter × 15.5 cm base diameter) with four small holes on the side wall, only a few centimeters away from the top lid. Each MRC was loaded with 1 L of water, 0.45 g of Tetramin Tropical Tablets (i.e. one and a half tablet) (Tetra GmbH, Herrenteich, Germany) and an ovistrip containing approximately 150–300 eggs. After six to seven days, about 80% of the immatures were pupae, and after 11 to 12 days, most of the adults had already emerged and left the device from the wall holes. Every 15 days, MRCs were checked and reloaded so that another rearing and release cycle could take place. Release sites were spatially distributed as evenly as possible (Supplementary Fig. S1), so as to maximize the spread of Wolbachia-harboring individuals and promote mating with their wild peers. The release strategy was optimized by splitting the sites into two groups, A and B, with alternate MRC loading schedules. Thus, while MRCs from group A were releasing adults, those from group B were being loaded with new ovistrips, water and food. In the following week, an opposite situation occurred, with MRCs from group B releasing adults. The release schedules, as well as the number of allocated MRCs, varied according to each Jurujuba’s sector (Supplementary Table S1).Ethical considerationsAll methods were carried out in accordance with relevant guidelines and regulations. Study protocol for Wolbachia field release was approved by the National Research Ethics Committee (CONEP, CAAE 02524513.0.1001.0008) and three government agencies: IBAMA (Ministry of Environment), Anvisa (Ministry of Health) and MAPA (Ministry of Agriculture, Livestock and Supply) to obtain the RET (Special Temporary Registry, 25351.392108/2013-96). Prior to mosquito releases, an informed consent was obtained from 70% of Jujuruba households. Also, a written informed consent was obtained from households that hosted BG-sentinel mosquito traps.For the maintenance and mass-rearing of Wolbachia-infected Ae. aegypti, adult females were fed human blood from a donation center (Hospital Antonio Pedro, Rio de Janeiro State University), with supporting regulatory approval (CONEP, CAAE 59175616.2.0000.0008) We only used blood bags which would have been discarded by the donation center, mainly due to insufficient volume to meet their quality assurance policy. Samples had no information on donor’s identity, sex, age and any clinical condition, but were tested negative for several diseases, including Hepatitis B, Hepatitis C, Chagas disease, syphilis, HIV and HTLV, as part of the Brazilian Government routine screening.For vector competence assays, human blood was obtained from Fundação Hemominas as part of a research agreement with Instituto René Rachou (Fiocruz Minas) (OF.GPO/CCO-Nr224/16).
Wolbachia field monitoring and density level assessmentAe. aegypti field population was monitored with BG-Sentinel traps (Biogents AG, Regensburg, Germany), spread across Jurujuba in a semi-homogeneous fashion (Supplementary Fig. S2, Supplementary Table S2, Supplementary Datasheet S1). These monitoring sites were chosen among suitable households who formally agreed with hosting of a trap, and had to be reallocated according to necessity (i.e. household quits hosting the trap). Working traps were checked weekly by removing the catch bags (e.g. small meshed envelopes placed inside the BG-Sentinels to collect trapped insects) and bringing them to the laboratory for species identification and Wolbachia screening. Catch bags were barcoded according to the trap ID and site, so as to create a pipeline for field samples.Screening for Wolbachia in Ae. aegypti samples was undertaken by qPCR. Briefly, individual DNA was extracted by homogenizing head/thorax body parts in Squash Buffer (10 mM Tris–Cl, 1 mM EDTA, 25 mM NaCl, pH 8.2) supplemented with Proteinase K (200 ug/ml) and incubating at 56 °C for 5 min. Extraction ended by enzyme inactivation at 98 °C for 15 min. DNA amplifications were carried out with FastStart Essential DNA Probes Master (Roche), using specific primers and probes to Wolbachia pipientis WD0513 and Ae. aegypti rps17 markers (Supplementary Table S3). Thermocycling conditions were set on a LightCycler 96 Instrument (Roche), as follows: 95 °C for 10 min (initial denaturation), and 40 cycles of 95 °C for 15 s and 60 °C for 30 s. Samples were analyzed using absolute quantification, by comparison to serial dilutions of either gene product, cloned and amplified in the pGEMT-Easy plasmid (Promega). Negative control samples were normalized between plates, and were used as reference to determine a minimum threshold for positive samples.DENV and ZIKV isolation and replication in mosquito cellsZIKV was kindly provided by Instituto Aggeu Magalhães (IAM, Fiocruz) through viral isolation of a symptomatic patient sample from Recife (PE, Brazil) in 2015 (ZikV/H.sapiens/Brazil/BRPE243/2015). DENV was sourced following a viral isolate from a patient sample diagnosed with Dengue type 1 in Contagem (MG, Brazil), also in 2015 (Den1/H.sapiens/Brazil/BRMV09/2015). Both ZIKV and DENV samples were accompanied by patients’ written consent (CONEP, reference number 862.912), being further catalogued into the national database of genetic patrimony and associated knowledge (SISGEN, access number AA1D462).In vitro culture of viral particles were done as previously described36. Briefly, ZIKV and DENV were replicated in Aedes albopictus C6/36 cells, grown at 28 °C in Leibovitz L-15 medium (ThermoFisher) supplemented with 10% fetal bovine serum (FBS) (ThermoFisher). After seven days, supernatants were harvested and virus titers were assessed, first by Reverse Transcription (RT)-qPCR, and later by plaque assay with VERO cells grown under 37 °C, 5% carbon dioxide, in Dulbecco’s Modified Eagle Medium (DMEM) (ThermoFisher) supplemented with 3% Carboxymethylcellulose (Synth) and 2% FBS.Vector competence assaysTo perform vector competence assays with field samples of Ae. aegypti, ovitraps were mounted in both Ponto Final (Jurujuba) and Urca, a Wolbachia-free area in Rio de Janeiro. Ovitraps were collected from the field over 13 weeks, from April to June 2017, which corresponds to the time-frame between 14 and 16 months along the post-release phase in Ponto Final. Once in the insectary, eggs samples were reared to the adult stage in a controlled insectary environment (refer to ‘mosquito lines and maintenance’ for details).For virus challenging assays through oral-feeding, young females (4–6 days old) were starved for 20 to 24 h, and subsequently offered culture supernatant containing ZIKV or DENV mixed with human red blood cells (2:1 ratio), using an artificial membrane feeding system36. It is important to mention that, as for the colony maintenance protocol, blood samples used here were also submitted to quality control prior to its use in the assays, mainly due to putative arbovirus contaminations which could affect the experimental outputs. Likewise, all samples were tested negative for DENV, ZIKV, CHIKV, MAYV and YFV by multiplex qPCR assays36,68. Oral-infections were performed twice for each virus. ZIKV was offered first from fresh (initial virus titer of 4.8 × 108 PFU/mL) and second from frozen culture supernatant (initial virus titer of 7.6 × 106 PFU/mL). In contrast, DENV was offered from fresh supernatants only (virus titers of 2 × 106 PFU/mL and 6.5 × 107 PFU/mL), since frozen versions failed to infect. Specimens were allowed to feed for one hour, after which engorged females were selected and maintained with 10% sucrose solution ad libitum, during the whole extrinsic incubation period. At 14 days post-infection (dpi), viral loads were assessed in heads/thorax extracts by RT-qPCR (refer to ‘Viral diagnosis’ for more details).For saliva-mediated virus challenging assays, ZIKV and DENV pre-exposed females (14 dpi) from Jurujuba (Wolbachia +) and Urca (Wolbachia −) were starved for about 16 h (overnight) before being knocked down and kept at 4 °C for wings and legs removal. Salivation was induced by introducing a 10 µL sterile filter tip, pre-loaded with 5 µl of a solution [30% sucrose (w/v) diluted in 50% fetal bovine serum (FBS) and 50% DMEM medium], into the mosquito proboscis for 30 min. Saliva samples were individually collected, and 276 nL was intrathoracically injected into young naive females (Urca) using a Nanoject II hand held injector (Drummond), as previously described36,68. Each saliva sample was used to inoculate 8–14 naïve Wolbachia-free Ae. aegypti specimens, of which 8 were screened for infective particles. ZIKV and DENV were quantified by RT-qPCR at 5 dpi and 7 dpi, respectively (refer to ‘Viral diagnosis’ for more details). Overall Intrathoracic Saliva Infection index (OISI) was obtained by averaging the percentages (± SD) of infected individuals in each group.Viral diagnosisTo identify ZIKV and DENV particles in individual samples, whole specimens were processed into head/thorax homogenates for RNA/DNA extraction with the High Pure Viral Nucleic Acid Kit (Roche), according to manufacturer’s instructions30. Extracted samples were diluted in nuclease-free water to a concentration of 50 ng/μL. ZIKV, DENV and Wolbachia levels, in vector competence assays, were quantified by RT-qPCR using TaqMan Fast Virus 1-Step Master Mix (ThermoFisher) and specific primers and probes (Supplementary Table S3). Reactions were run on a LightCycler 96 Instrument (Roche), using the following thermocycling conditions: 50 °C for 5 min (initial RT step), 95 °C for 20 s (RT inactivation/DNA initial denaturation), and then 40 cycles of 95 °C for 3 s and 60 °C for 30 s. Each RNA/DNA sample was used in two reactions, one with ZIKV, DENV or Wolbachia primers, and another with Ae. aegypti rps17 endogenous control30. Absolute quantification was achieved by comparing amplification profiles with standard curves generated by serial dilutions of their respective gene products, amplified from a cloned sequence in pGEM-T Easy vector (Promega). Negative control samples (no virus RNA) served as reference to fix a minimum threshold for positive ones. ZIKV and DENV loads were defined as their copy number per sample (head/thorax or saliva), while Wolbachia loads were relative quantifications to the rps17 reference gene. Here, it is worth noting that, while Wolbachia titer is naturally variable and dependent on its whole-body density, the overall expression of rps17 is stable and particularly suitable for internal controls in assays with adult females69, as demonstrated previously by us and others30,62,68.Map creation and source codesThe satellite image map of Jurujuba was created with ArcGIS Desktop 10.7 (Esri Inc., https://www.esri.com/en-us/arcgis/products/arcgis-desktop/overview) using Google Earth (Google LLC) source code, under the license and in accordance with the fair use described in ‘https://about.google/brand-resource-center/products-and-services/geo-guidelines/’. Maps with geotagged MRCs and BG-Sentinel traps were created with ArcGIS Desktop 10.7 and OpenStreetMap source code (OpenStreetMap contributors), under the license CC-BY-SA 2.0.Statistical analysesGraphs and statistical analyzes were performed in GraphPad Prism 8 (GraphPad Software Inc., https://www.graphpad.com). Kruskal–Wallis test followed by Dunn’s post-hoc multiple comparisons were used to analyze Wolbachia density data from field-collected and colony samples. ZIKV and DENV loads in head/thorax extracts, from both oral and saliva-challenging samples, were compared using the Mann–Whitney U test. For all statistical inferences, ⍺ was set to 0.05. More

Ethics statementAll procedures were conducted in accordance with the current laws in Thailand on experimental animals and were approved by the safety management committee for experiments of the Laboratory Animal Center, Chiang Mai University (Project Number 2561/FA-0001). The study also followed the recommendations in the ARRIVE guidelines.Species-specific primer designAll the DNA tissue analysed originated from the mucus of the individual giant snakehead. Total DNA was extracted from the mucus sample using the Qiagen DNeasy Blood and Tissue Kit (Qiagen, Valencia, CA). Extracted DNA was used as a template for qPCR assay together with synthetic fragments. DNA samples were quantified using a Qubit fluorometer (Life Technologies) calibrated with the Quant-iT dsDNA HS Assay following the manufacturer’s instructions. For each replicate, 3 µL volumes were measured.Species-specific primers and a minor-groove binding (MGB) probe incorporating a 5′ FAM reporter dye and a 3′ non-fluorescent quencher were designed to amplify an 127 bp targeting within the 16S region for the giant snakehead (C. micropeltes), using Primer Express (V3.0, Life Technologies; Table 3). Probe and primer sequences were matched against the National Centre for Biotechnology Information (NCBI, http://www.ncbi.nlm.nih.gov/) nucleotide database with BLASTn (Basic Local Alignment Search Tool) to confirm the species’ specificity for the giant snakehead in silico assays.Table 3 Details of species-specific primers and the probe designed to amplify a 127 bp fragment of the 16S region of Channa micropeltes (Cuvier, 1831).Full size tableTo ensure that the assay only amplified the giant snakehead, it was deployed on a closely related species commonly found in Thai freshwater environments using conventional PCR amplification and visualization on a 1.5% agarose gel stained with SYBR Safe DNA Gel Stain (Life Technologies).qPCR assayThe qPCR assay was deployed using Environmental Master Mix (Applied Biosystems) on mucus samples from the giant snakehead and related species to ensure the species specificity to the qPCR assay. In addition, eDNA qPCR assay for the giant snakehead, a water sample collected from tank at Phayao Freshwater Aquarium (Phayao Inland Fisheries Research and Development Center) was known to have only the giant snakehead was included as a positive control for the presence of amplifiable eDNA in water samples. The tank contains around 4.5 m3 of water with one individual of giant snakehead resides in the tank (the fish is about 60–70 cm in length).All eDNA qPCR amplifications were performed in three replicates in a final volume of 20 µL, using 10.0 µL of 2 × TaqMan Environmental Master Mix 2.0 (Thermo Fisher Scientific), 2.0 µL of DNA template, 900 nM each of the F/R primers, and 125 nM of the probe. Samples were run under the following conditions: an initial 10 min incubation at 95 °C followed by 50 cycles of denaturation at 95 °C for 15 s and annealing/extension at 60 °C for 1 min. Negative controls with all PCR reagents but no template (three replicates) were run in parallel to assess potential contamination. The quantification cycle (Cq) was converted to quantities per unit volume using the linear regression obtained from the synthesized target gene standard curve (Integrated DNA Technologies Pte. Ltd., Singapore). The giant snakehead eDNA concentrations were then reported as copies/mL. The limit of detection (LOD) and the limit of quantification (LOQ) were also measured using the standard dilution series of synthesized target gene fragment with known copy numbers. A dilution series containing 1.5 × 101 to 1.5 × 104 copies per PCR tube were prepared and used as quantification standards. The calculation of LOD and LOQ was done using published R script by Klymus et al.26.DNA extraction from the filtersDNA trapped on the filters obtained from the aquarium experiments and field collections were extracted using Qiagen DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany) using a protocol modified from the manufacturer’s protocol with the following changes: the DNA from all samples were eluted twice with 25 µL AE buffer, in a total volume of 50 µL to obtain a more concentrated eDNA solution. The volume of ATL buffer (360 µL), Proteinase K (40 µL), AL buffer (400 µL) and Ethanol (400 µL) were doubled.Aquarium experimentAn aquarium experiment was used to test the extent to which qPCR of water samples can detect eDNA of giant snakehead at low simulated densities. The juvenile giant snakehead was obtained from the fish store and transported to a laboratory at Chiang Mai University. The giant snakeheads were then held in separate 120 L plastic holding containers in which the water was continuously filtered. The fish were fed frozen shrimp/commercially available flake fish food three times a week, and were held at 23 ± 1 °C.The sensitivity of eDNA detection in the aquaria was evaluated by conducting three aquarium experiments using plastic tanks (30 × 45 × 25 cm) filled with 120 L of aged-tap water. The water in the tanks was continuously aerated through a filter. In each experiment, the giant snakeheads were randomly assigned to the tanks (10 individuals per tank). The average size of the snakeheads was 9.7 cm (body length ranging from 9.1 to 10.6 cm). The average weight was 8.15 g (ranging from 6.7 to 10.6 g). The water in the tanks was maintained at 23 ± 1 °C. A 300 mL water sample from each tank was collected at each time point (0, 3, 6, 12, 24, 48, 72, 96, 120, 144, and 168 after removal of the fishes from the tanks) in triplicate. Collected water was filtered on a GF/F filter (0.7 μm Whatman International Ltd., Maidstone, UK). The eDNA from each sample solution was extracted using a Qiagen DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany) in a final volume of 50 µL, detailed in DNA extraction from the filters. To confirm the absence of the giant snakehead eDNA in the water prior to the experiments, three tanks without giant snakehead were prepared and water sample was collected and treated as described above.Real-time PCR was performed with the species-specific primers and probe set using a Rotor-Gene Q system (Qiagen, Hilden, Germany). The reaction conditions were the same as described in qPCR assay. Three replicates were conducted for each sample including the negative PCR control and positive control.eDNA field collectionWater samples were collected at 6 points within Kwan Payao according to the survey locations of the Inland Fisheries Research and Development Center. Additional water samples were collected from 11 and 6 locations in Ing River where water flowed into and out of Kwan Payao, respectively (Fig. 1). To avoid contamination, all field equipment was sterilized using 10% bleach, UV-Crosslinker or autoclaved and sealed prior to transport to the study site, and a separate pair of nitrile disposable gloves were used for each sample. At each site, water samples were collected three replicate in bucket that had been previously decontaminated with a 10% bleach rinse followed by two distilled water rinses.In total, water samples were collected from 6 sites (in Kwan Phayao) and from 17 sites (in the Ing River) from 15th February to 5th March 2019, the middle of the dry season. Each site was sampled in triplicate, 300 mL samples of water were collected and filtered on GF/F filter (0.7 μm Whatman International Ltd., Maidstone, UK). In total, 306 water samples were collected from the surface water of lakes and rivers. For every sampling day, deionised water (300 mL) was filtrated as a negative control. The water samples and real-time PCR was processed as described above in qPCR assay. More