Synergistic Effects of Viruses and Environmental Gradients on Carbon Cycling in a River Ecosystem

TL;DR

This study shows how viruses in rivers influence carbon cycling, especially in agricultural areas, by responding to land use and nitrogen levels.

Contribution

The study reveals viruses as key environmental signal transducers in river carbon cycling, integrating viral ecology into biogeochemical models.

Findings

Viral activity peaks in agricultural plains, redirecting host carbon metabolism.

Cropland coverage increases TCA cycle transcription by 1.8-fold via viral gene expression.

Nitrogen loading suppresses both viral activity and carbon fixation in river ecosystems.

Abstract

Rivers play a pivotal yet underexplored role in the global carbon cycle, particularly regarding viral regulation of eukaryotic carbon processing across landscapes. Through meta-transcriptomic profiling along the Yongding River, we uncovered significant spatial disparities, with viral genes—notably those encoding major capsid proteins of large DNA viruses—and eukaryotic genes involved in carbon fixation, conversion, and metabolism exhibiting peak activity in agricultural plains. Integrated analyses revealed that land use intensifies viral activity, which subsequently redirects host carbon metabolism. Structural equation modeling demonstrated that cropland coverage elevates viral expression, correlating with a 1.8-fold increase in TCA cycle transcription, while nitrogen loading suppresses both viral activity and carbon fixation. Phylogenetic evidence supports virus–host specificity as a mechanistic driver. These findings position viruses as critical environmental signal transducers that shape riverine carbon cycling, underscoring the necessity of incorporating viral ecology into predictive biogeochemical models under global change.

Riverine ecosystems represent critical nodes in the global carbon cycle, where the mechanistic role of viruses in modulating eukaryotic carbon cycling remains underexplored, particularly across heterogeneous landscapes. Here, we applied metatranscriptomics to dissect how multi-scale environmental factors and viral gene activity jointly regulate the spatial transcription of carbon cycling genes in riverine eukaryotic communities along the Yongding River, China. Our analyses reveal pronounced spatial heterogeneity in both viral gene expression—notably major capsid proteins of large eukaryotic DNA viruses—and carbon fixation, conversion, and metabolism pathways, peaking in agriculturally impacted plain regions. Multivariate statistics and network analyses demonstrate that land use enhances viral gene activity, serving as biological amplifiers that modulate host carbon metabolism and transformation. Structural equation modeling further identifies a cascade in which cropland coverage elevates viral gene expression, ultimately driving a 1.8-fold increase in TCA cycle gene transcription in plain regions, whereas nitrogen loading at the site scale suppresses viral activity and carbon fixation. Phylogenetic analysis corroborates that virus–host specificity underpins these spatial patterns. Collectively, these findings advance a new model in which viruses act as key intermediaries, transmitting multiscale environmental signals to shape riverine carbon cycling. Our study highlights the urgency of incorporating viral ecology into predictive frameworks of riverine biogeochemical cycling under accelerating environmental change.

Figures

Full text

1. Introduction

River ecosystems are vital hotspots in the global carbon cycle, mediating terrestrial–marine exchange and contributing approximately 50% of organic carbon exports from freshwater systems [1]. Riverine carbon cycling plays a key role in regional climate feedbacks and shapes the productivity of downstream and coastal ecosystems through the transport of dissolved organic carbon (DOC) and particulate organic carbon (POC) [2]. While traditional research on carbon fluxes has predominantly focused on physicochemical processes (e.g., hydrologic disturbance, redox conditions), the critical roles of microorganisms and viruses have been increasingly recognized as the “hidden engine” of carbon cycling [3].

Eukaryotes, including algae, fungi, protozoans, and meiofauna, play dual roles in riverine organic carbon dynamics, serving as both primary producers and decomposers. Planktonic and benthic algae account for a majority (60–80%) of autochthonous organic carbon via photosynthetic fixation, a process strongly regulated by light availability, nutrient abundance and stoichiometry (e.g., N/P ratio), and water temperature [4]. Fungi drive the decomposition of recalcitrant terrestrial organic matter through secreting extracellular enzymes (e.g., cellulases, lignin peroxidases), thereby enhancing the bioavailability of organic matter [5]. Protozoa (e.g., ciliates, flagellates) can accelerate carbon remineralization and energy transfer within the microbial loop by grazing on bacteria and microalgae [6]. Studies indicate that microorganisms can influence ecosystem carbon cycling through strategies such as photobeterotrophy and mixotrophy to acquire energy and carbon. For instance, in oligotrophic glacial foreland soils, approximately 8% of the microbial genomes are predicted to engage in photobeterotrophic growth using rhodopsins or photosystem II, serving as a metabolic supplement in variable environments [7]. This metabolic flexibility allows micro-eukaryotes to adapt to resource fluctuations, directly impacting primary production and carbon transformation processes. Additionally, metazoans (e.g., rotifers, cladocerans) influence carbon sedimentation and resuspension via vertical migration and bioturbation.

Recent research demonstrates that species interactions—such as diatom–bacteria consortia—can directly modulate rbcL gene expression and thereby influence carbon fixation [3]. In addition, fungal communities drive diel fluctuations in dissolved organic carbon (DOC) production through the rhythmic expression of lignin-degrading genes, such as laccases [4]. However, most current studies focus on isolated taxa or environmental gradients, leaving a critical gap in our understanding of how eukaryotic metabolic functions involved in carbon cycling respond to multi-scale spatiotemporal heterogeneity. Our previous findings demonstrated that environmental variables across basin, reach, and site scales collectively shape the spatial variation in eukaryotic functional gene expression in rivers [8]. Yet, the regulatory mechanisms linking these gradients to carbon cycling genes remain unclear. Viruses, increasingly recognized as major modulators of host metabolism and biogeochemical cycling, impact host carbon pathways through metabolic reprogramming, cell lysis, and auxiliary metabolic genes (AMGs) [3,9,10].

Despite their importance, it is not yet well understood how viral gene activity, together with multi-scale environmental drivers, regulates eukaryotic carbon cycling in rivers [11]. Metatranscriptomics now enables simultaneous profiling of viral and eukaryotic gene expression across field environmental gradients, but field-based evidence for coupled viral and host carbon gene regulation under multiple environmental drivers remains scarce [12]. Recent multi-omics investigations of deep-sea ecosystems demonstrate the power of integrating metagenomic, viromic, and metatranscriptomic data to directly link the composition, life strategy, and functional gene expression of viral communities to the metabolism of complex organic matter under distinct geochemical conditions, providing a methodological framework for similar studies in riverine environments [13].

Here, we leverage metatranscriptomic data from the Yongding River’s mountain–plain–coast gradient to resolve how viral genes and environmental variables shape eukaryotic carbon cycling pathways. Specifically, we ask the following questions: (1) Do spatial patterns in viral and carbon cycling gene expression reflect heterogeneity in multi-scale environmental variables? (2) What are the mechanistic links between environmental drivers, viral gene activity, and host carbon cycling gene regulation? Our results clarify how viruses bridge environmental change and eukaryotic carbon cycling in riverine ecosystems [14,15].

2. Materials and Methods

2.1. Basin Description and Sampling

Water samples were collected in September 2019 from six stations (115.8253° E–117.7202° E, 39.1179° N–40.0371° N) along a northwest–southeast transect of the Yongding River (Figure 1a,b). Based on landscape characteristics, the sampling area was divided into three sections: the mountain section (M1, M2), the plain section (P3, P4), and the coastal section (C5, C6). At each station, water samples were collected at a depth of 0.5 m below the surface and split into two subsamples. One subsample (3 L) was used for metatranscriptomics analysis, and the other (2.5 L) was reserved for water chemistry analysis. Three replicates were collected at each station. All water samples were immediately pre-filtered through a 200 μm mesh sieve to remove large particles and then through 0.22 μm pore-size polycarbonate membranes (diameter of 142 mm, Millipore, Burlington, MA, USA). Filtration was performed within 15 min to minimize gene expression changes and RNA degradation. RNA samples were stored at −80 °C until extraction.

2.2. Determination of Environmental Variables

A multi-scale approach was used to collect environmental data at the basin, reach, and site levels. At the basin scale, a 90 m resolution Digital Elevation Model (DEM) was processed using a Geographic Information System (GIS), and the D8 flow-direction algorithm was applied to delineate the catchment area for each sampling station. Land use types were derived from the China Multi-Temporal Land Use Remote Sensing Monitoring Dataset (CNLUCC) and classified into six categories: cropland, forest, shrub, grassland, water, barren, and construction. The areal proportion of each category within the catchment was subsequently calculated. Mean annual precipitation (MAP) and mean annual temperature (MAT) were obtained from the National Earth System Science Data Center, with spatial resolution matched to site coordinates.

At the reach scale, each sampling site served as the center for delineating 0.5~5 km river segment units along the channel. Two morphological metrics, sinuosity (S) and slope gradient (G), were computed in ArcGIS (v10.8). S was defined as the ratio of actual river length to the straight-line distance between endpoints, reflecting channel curvature. G was calculated as the ratio of the elevation difference between endpoints divided by the river length, characterizing the energy gradient of the water flow. All calculations were manually verified to ensure accuracy. At the site scale, environmental data were obtained through a combination of field measurements and laboratory analyses. A Manta+ multi-parameter water quality analyzer (Eureka, CA, USA) was used to measure in situ physicochemical variables, including water temperature (WT), dissolved oxygen (DO), pH, oxidation–reduction potential (ORP), electrical conductivity (EC), and turbidity.

Laboratory analyses were performed to determine total nitrogen (TN), ammonia nitrogen (NH_4_^+^-N), total phosphorus (TP), biochemical oxygen demand (BOD_5_), permanganate index (COD_Mn_), chemical oxygen demand (COD_Cr_), and heavy metal concentrations (As, Se, Pb, Cu, Zn, Cr, and Cd) following established methods [10]. The final environmental database contained 28 indicators across three scales, effectively capturing the spatial heterogeneity along the “mountain–plain–coastal” gradient.

2.3. RNA Extraction, cDNA Synthesis, and Metatranscriptomic Sequencing

For RNA extraction, samples were supplemented with a 1:1 mixture of 0.5 mm and 0.1 mm glass beads and subjected to vortex oscillation to ensure thorough cell disruption. Total RNA was extracted using the RNeasy Mini Kit (Qiagen, Hilden, Germany). RNA quality was verified using gel electrophoresis, A260/A280 ratios measured on a Nanodrop One spectrophotometer (Thermo Fisher Scientific, Walthman, MA, USA), Qubit 2.0 fluorometric (Invitrogen, MA, USA) measurement, and Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) measurement of RNA integrity numbers (RINs) [8]. Eukaryotic mRNA was enriched by the oligo(dT) magnetic bead method [16], which selectively binds poly(A)-tailed transcripts.

The enriched mRNA was fragmented using fragmentation buffer, and first-strand cDNA was synthesized with random hexamer primers. Second-strand cDNA was subsequently generated using DNA polymerase I and RNase H. The double-stranded cDNA was purified with AMPure XP magnetic beads. Subsequent library construction involved end repair, A-tailing, and adaptor ligation, followed by a size selection of ~400 bp fragments. Library quality was assessed at three stages: (1) preliminary quantification with a Qubit 2.0 fluorometer, (2) the evaluation of insert size using an Agilent 2100 Bioanalyzer, and (3) the accurate determination of effective concentration (>2 nM) via qPCR. Qualified libraries were sequenced on the Illumina HiSeq platform (Illumina, San Diego, CA, USA) using the PE150 mode, yielding approximately 10 Gb of clean data per sample on average. All procedures were performed in a laminar flow clean bench, and consumables were treated with DEPC water to prevent RNase contamination.

2.4. Bioinformatic Analysis

Raw metatranscriptomic sequencing reads were processed using the following workflow to characterize both viral and eukaryotic host activities, with a focus on carbon cycling. First, adapter sequences and low-quality bases were trimmed. Ribosomal RNA (rRNA) was then depleted bioinformatically using SortMeRNA (version 2.1b) (https://github.com/biocore/sortmerna/, accessed on 1 January 2025). The remaining high-quality mRNA reads were assembled de novo into transcripts using idba_tran (v1.1.3) (https://github.com/loneknightpy/idba, accessed on 1 January 2025). These putative viral contigs were clustered to create a non-redundant viral gene catalog, as detailed in the Supplementary Materials (Table S1).

For eukaryotic host analysis, it is crucial to note that eukaryotic mRNA was experimentally enriched during library preparation using Oligo(dT) bead capture, which selectively targets polyadenylated transcripts characteristic of eukaryotes. Protein-coding regions within the eukaryotic-enriched transcriptome were predicted using Prodigal (https://github.com/hyattpd/Prodigal.git, accessed on 1 January 2025). A comprehensive non-redundant gene catalog was constructed and functionally annotated against the eggnog (http://eggnog6.embl.de/, accessed on 1 January 2025), KEGG, and GO databases [17]. To ensure the specificity of our eukaryotic host analysis, genes annotated with taxonomic assignments corresponding to bacteria or archaea were systematically filtered out from downstream ecological interpretations.

To ensure comprehensive profiling of carbon cycling pathways, our annotation strategy integrated multiple databases. General metabolic mapping was performed using KEGG (https://www.kegg.jp/, accessed on 1 January 2025) and eggNOG (https://github.com/eggnogdb/eggnog-mapper, accessed on 1 January 2025). Specifically, for organic carbon degradation pathways, we supplemented our analysis with annotation against the carbohydrate-active enzymes (CAZy) database. Key carbon fixation and conversion pathways (e.g., reductive citrate cycle, Calvin cycle) were manually reviewed by tracing enzyme commissions (ECs) and cross-referencing with the MetaCyc pathway database. Gene expression abundance for all analyzed genes (viral and host) was finally quantified and normalized using StringTie (https://www.gencodegenes.org/, accessed on 1 January 2025). Following assembly and gene prediction, open reading frames (ORFs) were taxonomically and functionally annotated using the eggNOG database. ORFs exclusively assigned to bacterial or archaeal taxa were excluded from subsequent analyses targeting eukaryotic and viral functional genes.

2.5. Statistical Analysis

The analysis of group differences in environmental parameters was performed using paired-sample t-tests. Statistical significance was defined as meeting both criteria: p < 0.05 and |t| > 2.776. Principal coordinates analysis (PCoA) was applied to environmental variables. Viral functional genes and carbon cycling-related functional genes were compared, and PCoA was subsequently conducted based on heatmap profiles of these functional genes. Spearman correlation analysis was used to examine relationships among environmental factors, viral gene expression, and functional gene expression, with statistical significance assessed using Student’s t-tests in SPSS 20. Co-occurrence network analysis was used to explore the synergistic regulation of host carbon cycling by environmental factors and viruses. Random forest analysis was employed to evaluate the influence of environmental and viral variables on carbon cycling gene expression. In addition, structural equation modeling (SEM) was applied to investigate the pathways linking carbon cycling processes to environmental drivers. A phylogenetic analysis of specific gene families was performed to elucidate their evolutionary origins and relationships with homologous genes.

3. Results

3.1. Spatial Heterogeneity of Environmental Variables Across Multi-Scales

PCoA revealed distinct environmental characteristics among the different segments of the Yongding River Basin, specifically, the mountain, urban plain, and coastal sections (Figure 1c). Annual precipitation was higher in the mountain (592.3608.0 mm) and coastal (570.1649.8 mm) sections compared to the plain section (548.0~554.3 mm). Mean annual temperature increased steadily from upstream to downstream. Land use patterns varied significantly: the plain section had the highest portion of cropland (56.2% to 78.9%), whereas the mountain region was dominated by forest cover (82.2% to 85.5%) and contained more shrub and grassland. Construction land was lowest in the mountain section, comprising only 0.41% to 3.28%. In contrast, the coastal section had the highest proportions of water and barren land.

Slope gradient decreased from the mountain section toward the downstream region, with the steepest slopes occurring in the mountain region. Sinuosity was high in both plain and coastal sections and increased sharply in the mountain section at the 2~5 km scale. Water quality exhibited distinct regional variations: ORP was the highest in the plain section, while EC peaked in the coastal section; turbidity, DO, and NH_4_^+^-N all increased downstream. TN, TP, and heavy metals (especially As, Pb, Zn, and Cr) were higher in the coastal and plain sections.

3.2. Spatial Heterogeneity in the Transcription of Eukaryotic-Related Viral Functional Genes

Sequencing of the 18 libraries yielded approximately 600 million raw paired-end reads (~360 Gb of data, with Q30 > 90.19%). After stringent quality control, adapter trimming, and removal of rRNA reads, an average of 2 million clean reads per sample were obtained for assembly. The total number of functionally annotated genes analyzed throughout the study is 41,936 (FPKM > 10). Our metatranscriptomic data captured viral transcripts infecting eukaryotes. The viral functional genes identified included uncoupled ATPase activity (77 sequences), geminivirus coat protein/nuclear export factor BR1 (8), helical virus capsid (47), major capsid proteins (MCPs, 54), RNA helicase (62), RNA processing (93), and viral S domain (60). Heatmap analysis indicated that viral gene abundance followed a unimodal (increase-then-decrease) pattern along the mountain–plain–coastal gradient (Figure 2a). PCoA revealed significant spatial clustering: coastal sites separated from the mountain and plain regions in their viral gene profiles (Figure 2b).

3.3. Spatial Heterogeneity in the Transcription of Eukaryotic Carbon Cycling Functional Genes

Based on eggNOG, KEGG, and CAZyme annotations (Table S2), the identified eukaryotic carbon cycling genes were grouped into carbon fixation (4 categories), conversion (13 categories), and metabolism (12 categories). The genes for fixation (e.g., rbcL, PEPC) showed the lowest abundance, but all three functional groups had elevated transcriptional levels in the plain section (especially glycine decarboxylase; 27 sequences), implying enhanced photorespiration.

Carbon conversion genes (carbonic anhydrase, UDP-glucuronate decarboxylase) and metabolism genes (malate dehydrogenase, acetyl-CoA carboxylase) peaked in the plain section (Figure 3a), supporting active carbon skeleton modification, glycolysis, and TCA cycle metabolism.

Abundance of all three gene groups followed an increase-then-decrease pattern (mountain/plain/coastal gradient; Figure 3a). PCoA and ANOSIM confirmed spatial differences across the sections (Figure 3b): the mountain region was dominated by terrestrial organic matter degradation, the plain by energy metabolism, and the coast by environmental adaptation.

3.4. Linking Multi-Scale Environmental Gradients to Viral Gene Expression Landscape

Spearman correlation analysis revealed the relationships between viral gene expression and multi-scale environmental factors (Figure 4). Uncoupled ATPase activity genes positively correlated strongly with cropland, water temperature (WT), ORP, and BOD5 (p < 0.01). The geminivirus coat protein was linked to 500 m sinuosity (p < 0.01). The helical virus protein gene was positively correlated with construction area and the concentrations of heavy metals (Pb, Zn, Cr) (p < 0.01), but negatively to pH (p < 0.01). MCP genes expressed by large eukaryotic DNA viruses were highly positively correlated with cropland area, WT, and heavy metals (Pb, Zn) (p < 0.01), but negatively to MAP (p < 0.01). The RNA helicase gene correlated with forest area, shrub area, 3 km slope, and 500 m sinuosity (p < 0.01), but negatively with cropland area and turbidity (p < 0.01). RNA processing-related genes were positively correlated with NH_4_^+^-N, ORP, BOD5, heavy metals (Pb, Zn, Cr), and construction area, but negatively to pH (p < 0.01). The expression of the viral capsid protein (S domain) gene was positively correlated with cropland area and WT (p < 0.01), but negatively with EC, TP, Se, and Cu (p < 0.01).

Of the 176 correlation tests performed between environmental factors and virus gene expression, 35 (19.89%) were significant (p < 0.01), 26 were positive, and 9 were negative. Cropland area showed the most frequent positive correlations (with uncoupled ATPase, MCP, and S domain genes) and a negative correlation with RNA helicase.

3.5. Co-Occurrence Network Linking Environmental Factors to Gene Transcriptional Landscapes

MCP genes from large eukaryotic DNA viruses correlated positively with carbon metabolism (TCA) and conversion (decarboxylases) genes (p < 0.05), especially at cropland-dominated sites (Figure 5 and Figure S1). Other viral carbon cycling links included geminivirus coat protein/BR1 family and carbonate metabolism (p < 0.05). RNA helicase expression was negatively correlated with carbon fixation genes in natural vegetation (e.g., forest, shrub, grassland), but positively correlated with sugar transporter genes. Viral activity appears to facilitate DOC release by inhibiting host photosynthesis (Figure 4 and Figure 5). Site-scale environmental factors also positively affected RNA processing genes and genes involved in carbon conversion and metabolism (Figure 5).

3.6. Response Pathways of Carbon Cycle Gene Expression to Environmental and Viral Effects

Random forest analysis identified significant associations among environmental factors, viral gene expression, and eukaryotic carbon cycling gene transcription (Figure 6). For carbon fixation genes, MAP, cropland area, 5 km river sinuosity, TN concentration, and MCP of large eukaryotic DNA viruses showed strong influences (importance score > 6, p < 0.01). For carbon conversion genes, the aforementioned environmental factors similarly demonstrated high relevance (importance score > 6). In the case of carbon metabolism genes, a broader suite of factors contributed to the model; nevertheless, cropland area, 5 km sinuosity, MAP, TN, and MCP consistently emerged as the primary drivers of spatial variation in transcriptional activity.

Based on the random forest results, structural equation modeling (SEM) indicated no significant multicollinearity among the selected environmental predictors (i.e., MAP, cropland area, 5 km river sinuosity, and TN concentration). We developed piecewise structural equation models (SEMs) for the MCP of large eukaryotic DNA viruses, uncoupled ATPase, and carbon cycling genes (Figure 7). These models retained statistically significant and biologically meaningful pathways while maintaining explanatory power through parsimonious construction. Our analysis revealed distinct pathways through which different factors influenced carbon cycling genes and virus functional genes. None of the factors showed a significant relationship (coefficients < 0.02).

For the major capsid protein (MCP) of large eukaryotic DNA viruses, cropland area exerted the strongest influence, with cropland expansion showing a positive effect on MCP (coefficient = 0.96; Figure 7a), which in turn indirectly enhanced multiple aspects of carbon cycling genes (Figure 7b; indirect effects: metabolism = 0.53, conversion = 0.62, fixation = 0.27). In contrast, the remaining factors showed opposing influences: TN exerted a negative effect on MCP (−0.36), while 5 km river sinuosity and MAP had relatively minor negative effects on MCP (coefficients < 0.1).

For uncoupled ATPase activity, neither TN nor cropland area displayed direct significant pathways (coefficients < 0.15). However, 5 km river sinuosity significantly suppressed the increase in its functional gene (coefficient = −0.80). MAP exhibited a negative total effect on the functional aspect of uncoupled ATPase (−0.73). Nevertheless, because uncoupled ATPase served as a control group for MCP, it had only a limited direct influence on carbon cycling genes. MCP demonstrated stronger effects on carbon fixation (0.28), conversion (0.64), and metabolism (0.56) compared to uncoupled ATPase, whose corresponding values were 0.02 (fixation), 0.08 (conversion), and −0.06 (metabolism), respectively (Figure 7a).

In summary, the multi-diversity structural equation model (Fisher’s C = 46.82, p = 0.54) incorporated key pathways from both viral functional gene groups and exhibited superior explanatory power relative to single-taxon models (R^2^ metabolism = 0.56, R^2^ conversion = 0.75, R^2^ fixation = 0.78; Figure 7a). Cropland expansion significantly enhanced the expression of carbon cycling genes (through both direct and indirect positive regulation) and maintained a dominant role in carbon metabolism, conversion, and fixation, thereby promoting a positive carbon cycle. In contrast, 5 km river sinuosity exhibited negative effects on carbon cycling genes, reflecting the opposing influence of reach-scale factors on carbon cycling genetic potential.

3.7. Phylogenetic and Spatial Distribution Analysis of Viral Genes

As shown in Figure 8, uncoupled ATPase activity genes were distributed in viruses infecting diverse hosts, especially crustaceans and aquatic invertebrates. Phylogenetic analysis revealed that the uncoupled ATPase activity genes most similar to homologs from viruses infecting these hosts (Figure 8a, Table S3) displayed high expression in the mountain section. MCP genes predominantly clustered with large DNA viruses (Iridoviridae and Ascoviridae) and giant viruses. They infect a broad range of hosts, including planktonic eukaryotic algae, protozoa, invertebrates, fish, and amphibians (Figure 8b and Table S3). The expression of these MCP genes peaked in the plains.

4. Discussion

To our knowledge, no previous study to date has systematically dissected how eukaryotic viral gene expression mechanistically couples environmental gradients to carbon cycling in riverine ecosystems using omics approaches. Our metatranscriptomic analysis provides the first in situ demonstration of how environmental drivers and viral genes interact in rivers to orchestrate carbon cycling gene expression. We show that multi-scale environmental heterogeneity (notably, land use and channel morphology) substantially modulates viral activity, which in turn acts as a biological amplifier reprogramming host carbon metabolic pathways.

4.1. Sectional Differences in Viral and Carbon Cycling Genes Along Environmental Gradients

We found clear spatial heterogeneity in both viral functional gene and eukaryotic carbon cycling gene expression across the “mountain–plain–coastal” continuum of the Yongding River Basin (Figure 1c, Figure 2 and Figure 3). Viral genes, especially uncoupled ATPase and major capsid proteins (MCPs) of large DNA viruses, were sharply upregulated in plain sections (P3, P4), showing a strong positive association with TCA cycle and decarboxylase gene expression (Figure 3). This suggests viral amplification in plain regions, boosting lytic infection and host metabolic reprogramming. By contrast, the mountain section (M1, M2) showed a dominance of helical virus genes and RNA processing genes, patterns that may facilitate the release of terrestrially derived phenolic DOC by lysing fungal hosts. Coastal zones (C5, C6) displayed peak RNA helicase gene expression, which likely suppresses photosynthetic carbon fixation by inhibiting the expression of related genes (e.g., rbcL).

These patterns reveal that viruses employ distinct ecological strategies, ranging from lytic infection in agricultural plains to persistent infection or metabolic suppression in coastal zones, thus tailoring their impact on carbon cycling to specific environmental contexts. This echoes findings from Baltic Sea virome studies reporting ecological functional turnover between viral types along salinity gradients [18] and the concept of aquatic RNA viruses as regulators of carbon flow [19].

Sites dominated by cropland exhibited the highest MCP gene abundance (e.g., Phycodnaviridae) (Figure 2), likely linked to phytoplankton blooms that favor viral replication. These MCPs correlated positively with host TCA and decarboxylase gene expression (Figure 3), suggesting the role viruses play in accelerating the conversion of POC to DOC and enhancing carbon turnover instead of sedimentation. At the basin scale, cropland robustly drives MCP transcription via promotion of algal hosts, mirroring the urbanization effects on microbial function observed in a previous study [20]. The lower abundance of RNA viruses (e.g., Marnaviridae) in vegetation-dominated sections (forest, shrub, grassland) (Figure 2a) signifies their adaptation to coastal and estuarine environments. Negative correlations between helicase and rbcL expression reinforce the idea that the viral suppression of photosynthetic hosts augments organic carbon release. The trend of elevated carbon conversion genes (e.g., carbonic anhydrase) in the plain section (Figure 3a) reflects a shift in energetic allocation toward decomposition rather than fixation, paralleling the demands placed on fungi and protists as they break down complex substrates [21].

4.2. Environmental Gradients Drive Viral Amplification and Modulate Carbon Fate

Our multi-scale analysis demonstrates that environmental factors acting across scales—from basin (e.g., cropland) and reach (e.g., sinuosity) to site (e.g., TN)—exert significant, predominantly positive effects on viral gene expression, which in turn regulates carbon cycling pathways (Figure 4 and Figure 5). The elevated expression of genes encoding the major capsid protein (MCP) of large eukaryotic DNA viruses (e.g., Phycodnaviridae), critical for virion stability and indicative of broad host ranges (algae, protists), signifies intense lytic activity. This supports the viral shunt theory: lysis diverts organic carbon from particulate sedimentation into the dissolved pool, potentially enhancing downstream respiration and reducing local carbon storage [15,22]. The co-occurrence of MCP genes with ATPase genes on the plains, along with RNA helicases and coastal RNA processing genes, points to a suite of viral strategies—from classic lytic to metabolic hijacking or even chronic suppression—that collectively structure riverine carbon fate [12,13]. This process aligns with riverine systems’ status as CO_2_ sources [1,23] and contrasts with the carbon-sequestering “viral pump” found in oceans [19,20].

Beyond lysis, viral-driven metabolic reprogramming, e.g., ATPase genes redirecting host energy, further shifts river carbon cycling toward enhanced DOC production and CO_2_ release in agricultural regions. This is supported by the strong upregulation of glycolytic and TCA pathways and the repression of rbcL and carbon fixation genes (Figure 3 and Figure 5). The metabolic flexibility observed here mirrors recent findings under nutrient-enriched conditions [3]. Such metabolic flexibility, while potentially conferring adaptive advantages to organisms in fluctuating environments, comes at the cost of net carbon loss from the system [23,24,25].

4.3. Mechanistic Pathways: Environmental Drivers Reshaping Carbon Cycling via Viruses

Random forest and SEM analyses (Figure 6 and Figure 7) reveal a mechanistic cascade: multi-scale landscape features (cropland area, sinuosity, TN) transmit signals through the upregulation or repression of viral gene expression, which, in turn, drives host carbon cycle gene transcription. At the basin scale, cropland strongly enhanced the expression of the major capsid protein (MCP) of large eukaryotic DNA viruses (SEM path coefficient = 0.96, Figure 7a), which subsequently induced a multi-faceted upregulation of carbon cycling genes (indirect effects: metabolism = 0.53, conversion = 0.62, fixation = 0.27), resulting in a 1.8-fold rise in TCA transcription in the plain section. Notably, as a capsid protein of lytic viruses, MCP promotes carbon metabolism more strongly than fixation, likely due to the rapid lysis of algal hosts [26,27,28]. However, the nutrient inputs associated with cropland expansion stimulate algal growth, leading to a direct effect where carbon fixation (0.34) exceeds carbon metabolism (0.06) [29,30,31].

In contrast, increased reach-scale river sinuosity (5 km sinuosity) exerted a pronounced negative effect on viral functional genes and consequently suppressed carbon cycling potential, likely mediated by prolonged water residence time that alters viral–host encounter rates, especially uncoupled ATPase (SEM path coefficient = −0.80) [32,33]. At the site scale, elevated TN coincided with suppressed expression of both MCP and key carbon fixation genes, underscoring the nutrient-mediated control of viral–host interactions [34,35].

Anthropogenic land use (notably cropland) amplifies viral abundance and replication, while natural vegetation and water quality filter these effects—together reshaping virus–host interaction networks (Figure 7) [36]. These findings position viruses as rapidly responsive, agentic biological nodes transmitting environmental information to eukaryotic carbon cycling machinery, a concept gaining traction in both marine and soil microbiome studies [9,37,38]. Our multi-diversity SEM (Fisher’s C = 46.82, p = 0.54) further integrates these pathways with improved explanatory power (R^2^ > 0.5), highlighting the utility of a functional-gene framework for deciphering how terrestrial and aquatic landscapes jointly regulate carbon cycling dynamics through viral intermediaries.

4.4. Host Specificity and Ecological Implications of Viral Gene Distributions

Phylogenetic mapping (Figure 8) affirms that the spatial variation in viral functional gene abundance is tightly linked to the biogeography of their eukaryotic hosts, including algal, fungal, protozoan, and invertebrate communities. Uncoupled ATPase activity gene enrichment in mountain sections (Figure 8a) echoes the dominance of crustacean and invertebrate hosts in benthic communities, fitting a regime dominated by the degradation of allochthonous organic matter (e.g., leaf litter and terrestrial detritus) and sediment resuspension. Reduced ATPase gene prevalence in forested regions may reflect reduced runoff and host scarcity. MCP genes of large DNA viruses peaking in plains (Figure 8b) correspond to phytoplankton blooms, with lysis shifting carbon from POC to DOC and ultimately stimulating CO_2_ release (see TCA cycle and decarboxylase upregulation; Figure 3). These results demonstrate viral control over primary productivity and CO_2_-concentrating mechanisms [39], processes that are tightly linked to land use.

In estuarine and coastal transition zones, MCP and ATPase viral activities declined (Figure 8), likely due to host constraints (freshwater taxa in saline water; [40]), while persistent viral infection strategies (RNA helicase activity) take precedence, suppressing carbon fixation and reshaping local biogeochemical flows.

Overall, viral functional gene distribution and, by extension, their “amplification” of carbon cycle signals, is dictated not by environmental factors alone but by multi-scale landscape-driven host filtering (Figure 8). This extends marine paradigms of virus–host specificity [41] to the riverine and estuarine domains.

5. Conclusions

Our results highlight viruses as keystone intermediates, directly transmitting and amplifying multi-scale environmental information to reshape eukaryotic riverine carbon cycling. Our main findings include the following: (1) Viral genes, particularly those encoding major capsid proteins (MCPs) and ATPases, are highly responsive to multi-scale environmental variations, serving as potential bioindicators of anthropogenic and natural influences. (2) Virus–host interactions can enhance carbon conversion and metabolism in plains while suppressing carbon fixation in coastal zones, reflecting contrasting viral strategies under different environmental contexts. (3) Multi-scale environmental factors, especially croplands (basin scale), 5 km channel sinuosity (reach scale), and total nitrogen (site scale), modulated viral MCP gene expression, which in turn regulated carbon cycling pathways, potentially reducing the carbon sequestration capacity of river ecosystems. Forecasting the fate of carbon in global river networks demands the integration of viral ecological processes into biogeochemical modeling frameworks, especially in the face of rapidly intensifying anthropogenic land use and climate-driven hydrological change.