Environmental RNA-Based Metatranscriptomics as a Novel Biomonitoring Tool: A Case Study of Glyphosate-Based Herbicide Effects on Freshwater Eukaryotic Communities

Abstract

Traditional morphology- and molecular-based biodiversity surveys provide essential information on species composition and diversity, but they rarely provide information about the physiological states of organisms, which are key indicators of ecosystem health. Environmental RNA (eRNA) has the potential to significantly enhance biomonitoring by providing insights beyond species detection. Recent studies suggest that extra-organismal RNA released into the environment could help identify differentially expressed genes of single species. However, the feasibility of eRNA-based metatranscriptomics on complex environmental samples, containing both extra-organismal and organismal eukaryotic RNA, remains untested due to numerous experimental and analytical challenges. In this study, we explored the potential of eRNA-based metatranscriptomics, enriched for eukaryotes, as a tool to monitor environmental stress. We used outdoor mesocosms to examine the acute effects of a glyphosate-based herbicide (GBH) on gene transcription across diverse freshwater eukaryotic taxa. Our metatranscriptomics data revealed diverse eukaryotic taxa spanning multiple trophic levels, including phytoplankton, zooplankton, ciliates, and aquatic insects. GBH treatment significantly altered the relative transcript abundances of most eukaryotic classes, with longer-lived taxa demonstrating greater tolerance compared to shorter-lived taxa. Differential expression analysis showed more gene downregulation than upregulation in response to GBH, likely due to its acute toxicity. Many differentially expressed genes were involved in molecular pathways associated with responses to GBH exposure, such as oxidative stress response and detoxification. Our results demonstrate that eRNA-based metatranscriptomics captures transcriptional signals from diverse aquatic eukaryotic taxa, providing insights into functional gene expression. As such, its application to support environmental monitoring of aquatic ecosystems warrants further exploration.

Keywords: biodiversity survey; environmental biomonitoring; environmental transcriptomics; gene expression; oxidative stress response; plankton.

FIGURE 1

Flowchart illustrating the general bioinformatic and post‐bioinformatic analysis workflow used in this study.

FIGURE 2

Effect of GBH treatment on relative transcript contributions of freshwater eukaryotic classes. (A) Relative transcript abundances of classes in each sample. “Others” represents the sum of eukaryotic classes with < 1% relative transcript abundance in every sample. “Unassigned_Eukaryota” represents the reads assigned to Eukaryota but not to the class level. (B) Comparison of relative transcript abundances between treatment and control for nine major classes. Asterisks (‘*’) indicate significant differences revealed by either LMM or rLMM, while “NS” denotes non‐significant difference.

FIGURE 3

Comparison of samples based on transcript‐derived taxonomic composition and overall gene expression. (A) Principal coordinate analysis (PCoA) of Bray–Curtis dissimilarity based on transcript read counts assigned to eukaryotic classes (rarified to 600,000 reads per sample). (B) Heatmap of sample‐to‐sample distances based on gene expression data. Distances were calculated using eukaryotic contigs with at least 10 reads in at least four samples, and hierarchical clustering was performed on the sample distances.

FIGURE 4

Heatmap showing the number of differentially expressed contigs assigned to KEGG Level 2 pathways across five taxonomic categories. Each row represents one of the KEGG Level 2 pathways, and each column represents the number of contigs that were either downregulated (left five columns) or upregulated (right five columns) in response to the GBH treatment.

FIGURE 5

Heatmaps showing expression patterns of differentially expressed contigs associated with oxidative stress response and xenobiotic metabolism pathways across nine major taxonomic classes. Z‐scores were calculated for each contig for clustering. (A) Downregulated contigs associated with three KEGG pathways related to oxidative stress response. (B) Upregulated contigs associated with three KEGG pathways related to oxidative stress response. (C) Downregulated contigs associated with the “Metabolism of xenobiotics by cytochrome P450” pathway. (D) Upregulated contigs associated with the “Metabolism of xenobiotics by cytochrome P450” pathway.