Philippa Kaur

Women in northern Mumbai, India, have planted mangrove saplings to protect the area against rising sea levels.Credit: Mahendra Parikh/Hindustan Times via Getty

Projects that manage, protect and restore ecosystems are widely viewed as win–win strategies for addressing two of this century’s biggest global challenges: climate change and biodiversity loss. Yet the potential contribution of such nature-based solutions to mitigating climate change remains controversial.Decision-makers urgently need to know: what role do nature-based solutions have in the race to net-zero emissions and stop further global temperature increases?Analyses of nature-based solutions often focus on how much carbon they can remove from the atmosphere. Here, we provide a new perspective by modelling how these solutions will affect global temperatures — a crucial metric as humanity attempts to limit global warming.Our analysis shows that nature-based solutions can have a powerful role in reducing temperatures in the long term. Land-use changes will continue to act long past the point at which net-zero emissions are achieved and global temperatures peak (known as peak warming), and will have an important role in planetary cooling in the second half of this century. Before then, nature-based solutions can provide real but limited mitigation benefits. Crucially, the more ambitious the climate target, the shorter the time frame for such solutions to have an effect on peak warming.In other words, nature-based solutions must be designed for longevity. This means paying closer attention to their long-term carbon-sink potential, as well as their impacts on biodiversity, equity and sustainable development goals. It also means continuing to limit global warming through other methods, from decarbonization to geological storage of carbon dioxide.Our model reinforces the conclusion that an ambitious scaling-up of nature-based solutions needs to be implemented fast and thoughtfully — and not at the expense of other measures.Win–winsThe world is currently likely to hit 3 °C of warming above pre-industrial levels by 2100 (although recent policy announcements from the United States and China could reduce this). The 2015 Paris climate agreement aims to limit the global temperature rise this century to well below 2 °C, and, ideally, to 1.5 °C. There is no date for either goal, beyond the “end of this century”. The metric that matters most is the peak temperature, with more-aggressive efforts required to stay below 1.5 °C of warming than for the 2 °C target.
Emissions: world has four times the work or one-third of the time
It is impossible to achieve the needed reduction in peak warming solely through cuts to greenhouse gases, because emissions from certain sectors, such as agriculture and some heavy industry, cannot be driven to zero any time soon. For this reason, we also need to remove greenhouse gases from the atmosphere on an unprecedented scale1.There are various options for doing this. For example, when biomass vegetation is burnt for energy, the emitted CO2 can be retained and stored underground. This process, known as bioenergy with carbon capture and storage (BECCS), requires vast areas of land — compromising food security and biodiversity — as well as time to develop on a large scale. Other options involve industrial machines that capture CO2 from the air; these are currently nascent, expensive technologies.A subset of nature-based solutions can be used specifically to limit warming. These ‘natural climate solutions’ aim to reduce atmospheric greenhouse-gas concentrations in three ways. One is to avoid emissions by protecting ecosystems and thus reducing carbon release; this includes efforts to limit deforestation. Another is to restore ecosystems, such as wetlands, so that they sequester carbon. The third is to improve land management — for timber, crops and grazing — to reduce emissions of carbon, methane and nitrous oxide, as well as to sequester carbon (see ‘Three steps to natural cooling’).

Source: S. Jenkins et al. Geophys. Res. Lett. 45, 2795–2804 (2018).

Decades of work provide strong evidence that nature-based solutions can deliver many local ecological and socio-economic benefits2. Restoring a forest next to a stream, for example, might reduce flooding, improve carbon storage and support fisheries. Growing recognition of such benefits means that interest in nature-based solutions is soaring: they can help people adapt to climate change, achieve sustainable development goals, protect biodiversity and mitigate climate change3.Quantifying nature’s roleThere is still debate around how much nature-based solutions can contribute to achieving net-zero targets by mid-century. This is because results have been estimated across a range of objectives, time frames and model assumptions4,5 (see Supplementary information; SI). Some researchers say that tree restoration is the most effective climate-change solution we have available6 (this in itself has been robustly contested); others argue that nature-based solutions won’t be nearly as fast or as effective as is often stated7.Part of the reason for the impasse is this: many well-known papers discuss the annual carbon uptake possibilities of nature-based solutions; they do not discuss their cooling impact year on year. Because the Paris agreement is framed in terms of temperature, we argue that this gap is critical: researchers need to know how nature-based solutions will affect global temperature.To model this, we consider an ambitious but realistic scenario — an update to previous estimates by one of our co-authors (B.W.G)4,8,9. This scenario considers only those projects for nature-based solutions that are constrained by many factors: they are cost-effective (costing less than US$100 per tonne of CO2 equivalent); ensure adequate global production of food and wood-based products; and involve sufficient biodiversity conservation. They also respect land tenure rights and don’t change the amount of sunlight reflected from Earth, or albedo (see SI). In our scenario, nature-based solutions that avoid emissions ramp up quickly — by 2025 — and absorb carbon while avoiding emissions at a rate of 10 gigatonnes of CO2 per year (Gt CO2 yr−1). This rises to 20 Gt CO2 yr−1 in the most ambitious scenario (peak warming of 1.5 °C by 2055), in which we assume a higher price of carbon. The 10-Gt value is cost-contained. But we also account for 30 years of higher-priced nature-based solutions in the 1.5 °C scenario (up to $200 per tonne of CO2 equivalent; see SI). For comparison, 10 Gt CO2 yr−1 is more than the emissions from the entire global transportation sector.

Instituto Terra, an initiative in Aimorés, Brazil, is restoring a devastated ecosystem.Credit: Christian Ender/Getty

Achieving 10 Gt CO2 yr−1 of mitigation in this way would involve stopping the destruction of ecosystems worldwide (including 270 million hectares of deforestation); restoring 678 million hectares of ecosystems (more than twice the size of India); and improving the management of around 2.5 billion hectares of land by mid-century4. This is ambitious, but it is important to note that the bulk of land required (85%) comes from improving management of existing lands for agriculture, grazing and production forest without displacing yields of food, wood-based products or fuel (see ‘Three steps to natural cooling’).These estimates come with caveats (see SI). The role of nature-based solutions could be larger if one considers, for example, their impacts on other greenhouse gases besides CO2. This could represent an additional amount of roughly 1–3 Gt CO2 equivalent yr−1 of climate mitigation. Alternatively, the contribution of such solutions might be smaller in the long term, if the carbon drawdown from land-based interventions decreased over time. This could happen if these natural sinks became saturated or were affected by climate impacts such as forest fires. These caveats are not included in our estimates.We then modelled how this level of nature-based solutions would affect global temperature up to 2100 (see ‘The long game’ and SI). We looked at illustrative pathways from the Intergovernmental Panel on Climate Change, in which peak warming is constrained to 1.5 °C or 2 °C, and ran these scenarios with the added contribution of nature-based solutions as described. These pathways include BECCS, but no nature-based solutions beyond some avoided deforestation.Taking the temperatureOur analysis shows that implementing this level of nature-based solutions could reduce the peak warming by an additional 0.1 °C under a scenario consistent with a 1.5 °C rise by 2055; 0.3 °C under a scenario consistent with a 2 °C rise by 2085; and 0.3 °C under a 3 °C-by-2100 scenario (see ‘The long game’).

Adapted from Fig. SPM.1 of Ref. 1

The most significant contribution nature-based solutions can make to mitigating the peak temperature is in the 2 °C scenario. In a more ambitious 1.5 °C scenario, there isn’t enough time for nature-based solutions to have as great an impact on peak warming. In the 3 °C scenario, several issues constrain the impact of nature-based solutions, including the limited ability of ecosystems to absorb carbon in a warmer world.Overall, the mitigation potential of nature-based solutions remains small compared to what can be achieved by decarbonizing the economy. Yet, assuming that decarbonization takes place, nature-based solutions can still suppress a chunk of the warming (see SI).Crucially, nature-based solutions cool the planet long after the peak temperature is reached. In the 1.5 °C scenario, they take a total of 0.4 °C off warming by 2100 — four times their suppression to the 2055 peak temperature (see SI, Table S2).
Restoring natural forests is the best way to remove atmospheric carbon
Achieving these significant long-term benefits requires several things. Nature-based solutions of good quality must be scaled up rapidly — and not at the expense of other robust strategies. Long-term geological storage of CO2, for example, will need to be ramped up significantly in the next decade as technologies mature and prices fall. The long-term benefits of nature-based solutions also depend on warming being held in check. The increased frequency and intensity of impacts such as wildfires can undermine ecosystems and their capacity to store carbon or provide other benefits to society.Ecosystems that are protected and carefully managed — such as intact peatlands and old-growth tropical rainforests — are very likely to continue to store carbon for thousands of years. These are also more resilient to climate extremes and pathogens.The right metricsRestoration of forest cover is widely considered the most viable near-term opportunity for carbon removal. Unfortunately, some of this enthusiasm has been used to promote plantation forestry — growing trees of a limited variety of ages and species (for example, in monoculture plantations) does not have the same carbon benefits as maintaining an intact forest ecosystem10.One serious problem is that some nature-based solutions, as currently implemented, can have unintended and unwanted consequences. For example, an area of 34,007 hectares of intact forest ecosystem in Cambodia became a logging concession, with much of it replaced with an acacia monoculture. This was the first large-scale reforestation project to be funded in Cambodia in the context of climate-change mitigation. The project resulted in unethical ecological devastation, affecting 1,900 families in the area11.Similarly, Chilean government subsidies for new plantations of pine and eucalyptus have resulted in plantations expanding by 1.3 million hectares since 1986, with an associated sequestration of about 5.6 million tonnes of carbon. However, regulations stating that expansion cannot happen at the expense of native biodiverse forests were not enforced, resulting in large-scale reductions in native forest cover. Clearing of the original forest has resulted in a net decrease of approximately 0.05 million tonnes of stored carbon since 198612.These examples show how a singular focus on rapid carbon sequestration as the metric of success for land-based climate mitigation can result in perverse outcomes. Activities should be evaluated and monitored with the right metrics, to account for the multitude of benefits they provide in the long term.
Adopt a carbon tax to protect tropical forests
To ensure long-term resilience, projects involving nature-based solutions should adhere to four high-level principles (see nbsguidelines.info). First, nature-based solutions are not an alternative to decarbonization; second, they need to involve a wide range of ecosystems; third, they should be designed in partnership with local communities while respecting Indigenous and other rights; and, finally, they must support biodiversity, from the level of the gene to the ecosystem. In addition, the Oxford principles13 for high-quality offsets call for safe and durable CO2 removal and storage for every tonne of CO2 emitted. Metrics of success should include those for carbon dynamics, biodiversity across multiple trophic levels, and socio-economic factors such as women’s empowerment and youth employment.There are many examples of good-practice projects (see also case studies by the University of Oxford’s Nature-based Solutions Initiative, where N.S. and C.A.J.G. work). For example, mangrove forests in eastern India that have been protected from deforestation since 1985 have been shown to protect coastal regions from the negative impacts of cyclones much better than artificial defences do, while also soaking up carbon14. In Sierra Leone’s tropical rainforest, cocoa agroforestry — where cocoa is planted with trees for shade, alongside pineapples, chillies and maize (corn) as an additional source of food and income — has been shown to produce cocoa sustainably while diminishing forest clearance. One agroforestry project in the Gola Rainforest National Park, initiated 30 years ago, has increased biodiversity and the profitability of crops while saving an estimated 500,000 tonnes of carbon each year through sequestration and avoiding deforestation.Invest wiselyThis much is clear: we urgently need to increase investment in high-quality nature-based solutions. They currently receive a small proportion of existing climate-mitigation financing4,15, which does not reflect their potential.Carbon markets are increasingly relied on to finance nature-based solutions. But carbon offsets on the voluntary market are of variable quality. It can be unclear whether projects really represent a carbon sink, whether they are permanent or if they safeguard social and ecological factors. Offsets that adhere to standards can allow organizations to deliver lower-cost and hence larger climate-mitigation outcomes through nature-based solutions; however, budgets to emit fossil fuels should be ratcheted down rapidly to avoid delaying decarbonization from continued greenhouse-gas emissions.
Account for depreciation of natural capital
Nature-based solutions need both public and private finance; in particular, governments need to reward ecosystem stewardship while taxing polluters and ramping up regulation to ensure that companies meet strict social and environmental safeguards.The United Nations Framework Convention on Climate Change (UNFCCC) needs to provide clear guidelines on national-level accounting for nature-based solutions. This will guide the targets set in the Paris agreement’s Nationally Determined Contributions, and the monitoring, reporting and verification methodologies required to comply with these targets.The next UNFCCC meeting, COP26, is due to be held in Glasgow, UK, this November and provides an opportunity for national reporting systems to tighten national carbon accounting related to nature-based solutions. This would ensure that such solutions make a real, long-term contribution to carbon mitigation and could set metrics to ensure high biodiversity levels and maximize human well-being. One pressing issue for COP26 is Article 6 of the Paris agreement, which established a “mechanism to contribute to the mitigation of greenhouse gas emissions and support sustainable development”. A tightly regulated compliance market defined in Article 6 will provide the grounding for a tightly regulated voluntary offsetting market.COP26 also presents the chance to harmonize the goals of the UNFCCC and those of the Convention on Biological Diversity. For example, nature-based solutions projects are likely to be required to adhere to the principle of free prior informed consent of local people: local communities need to be involved at all stages of project planning and management. Similarly, nature-based solutions should be required to protect and enhance biodiversity. This work can build on existing social and biodiversity standards3.Our economy must be decarbonized at unprecedented rates to achieve net-zero targets by mid-century. Carbon must also be removed from the atmosphere to counter emissions that are hard to eliminate, using nature-based solutions and other means. To transform social and economic systems to deliver resilience in the face of ongoing climate impacts, the world must invest now in nature-based solutions that are ecologically sound, socially equitable and designed to pay dividends over a century or more. Properly managed, these could benefit many generations to come. More

Isolation and characterisation of viruses infecting the picoeukaryote Bathycoccus Clade BIIBathycoccus BII isolates RCC716 and RCC715 used in our experiments were originally cultured from a nutrient-limited region in the Indian Ocean. Clade BII as a whole has been reported extensively in warm oligotrophic ocean gyres based on metagenome analyses [22,23,24]. Peak abundances occurr when well-developed deep chlorophyll maxima are present, or throughout the photic zone during mixing periods at Station ALOHA of the Hawaii Ocean Time-series [12]. We targeted BATS for viral isolation in springtime because Bathycoccus has been observed at relatively high abundance during this period using qPCR [74]. Here, three viruses were isolated against RCC716 [12] using seawater flown from BATS/Bermuda to the laboratory, obviating bringing this finicky strain into the field for use as a viral host. We then purified the viruses by serial dilutions and sequenced the partial PolB gene to determine whether they were evolutionarily different from other cultured viruses. BLASTn and preliminary phylogenetic analysis using GenBank nr sequences indicated they were distinct from described viruses with deposited sequences, with best BLASTn hits to Bathycoccus prasinos viruses (62–74% nucleotide identity). Transmission electron microscopy (TEM) revealed that all three viruses have similar morphology to other characterised prasinoviruses [75], with icosahedral capsids diameter ranging between 120 and 140 nm (Fig. 1A).Fig. 1: Morphology and evolutionary relationships of newly discovered Bathycoccus viruses.A Transmission electron micrographs of BII-V1, BII-V2 and BII-V3 (scale bar, 50 nm). The capsid diameters (n = 6 virions) measured 138 ± 2 nm (BII-V1), 150 ± 5 nm (BII-V2) and 152 ± 11 nm (BII-V3). B Maximum Likelihood (ML) phylogenetic reconstruction of green algal viruses inferred from a concatenated alignment of 22 core proteins shared among the viruses (7,001 positions) under the LG + G + F model. Node support was calculated from 1000 bootstrap (BS) replicates, with all branches acquired support values of 100% (white dots). Viruses infecting Chlorella were used as an outgroup and the branch connecting the prasinoviruses to the outgroup was truncated for display purpose. The new Bathycoccus viruses isolated against Bathycoccus Clade II (sensu [12]) isolate RCC716 (named as species Bathycoccus calidus herein, see below) are in bold. Colours reflect different host species within each genus. Letters alongside vertical lines (a and b) correspond to Bathycoccus viral clades. C Venn diagram of the shared and unique protein-encoding genes in the genome sequences of the new Bathycoccus viruses.Full size imageGenomic sequencing and multi-gene evolutionary analysesAssembly of DNA sequences from the viral isolates after deep sequencing by Illumina rendered one complete dsDNA genome sequence (BII-V3), and two others may still be partial (Table 1). The BII-V2 genome, which was in one contig, was similar in size (~208 kb) to that of BII-V3 (~212 kb). The BII-V1 genome assembly was ~174 kb and comprised of four linear dsDNA scaffolds. The viral concentrate was deeply sequenced ( >50x coverage) and minor fragmentation of the genome was partially related to repeats that were not resolved during assembly. The total number of putative open reading frames (ORFs) varied from 220 in BII-V1 to 235 in BII-V2 (Table 1). Gene synteny was globally well-conserved across the BII-Vs and the BpV1 and BpV2 viruses of B. prasinos (Fig. S1), with limited genomic rearrangements. Other genome characteristics such as the coding proportion (~90%) and G + C % (~36%) were similar to other described prasinoviruses infecting Mamiellophyceae [64, 75], for which the reported number of proteins range from 203 to 268 and G + C % from 37 to 45%. However, the full-length PolB gene from the genome assemblies differed for BII-V3 from the other two, in having a 329 amino acid intein (Supplementary information table S3). Likewise, inteins have been reported at the same PolB position in uncultivated prasinoviruses from the subtropical Pacific Ocean [76], where Bathycoccus BII is abundant [12].Table 1 Genomic characteristics of the three Bathycoccus viruses (BII-Vs) isolated against Clade BII isolate RCC716.Full size tableTo reconstruct a robust phylogeny for the new viruses, we employed 22 proteins previously identified as being shared across all available green algal virus genomes, including both prasinoviruses and chloroviruses [65]. We found all 22 in the predicted coding sequences of BII-V1; however, DNA helicase (SNF2) was not found in BII-V2 or -V3, FAD-dependent thymidylate synthase (thy1) and the topoisomerase IV were not found in BII-V2, nor was the prolyl 4-hydroxylase in the BII-V3 genome. Additional searches with tBLASTn did not recover these genes or fragments of them, suggesting they have been lost. Phylogenomic reconstruction grouped the three BII-Vs with the two BpVs [32], in a fully supported clade that branched adjacent to a large group of viruses that infect various species of Ostreococcus and Micromonas (Fig. 1B). The clade of Bathycoccus viruses was segregated in two subclades with BII-V2 and BII-V3 clustering together adjacent to BII-V1 and BpVs (Fig. 1B). While better resolution of the position of BII-V1 awaits greater taxonomic sampling, our results demonstrated that the three new viruses branch adjacent or basally to BpVs.Variation in prasinovirus gene content and functions encodedThe three Bathycoccus Clade BII viruses had 72–77% of their proteins held in common, and ~30 unique proteins as well as a few proteins shared by just two of the three viruses (Fig. 1C). The 170 shared proteins had higher amino acid identities between BII-V2 and BII-V3 (73% aa identity) than to BII-V1 (69% and 68%, respectively). Generally, only 19–21% of Bathycoccus viral genes could be assigned a functional category, based on EggNOG classification. Similar functional category distributions were observed across both prasinoviruses and chloroviruses, including lipid metabolism, RNA processing and modification, and nucleotide metabolism and transport (Fig. 2A). Other functional categories were more variable, such as cell wall/membrane/envelope biogenesis genes prevalent in chloroviruses (potentially related to their enveloped nature), as well as genes involved in modification of the capsid with compounds such as with chitin and hyaluronan [77, 78] that are absent from prasinoviruses sequenced to date (Fig. 2A). Within prasinoviruses, most of the unique proteins in the Bathycoccus viruses lack defined functional categories. Among those with functional assignments, all five Bathycoccus viruses encoded a P2X receptor in the intracellular trafficking and secretion category, and both BII-V2 and -V3 encode two proteins putatively involved in degrading the aromatic compound 4-hydroxy-2-oxopentanoate to acetyl-CoA (secondary metabolite category), that otherwise are only encoded by one other prasinovirus, MpV1 [32]. Similar to the phylogenetic relationships, the functional category distributions of BII-V1 were closer to those of BpVs than to BII-Vs. The primary difference was in carbohydrate metabolism, where BII-V2 and -V3 each encodes ribulose-phosphate 3-epimerase (involved in the pentose phosphate pathway and carbon fixation; not found in any other available virus genomes, but encoded by B. prasinos) and TDP-glucose 4,6-dehydratase (involved in biosynthesis of rhamnose and encoded by most other chloroviruses and prasinoviruses [79]). Notably, the putative high-affinity phosphate transporter (PHO4, also termed HAPT) was only present in BII-V1 and BpV1, as well as OtV2 (isolated against the Ostreococcus Clade OII ecotype), and most sequenced viruses of O. lucimarinus (Supplementary information table S3). This gene is hypothesised to enhance phosphate uptake during infection under phosphorus‐limited host growth [25], as observed for the PstS phosphate transport system expressed by cyanophages [80], mitigating limitation of this key component of viral genomes. However, most isolated prasinovirus genomes come from waters that are not considered phosphate-limited, hence presence of this gene may connect to poising the host for responding to sudden availability of other nutrients, such as nitrogen, which is often limiting in the ecosystems from which these viruses were isolated. Studies of virus-cell responses under various limiting nutrients are required to understand the retention of this host-derived HGT.Fig. 2: Distribution of functions and orthologous protein families across genome-sequenced prasinoviruses.A Functional category distributions across 21 genome-sequenced prasinoviruses and chloroviruses based on EggNOG categorisation. Viruses are clustered by similarity in their distribution of the functional categories on the y-axis and the frequency of each category across the viral genomes determines clustering along the x-axis ordering. Genes with homology to proteins in the EggNOG database but could not be assigned a function are in the “function unknown” category. B Orthogroups presence/absence patterns ordered along the x-axis by ranking according to the total number of genes in the orthogroup. For inclusion, the orthogroup was required to include protein sequences from at least two different viral genomes. Viruses are ordered along the vertical by their presence/absence pattern reconstructed by hierarchical clustering (topology on the left). Top histogram: frequency of each orthogroup in sequenced prasinoviruses. C Genes in each virus (number) not assigned to any orthogroup, with viruses in the same vertical order as B.Full size imageHierarchical clustering of orthologous proteins revealed patterns across prasinoviruses that generally corresponded with phylogenetic relationships. The BII- and Bp-viruses shared 130 orthologous proteins and hierarchical clustering (Fig. 2B) followed the clade structure of the phylogenomic reconstruction (Fig. 1B) with the exception of BII-V1 that grouped with BII-Vs, as well as OtV6, which grouped with Micromonas viruses. These orthologous proteins had on average 72% amino acid identity between BII-V2 and BII-V3, and 88% between the two B. prasinos viruses, but between 65 to 67% when comparing members of these two groups (Table 2). BII-V1 orthologs also had 67% and 66% amino acid identity to BII-V2 and BII-V3, respectively, while they had 83% and 80% identity to BpV1 and BpV2, respectively. Collectively, these results indicate that BII-V2 and -V3 diverged from BpVs prior to the divergence of BII-V1.Table 2 Average percent amino acid identity of the orthologous proteins between the five Bathycoccus viruses.Full size tableOf the 130 orthologous Bathycoccus virus proteins, 37% were assigned putative functions revealing core components of this viral group (Supplementary information table S3). These included genes involved in DNA replication and transcription, including PolB (type II), a DNA topoisomerase, a transcription factor S-II, mRNA capping enzymes, ribonucleases, a ribonucleotide reductase, and a dUTPase. Several others are necessary for viral particle synthesis, such as genes encoding structural elements for assembling the virion, including capsid proteins (5–6 copies per genome), as well as transcriptional regulators connected to the replication cycle. The BII viruses showed a number of differences among orthologous protein families. In addition to each having “unique” protein sets, there was a set of BII-V specific orthogroups, as well as some shared with BpVs, and/or other prasinoviruses (Fig. 1C and Supplementary information table S3). First, six predicted proteins showed orthologs across the three BII-Vs, but were not present in other prasinoviruses sequenced to date. Only one of these six was assigned putative function, belonging to the XRE family of transcriptional regulators. Additionally, all BII viruses harboured a tandem duplication of the FstH gene, while other sequenced prasinoviruses (including the two Clade BI viruses) have one copy (Supplementary information table S3). This ATP-dependent metalloprotease has been shown to be involved in photosystem II repair in cyanobacteria [81], and is present in genomes of photosynthetic eukaryotes, including all Mamiellophyceae [15, 16]. In Arabidopsis and Chlamydomonas it has been shown to be involved in protein quality control in the thylakoid membranes [82]. A gene of unknown function was also duplicated in the BII-virus genomes, that is a single copy in BpVs and absent from other sequenced prasinoviruses. Genes putatively encoding a glucose-1-phosphate adenylyltransferase, a glycosyltransferase and a thiamine pyrophosphate-requiring enzyme involved in amino acid biosynthesis were sporadically found in BII-viruses.Considering the two Bathycoccus virus subclades (Fig. 1B), there is one predicted protein of unknown function exclusive to BpV1, BpV2 and BII-V1 and six predicted proteins shared only by BII-V2 and BII-V3. Among the latter, one belonged to the Ribulose-5-Phosphate-3-Epimerase (RPE) family, which catalyses the interconversion of D-ribulose 5-phosphate (Ru5P) into d-xylulose 5-phosphate, as part of the Calvin cycle (although no transit peptide was detected using TargetP) and in the oxidative pentose phosphate pathway. The ortholog analyses further showed that among prasinoviruses, 9, 17 and 18 genes were unique to BII-V1, BII-V2 and BII-V3, respectively (Fig. 2B). Apart from one nucleotidyltransferase and one glycosyltransferase (group 1) in BII-V1, none of these unique genes had known functions.To study the evolutionary aspects of the shared prasinovirus proteins, we constructed and examined 130 phylogenies of orthogroups shared between Bathycoccus viruses. Nine showed a topology where all three BII-Vs grouped together with full support (100% bootstrap support), separate from the BpV orthologs, and in contrast to the multi-gene phylogeny where BII-V1 grouped with BpVs (Fig. 1B). The average amino acid similarities within each of these nine protein ortholog groups ranged from 85 to 88% between BII-Vs proteins, while they were 77 to 81% between BII-Vs and BpVs, different from overall amino acid similarity averages (Table 2). Interestingly, proteins from three of these nine ortholog groups, all lacking known functions, were adjacent to each other in the genome, or separated by only one gene. This synteny and co-location likely reflects the acquisition of these genes before co-infecting viruses diverged via viral recombination [83].Infection dynamics of Bathycoccus virusesGeneral host specificity of BII-viruses was assessed using two B. prasinos isolates (CCMP1898 and RCC4222, Clade BI), the two available Clade BII isolates (RCC715 and RCC716), four Ostreococcus species and one Micromonas species (Table 3). None were able to infect the B. prasinos, Ostreococcus or Micromonas isolates tested, suggesting BII-V specificity for Bathycoccus Clade BII. Similar host specificity has been observed in O. lucimarinus viruses, none of which infect O. tauri [64], and other viruses of eukaryotic and prokaryotic algae [84, 85]. Some other prasinoviruses appear to have broader host ranges [85,86,87], or their host species are less divergent than the two known Bathycoccus clades. For example, generalist viruses isolated against Micromonas commoda can infect M. bravo [85]. Further investigations are necessary to determine the extent to which the six shared proteins in BII-Vs (absent from BpVs), are responsible for the differences in host and virus specificity of interactions, versus variations in the shared Bathycoccus virus proteins (65–83% similarity). Importantly, host specificity tests for the new viruses described herein were limited by weak sampling of Bathycoccus diversity (in culture; all that we could acquire were tested).Table 3 Results of cross infectivity tests of BII-V1, BII-V2 and BII-V3 against isolates representing various picoprasinophyte species within the Class Mamiellophyceae.Full size tableAlthough specific for the BII clade, the three BII-Vs exhibited variations in infectivity of the two cultured BII strains available, despite their isolation from the same sample and having identical ITS1 and ITS2 sequences. BII-V1 lysed and cleared RCC715 and RCC716 cultures after four days (Table 3). The same was true for BII-V2 and BII-V3, when incubated with RCC716. Different from results for BII-V1, we found that while BII-V2 and -V3 initially lysed RCC715 cultures, resistant populations became evident at day 7 of infectivity tests, and measureable lysis of RCC715 could not be achieved thereafter. These results underscored the need to further examine host-virus interactions for the three new viruses.Infection dynamics over time course experiments further illuminated differences in BII-V impacts on hosts. In these experiments, growth rates of the uninfected (control) RCC715 and RCC716 cultures were 0.45 ± 0.04 day−1 and 0.49 ± 0.06 per day, respectively, similar to rates during the pre-experiment acclimation period (T-test, p  > 0.05). Host and virus dynamics were similar for RCC715 and RCC716 infected with BII-V1 (Fig. S2 and Fig. 3), with cell numbers starting to diverge from control abundances 10 h after inoculation (T-test, p  More

Data setsRadioactivity measurement data for several species of wild game mammals and birds in Fukushima Prefecture from May 2011 to March 2018 were released to the public by the Fukushima Prefecture Government (https://emdb.jaea.go.jp/emdb/en/portals/1040501000/). We extracted the data for wild boar (Sus scrofa), 1404 samples, and Asian black bear (Ursus thibetanus), 422 samples. The resulting boar and bear data sets contained total radioactive cesium activity (134Cs + 137Cs isotopes) values (in Bq/kg) from animals captured at different times and locations within Fukushima Prefecture. The data were imported for analysis into R 4.0.3 software21.We ln-transformed the cesium activity values to bring their distribution closer to normal, creating the variable LnCsTot. To facilitate regression analyses (described below), we removed instances of missing data and cesium levels below detection: 20 samples (1.4%) for boar and 15 samples (3.3%) for bears. The time when each sample was taken (labeled “Day of collection” in the Fukushima Prefecture Government data set) was converted to years since the Fukushima accident (since March 11, 2011), assuming that 1 year = 365.25 days. This time of sample collection in years was called variable T.Since for each sample some time passed between sample collection and radioactivity measurement (labeled “Result found Date”, called Tr in our notation), we needed to correct the reported LnCsTot values for physical decay over this time, which was different for different samples. The procedure used to perform this correction is described in Supplementary methods. The data with corrected total cesium values (LnCsc) are provided in Supplementary data (Supplementary_Dataset_File_Full).Mathematical modelTo describe the data on ln-transformed total radioactive cesium levels (LnScc) in each species as function of time after the accident (T), we developed the following simple mathematical model (Eqs. 1A, 1B):$${LnCs}_{c}=X+Q-mu times {T}^{nu }+Atimes mathrm{sin}left[2times pi times left(T+Pright)right], $$
(1A)
$$X=mathrm{ln}left[mathrm{exp}left(LnCs{134}_{t{0}_{r}}right)times {2}^{-frac{T}{{Th}_{Cs134}}}+mathrm{exp}left(LnCs{137}_{t{0}_{r}}right)times {2}^{-frac{T}{{Th}_{Cs137}}}right]$$
(1B)
Here the term X represents the estimated average radioactive cesium level in the studied area, based on the intercepts (LnCs134t0r for 134Cs and LnCs137t0r for 137Cs, respectively) from robust regression discussed in Supplementary methods, and taking into account only physical decay for each isotope (with half-lives of ThCs134 for 134Cs and ThCs134 for 137Cs, respectively). The terms Q, µ, ν, A and P represent adjustable model parameters. Parameter Q represents the fitted relationship between radioactive cesium levels in the animal (Bq/kg), relative to the external environment (Bq/m2). Parameter µ represents the net rate of radioactive cesium reduction in animal tissues over time due to all processes except physical decay (e.g. decrease in bioavailability due to migration of cesium into deeper soil layers, human-mediated cleanup efforts, etc.). Parameter ν is a potential power dependence for these processes. By default, ν was set to ν = 1, but exploratory calculations using ν = 2 or treating ν as a freely adjustable parameter (≥ 0.1) were performed as well. Parameters A and P in the sine function represent a sinusoidal approximation for seasonal changes in radioactive cesium levels in animal tissues (e.g. due to seasonal variations in diet and life style), where A is the amplitude of the oscillations, P is the phase shift, and the period is set to 1 year. For simplicity, these parameters were assumed to be the same for both studied cesium isotopes. The descriptions of each parameter are also presented in Table 1.Table 1 The meanings of all parameters used in our mathematical model (Eq. 1A, 1B) for radioactive cesium levels in wild boar (Sus scrofa) and Asian black bear (Ursus thibetanus).Full size tableModel fitting approachesInitially, we used nonlinear ordinary least squares (OLS) regression (nls R function) to fit the model (Eq. 1A, 1B) to the data. To find the global optimum fit, we repeated the fitting procedure 2000 times with slightly different random initial parameter values and recorded the solution with the smallest root mean squared error (RMSE). Diagnostics on this regression included checking of convergence criteria and analyses of residuals (by scatter plot and histogram, regressing residuals as function of T, visualizing the QQ plot, autocorrelation and partial autocorrelation functions with 95% confidence intervals, performing the Shapiro–Wilk normality test, and calculating skewness and kurtosis). For boar data, diagnostics revealed problems with convergence (both X-convergence and relative convergence) and non-normality of residuals: e.g. Shapiro–Wilk p-value = 1.476 × 10–7, skewness = − 0.37, kurtosis = 3.50. For black bear data similar problems occurred with convergence, but residuals were closer to the normal distribution (perhaps due to smaller sample size): e.g. Shapiro–Wilk p-value = 0.0526, skewness = − 0.058, kurtosis = 2.45.Due to these issues, we used robust nonlinear regression (nlrob R package) to reduce the effects of “outlier” data points. To find the global optimum, we repeated the fitting procedure 2000 times with slightly different random initial parameter values and selected the solution with the smallest absolute value of median residuals. The best-fit parameters for OLS and robust regressions were somewhat different for both boar and bear data. For boar data, the minimum robustness weight was 0.339 and the median was 0.762, and the corresponding values for black bear data were 0.557 and 0.821, respectively.For each species, we compared the performances of model variants with different assumptions about parameter ν: (1) The default case with ν = 1, which represents an exponential rate of radioactive cesium decrease due to processes other than physical decay. (2) The case with ν = 2, which represents quadratic decay. (3) The case with ν being freely adjustable (≥ 0.1). The comparisons were based on Akaike information criterion (AIC)22,23. The purpose of these calculations was to better assess the shape of the time course for non-physical factors involved in radioactive cesium level decline in animal tissues over time after the accident.In addition to analyzing the full data set for each species, we also performed separate analyses on subsets of data from specific locations: from those districts of Fukushima Prefecture where the mean radioactive cesium levels in animal samples were the highest, and where a sufficiently large number of samples was present. For wild boar, the two selected districts for this subset analysis were Soso and Kenpoku (819 samples), and for black bear they were Kenpoku and Kenchu (163 samples).To further assess the sensitivity of model results to geographical and temporal factors, we also constructed a separate subset of data for each species. This subset excluded data from the Aizu and Minamiaizu districts, which are separated by mountains from the Fukushima Daiichi Nuclear Power Plant, and excluded data collected ≤ 6 months after the accident. These restrictions were intended to determine model performance on data collected in a more geographically contiguous area after the initial abrupt changes in contamination levels were completed and the system entered the phase of more stable kinetics. The purpose of all these analyses was to assess whether the time course of radioactive cesium levels in the bodies of each species differed between locations with high contamination vs. those with lower contamination, and as function of time after the accident.We were interested in quantifying not only the center of the distribution of radioactive cesium values in each species over time, but also in assessing the lower and upper tails of this distribution. For this purpose, we fitted the model (Eq. 1A, 1B) for each species using quantile regression (nlrq function in quantreg R package) for the median (50th percentile), and also for the 25th and 75th percentiles. Initial parameter estimates for the quantile regressions were taken from best-fit parameters from robust regression described above. The 25th and 75th percentiles were selected instead of more extreme values (e.g. 5th and 95th) because the latter resulted in poor quality fits due to limited amounts of data at the fringes of the distribution.To assess the variability of model parameters by location in more detail, we used mixed effects modeling (nlme R package) on the data from each species. Since original OLS fits suggested substantial deviations of residuals from the normality assumption, we performed mixed effects modeling on data with some outlier data points removed. The OutlierDetection package in R removed 43 boar samples and 5 bear samples. These outliers are listed in the Supplementary_outlier_data_points file. The remaining samples were used for mixed effects model fitting, but model performance metrics like coefficient of determination (R2) and RMSE were assessed on the full data set (with outliers included) for each species.Since the Fligner-Killeen test of homogeneity of variances by district generated low p-values for both species (4.6 × 10–14 for boar and 0.018 for black bear), we allowed modelled variances to differ by district (using the weights option in nlme). We investigated several random effects structures for some or all model parameters, with randomness by district only, or by district and municipality within district. Model diagnostics were the same as for fixed effects OLS modeling described above, and also included boxplots of model residuals by district. The mixed effects model variants with different random effects structures were compared using the anova function in R, and also by assessing convergence criteria, normality of residuals, skewness, and kurtosis. Consequently, preferred mixed effects model variants were selected for the full data as well as for the subset of two districts with high radioactive cesium levels, separately for each species.Model extrapolation from training to testing dataTo investigate how the robust and quantile regression fits of our model could extrapolate beyond the time range that was used for model fitting, we split the data for each species into “training” (early times) and “testing” (later times) parts. The split was done based on time since the accident (T variable), so that approximately ½ of the samples were assigned to the training and testing sets, respectively. For wild boar data, the training set included times between 0.20 and 3.45 years after the accident, and the testing set included times between 3.45 and 7.03 years. For black bear data, the training set included times between 0.42 and 3.46 years after the accident, and the testing set included times between 3.46 and 6.87 years.We also evaluated an alternative approach to splitting the data, where the split was done randomly instead of by time. In other words, any data point regardless of time had an equal probability of being assigned to either the training or the testing data set. Both the training and testing data subsets generated by this random split included the complete time range. This approach was implemented in context of the sensitivity analysis described above.For each species, robust and quantile regressions were fitted to training data, and their predictions were calculated for testing data. For robust regression, RMSE was calculated on testing data for two scenarios: (1) for the model fitted to training data only, and (2) for the model fitted over the entire data range (training + testing). These RMSE values for conditions 1 and 2 were compared to assess the quality of model extrapolation. Extrapolation performance for robust and quantile regressions was also assessed visually by plotting the model predictions and data.Application of the model to wild boar data from the Chernobyl accident areaTo compare the results of our analysis of wild boar contamination with radioactive cesium in the area affected by the Fukushima accident with data from another location, we also analyzed wild boar data from the Chernobyl accident area. These data were published by Gulakov14 and contain summaries of 137Cs contamination levels in the muscles of 188 boar collected between 1991 and 2008 (i.e. from 5 to 22 years after the 1986 accident). Sampling was carried out in three zones with different land contamination levels with 137Cs. This data set provides important information on radioactive cesium contamination in wild boar in the Chernobyl area. Unfortunately, 137Cs measurements in each sampled boar were not provided by Gulakov14, and only summary statistics are available for each zone and year after the accident (Tables 1–3 in reference14): number of animals, mean, minimum and maximum 137Cs levels.We could not apply the full model (Eq. 1A, 1B) to these summary data which lacked seasonality information and 134Cs data. However, we were able to perform a weighted linear regression to quantify the ecological half-life of 137Cs in Chernobyl boar and the relationship between 137Cs levels in the animals (Bq/kg), relative to the external environment (Bq/m2). The data used for this analysis, derived from Gulakov14, are provided in Supplementary data (Supplementary_Dataset_File_Full). They contain the following variables. Zone = location of sample collection (Alienation, Permanent control or Periodic control). Time = time in years after the Chernobyl accident. LnMeanCs = ln-transformed mean 137Cs level in boar muscle (Bq/kg). LnMeanCs_c = LnMeanCs − X, where X is ln-transformed 137Cs land contamination (Bq/m2) in the given zone, corrected for physical decay of 137Cs. Weight = weighting of each data point used for regression. Weight = number of animals/(ln[maximum 137Cs level] − ln[minimum 137Cs level])2. These approximately inverse-variance weights were normalized by the overall mean across all data points, so that the mean weight across all data points was set to 1.These data were analyzed by weighted linear regression in R, where LnMeanCs_c was allowed to depend on Time and Zone variables. Model variants containing all possible combinations and pairwise interactions between these predictor variables were fitted and their performances were compared using the Akaike information criterion with correction for small sample size (AICc). These calculations were performed using the glmulti R package. Multimodel inference (MMI) was performed on this collection of fitted model variants. It resulted in the calculation of model-averaged parameter estimates, 95% CIs and importance scores, corrected for model selection uncertainty. Of main interest here were the intercept parameter, which is analogous to parameter Q in the full model (Eq. 1A, 1B), and the Time parameter, which is analogous to parameter µ in the full model. The ecological half-life for 137Cs was calculated based on the Time parameter. More

1.Dong Y, Silbermann M, Speiser A, Forieri I, Linster E, Poschet G, et al. Sulfur availability regulates plant growth via glucose-TOR signaling. Nat Commun. 2017;8:1174.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
2.Freney JR, Melville GE, Williams CH. Soil organic matter fractions as sources of plant-available sulphur. Soil Biol Biochem. 1975;7:217–21.CAS 
Article 

Google Scholar 
3.Kopittke PM, Dalal RC, Finn D, Menzies NW. Global changes in soil stocks of carbon, nitrogen, phosphorus, and sulphur as influenced by long‐term agricultural production. Global Change Biol. 2017;23:2509–19.Article 

Google Scholar 
4.Ciaffi M, Paolacci AR, Celletti S, Catarcione G, Kopriva S, Astolfi S. Transcriptional and physiological changes in the S assimilation pathway due to single or combined S and Fe deprivation in durum wheat (Triticum durum L.) seedlings. J Exp Bot. 2013;64:1663.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
5.Ma Q, Luo Y, Wen Y, Hill PW, Chadwick DR, Wu L, et al. Carbon and sulphur tracing from soil organic sulphur in plants and soil microorganisms. Soil Biol Biochem. 2020;150:107971.CAS 
Article 

Google Scholar 
6.Piotrowska-Długosz A, Siwik-Ziomek A, Długosz J, Gozdowski D. Spatio-temporal variability of soil sulfur content and arylsulfatase activity at a conventionally managed arable field. Geoderma. 2017;295:107–18.Article 
CAS 

Google Scholar 
7.Fitzgerald JW, Watwood ME. Amino-acid metabolism in forest soil—Isolation and turnover of organic matter covalently labelled with 35S-methionine. Soil Biol Biochem. 1988;20:833–8.CAS 
Article 

Google Scholar 
8.Ma Q, Wen Y, Pan W, Macdonald A, Hill PW, Chadwick DR, et al. Soil carbon, nitrogen, and sulphur status affects the metabolism of organic S but not its uptake by microorganisms. Soil Biol Biochem. 2020;149:107943.CAS 
Article 

Google Scholar 
9.Jan MT, Roberts P, Tonheim SK, Jones DL. Protein breakdown represents a major bottleneck in nitrogen cycling in grassland soils. Soil Biol Biochem. 2009;41:2272–82.CAS 
Article 

Google Scholar 
10.Farrell M, Macdonald LM, Hill PW, Wanniarachchi SD, Farrar J, Bardgett RD, et al. Amino acid dynamics across a grassland altitudinal gradient. Soil Biol Biochem. 2014;76:179–82.CAS 
Article 

Google Scholar 
11.Wilkinson A, Hill PW, Farrar JF, Jones DL, Bardgett RD. Rapid microbial uptake and mineralization of amino acids and peptides along a grassland productivity gradient. Soil Biol Biochem. 2014;72:75–83.CAS 
Article 

Google Scholar 
12.Hill PW, Jones DL. Plant–microbe competition: does injection of isotopes of C and N into the rhizosphere effectively characterise plant use of soil N? N Phytol. 2018;221:796–806.Article 
CAS 

Google Scholar 
13.Godwin CM, Cotner JB. Aquatic heterotrophic bacteria have highly flexible phosphorus content and biomass stoichiometry. ISME J. 2015;9:2324–27.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
14.Hartman WH, Ye R, Horwath WR, Tringe SG. A genomic perspective on stoichiometric regulation of soil carbon cycling. ISME J. 2017;11:2652–65.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
15.Manzoni S, Jackson RB, Trofymow JA, Porporato A. The global stoichiometry of litter nitrogen mineralization. Science. 2008;321:684–6.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
16.Cui J, Zhu Z, Xu X, Liu S, Jones DL, Kuzyakov Y, et al. Carbon and nitrogen recycling from microbial necromass to cope with C:N stoichiometric imbalance by priming. Soil Biol Biochem. 2020;142:107720.CAS 
Article 

Google Scholar 
17.Wei X, Zhu Z, Liu Y, Luo Y, Deng Y, Xu X, et al. C:N:P stoichiometry regulates soil organic carbon mineralization and concomitant shifts in microbial community composition in paddy soil. Biol Fert Soils. 2020;56:1093–107.Article 
CAS 

Google Scholar 
18.Maria M, Wolfgang W, Sophie ZB, Andreas R. Stoichiometric imbalances between terrestrial decomposer communities and their resources: mechanisms and implications of microbial adaptations to their resources. Front Microbiol. 2014;5:22.
Google Scholar 
19.Qiao N, Xu XL, Hu YH, et al. Carbon and nitrogen additions induce distinct priming effects along an organic-matter decay continuum. Sci Rep-UK. 2016;6:19865.CAS 
Article 

Google Scholar 
20.Kirkby CA, Kirkegaard JA, Richardson AE, Wade LJ, Blanchard C, Batten G. Stable soil organic matter: a comparison of C:N:P:S ratios in Australian and other world soils. Geoderma. 2011;163:197–208.CAS 
Article 

Google Scholar 
21.Manzoni S, Čapek P, Mooshammer M, Lindahl BD, Richter A, Šantrůčková H. Optimal metabolic regulation along resource stoichiometry gradients. Ecol Lett. 2017;20:1182–91.PubMed 
Article 
PubMed Central 

Google Scholar 
22.Mooshammer M, Wanek W, Hämmerle I, Fuchslueger L, Hofhansl F, Knoltsch A, et al. Adjustment of microbial nitrogen use efficiency to carbon:nitrogen imbalances regulates soil nitrogen cycling. Nat Commun. 2014;5:3694.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
23.Delgado-Baquerizo M, Reich PB, Khachane AN, Campbell CD, Thomas N, Freitag TE, et al. It is elemental: soil nutrient stoichiometry drives bacterial diversity. Environ Microbiol. 2017;19:1176–88.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
24.Mayor JR, Mack MC, Schuur EAG. Decoupled stoichiometric, isotopic, and fungal responses of an ectomycorrhizal black spruce forest to nitrogen and phosphorus additions. Soil Biol Biochem. 2015;88:247–56.CAS 
Article 

Google Scholar 
25.Mooshammer M, Hofhansl F, Frank AH, Wanek W, Hämmerle I, Leitner S, et al. Decoupling of microbial carbon, nitrogen, and phosphorus cycling in response to extreme temperature events. Sci Adv. 2017;3:e1602781.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
26.Chen J, Seven J, Zilla T, Dippold MA, Blagodatskaya E, Kuzyakov Y. Microbial C:N:P stoichiometry and turnover depend on nutrients availability in soil: A 14C, 15N and 33P triple labelling study. Soil Biol Biochem. 2019;131:206–16.CAS 
Article 

Google Scholar 
27.Ma Q, Wu L, Wang J, Ma J, Zheng N, Hill PW, et al. Fertilizer regime changes the competitive uptake of organic nitrogen by wheat and soil microorganisms: An in-situ uptake test using 13C, 15N labelling, and 13C-PLFA analysis. Soil Biol Biochem. 2018;125:319–27.CAS 
Article 

Google Scholar 
28.Jiang G, Zhang W, Xu M, Kuzyakov Y, Zhang X, Wang J, et al. Manure and Mineral Fertilizer Effects on Crop Yield and Soil Carbon Sequestration: a Meta‐Analysis and Modeling Across China. Glob Biogeochem Cy. 2018;32:1659–72.CAS 
Article 

Google Scholar 
29.Liu S, Wang J, Pu S, Blagodatskaya E, Kuzyakov Y, Razavi BS. Impact of manure on soil biochemical properties: a global synthesis. Sci Total Environ. 2020;745:141003.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
30.Cheng L, Zhang N, Yuan M, Xiao J, Qin Y, Deng Y, et al. Warming enhances old organic carbon decomposition through altering functional microbial communities. ISME J. 2017;11:1825–35.PubMed 
PubMed Central 
Article 

Google Scholar 
31.Ma Q, Wen Y, Wang D, Sun X, Hill PW, Macdonald A, et al. Farmyard manure applications stimulate soil carbon and nitrogen cycling by boosting microbial biomass rather than changing its community composition. Soil Biol Biochem. 2020;144:107760.CAS 
Article 

Google Scholar 
32.Ma Q, Wen Y, Ma J, Macdonald A, Hill PW, Chadwick DR, et al. Long-term farmyard manure application affects soil organic phosphorus cycling: a combined metagenomic and 33P/14C labelling study. Soil Biol Biochem. 2020;149:107959.CAS 
Article 

Google Scholar 
33.Jones DL, Hill PW, Smith AR, Farrell M, Ge T, Banning NC, et al. Role of substrate supply on microbial carbon use efficiency and its role in interpreting soil microbial community-level physiological profiles (CLPP). Soil Biol Biochem. 2018;123:1–6.CAS 
Article 

Google Scholar 
34.Jones DL, Magthab EA, Gleeson DB, Hill PW, Sánchez-Rodríguez AR, Roberts P, et al. Microbial competition for nitrogen and carbon is as intense in the subsoil as in the topsoil. Soil Biol Biochem. 2018;117:72–82.CAS 
Article 

Google Scholar 
35.Mariano E, Jones DL, Hill PW, Trivelin PCO. Mineralisation and sorption of dissolved organic nitrogen compounds in litter and soil from sugarcane fields. Soil Biol Biochem. 2016;103:522–32.CAS 
Article 

Google Scholar 
36.Corre M, Brumme R, Veldkamp EF. Changes in nitrogen cycling and retention processes in soils under spruce forests along a nitrogen enrichment gradient in Germany. Glob Change Biol. 2010;13:1509–27.Article 

Google Scholar 
37.Spohn M, Kuzyakov Y. Phosphorus mineralization can be driven by microbial need for carbon. Soil Biol Biochem. 2013;61:69–75.CAS 
Article 

Google Scholar 
38.Wu J, Joergensen RG, Pommerening B, Chaussod R, Brookes PC. Measurement of soil microbial biomass C by fumigation-extraction—an automated procedure. Soil Biol Biochem. 1990;22:1167–9.CAS 
Article 

Google Scholar 
39.Glanville H, Hill PW, Schnepf A, Oburger E, Jones DL. Combined use of empirical data and mathematical modelling to better estimate the microbial turnover of isotopically labelled carbon substrates in soil. Soil Biol Biochem. 2015;94:154–68.Article 
CAS 

Google Scholar 
40.Greenfield LM, Hill PW, Paterson E, Baggs EM, Jones DL. Methodological bias associated with soluble protein recovery from soil. Sci Rep-UK. 2018;8:11186.Article 
CAS 

Google Scholar 
41.Näsholm T, Kielland K, Ganeteg U. Uptake of organic nitrogen by plants. N Phytol. 2010;182:31–48.Article 
CAS 

Google Scholar 
42.Farrell M, Prendergast-Miller M, Jones DL, Hill PW, Condron LM. Soil microbial organic nitrogen uptake is regulated by carbon availability. Soil Biol Biochem. 2014;77:261–7.CAS 
Article 

Google Scholar 
43.Niknahad-Gharmakher H, Piutti S, Machet JM, Benizri E, Recous S. Mineralization-immobilization of sulphur in a soil during decomposition of plant residues of varied chemical composition and S content. Plant Soil. 2012;360:391–404.CAS 
Article 

Google Scholar 
44.Vong PC, Piutti S, Slezackdeschaumes S, Benizri E, Guckert A. Effects of low-molecular-weight organic compounds on sulphur immobilization and re-mineralization and extraction of immobilized sulphur by hot-water and acid hydrolysis. Eur J Soil Sci. 2010;61:287–97.CAS 
Article 

Google Scholar 
45.Takriti M, Wild B, Schnecker J, Mooshammer M, Knoltsch A, Lashchinskiy N, et al. Soil organic matter quality exerts a stronger control than stoichiometry on microbial substrate use efficiency along a latitudinal transect. Soil Biol Biochem. 2018;121:212–20.CAS 
Article 

Google Scholar 
46.Zhu Z, Ge T, Luo Y, Liu S, Xu X, Tong C, et al. Microbial stoichiometric flexibility regulates rice straw mineralization and its priming effect in paddy soil. Soil Biol Biochem. 2018;121:67–76.CAS 
Article 

Google Scholar 
47.Xu X, Hui D, King AW, Song X, Thornton PE, Zhang L. Convergence of microbial assimilations of soil carbon, nitrogen, phosphorus and sulfur in terrestrial ecosystems. Sci Rep-UK. 2015;5:17445.CAS 
Article 

Google Scholar 
48.Fanin N, Fromin N, Buatois B, Ttenschwiler SH. An experimental test of the hypothesis of non-homeostatic consumer stoichiometry in a plant litter-microbe system. Ecol Lett. 2013;16:764–72.PubMed 
Article 
PubMed Central 

Google Scholar 
49.Peters MK, Hemp A, Appelhans T, Becker JN, Behler C, Classen A, et al. Climate–land-use interactions shape tropical mountain biodiversity and ecosystem functions. Nature. 2019;568:88–92.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
50.Kaiser C, Franklin O, Dieckmann U, Richter A. Microbial community dynamics alleviate stoichiometric constraints during litter decay. Ecol Lett. 2014;17:680–90.PubMed 
PubMed Central 
Article 

Google Scholar 
51.Xu X, Schimel JP, Thornton PE, Song X, Yuan F, Goswami S. Substrate and environmental controls on microbial assimilation of soil organic carbon: a framework for Earth system models. Ecol Lett. 2014;17:547–55.PubMed 
Article 
PubMed Central 

Google Scholar 
52.Apostel C, Dippold M, Kuzyakov Y. Biochemistry of hexose and pentose transformations in soil analyzed by position-specific labeling and 13C-PLFA. Soil Biol Biochem. 2015;80:199–208.CAS 
Article 

Google Scholar 
53.Manzoni S, Taylor P, Richter A, Porporato A, Ågren GI. Environmental and stoichiometric controls on microbial carbon‐use efficiency in soils. N Phytol. 2012;196:79–91.CAS 
Article 

Google Scholar 
54.Fitzgerald JW, Hale DD, Swank WT. Sulphur-containing amino acid metabolism in surface horizons of a hardwood forest. Soil Biol Biochem. 1988;20:825–31.CAS 
Article 

Google Scholar 
55.Akashi H, Gojobori T. Metabolic efficiency and amino acid composition in the proteomes of Escherichia coli and Bacillus subtilis. P Natl Acad Sci USA. 2002;99:3695–700.CAS 
Article 

Google Scholar 
56.Nozaki T, Ali V, Tokoro M. Sulfur-Containing Amino Acid Metabolism in Parasitic Protozoa. Adv Parasit. 2005;60:1–99.Article 

Google Scholar 
57.Takagi H, Ohtsu I. l-Cysteine Metabolism and Fermentation in Microorganisms. Adv Biochem Eng Biotechnol. 2016;159:129–51.
Google Scholar 
58.Bustos I, Miguel AM, Fouad A, Carmen P, Teresa R, CM M. Volatile sulphur compounds-forming abilities of lactic acid bacteria: C-S lyase activities. Int J Food Microbiol. 2011;148:121–7.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
59.Assefa MK, Tucher SV, Schmidhalter U. Soil sulfur availability due to mineralization: soil amended with biogas residues. J Soil Sci Enviro Manag. 2014;5:13–9.CAS 
Article 

Google Scholar  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

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

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

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