Integrated miRNA-mRNA transcriptomic analysis of hepatopancreas reveals molecular mechanisms in Macrobrachium rosenbergii under graded nitrite stress

TL;DR

This study explores how Macrobrachium rosenbergii prawns respond to nitrite stress at the molecular level, revealing key genes and pathways involved in their survival strategies.

Contribution

The first integrated mRNA-miRNA regulatory network analysis in Macrobrachium rosenbergii under nitrite stress.

Findings

Differential gene expression in hepatopancreas under low and high nitrite concentrations highlights distinct metabolic and signaling pathways.

miRNA profiling identified concentration-dependent expression patterns, including novel miRNAs and stress-responsive miRNAs.

qPCR validation confirmed key stress-responsive genes and miRNAs, supporting the sequencing data.

Abstract

Nitrite accumulation poses a significant threat to aquatic organisms in intensive aquaculture systems. Macrobrachium rosenbergii, a commercially vital freshwater prawn, exhibits adaptive responses to environmental stressors, yet the molecular mechanisms underlying nitrite tolerance remain poorly understood. This study employed integrated mRNA and miRNA transcriptomics to dissect the regulatory networks activated in M. rosenbergii hepatopancreas under acute nitrite stress (0, 40, and 87.25 mg/L nitrite-N over 48 h). High-throughput sequencing revealed 640 and 912 differentially expressed genes (DEGs) in low-concentration (LC) and high-concentration (HC) groups, respectively, compared to controls (CK). In the LC group, enrichment was predominantly observed in ribosome biogenesis (96 genes, p < 0.001). Conversely, the HC group was characterized by the significant modulation of the PPAR signaling (11 genes), glycerophospholipid metabolism, and the citrate cycle (p < 0.005). Calcium signaling and MAPK pathways may be central to stress adaptation across both groups. miRNA profiling revealed 17 downregulated and 2 upregulated miRNAs within the HC group relative to the CK group, with miR-193-y, miR-263-x, and miR-145-x implicated in metabolic regulation. Notably, novel miRNAs (e.g., novel-m0087-3p) showed concentration-dependent expression. qPCR validated the consistency of sequencing data, confirming stress-responsive genes (P53, HORMA, SLC25a28) and miRNAs. This study is the first to integrate the mRNA-miRNA regulatory network in Macrobrachium rosenbergii to elucidate the response mechanism to nitrite stress, emphasizing the role of metabolic reprogramming and signal pathway regulation as key survival strategies. It provides a foundation and novel perspectives for the molecular evolution of nitrite adaptability in aquatic animals and breeding programs.

Figures

Schematic diagram of experimental process.

Unigene annotation results.(A) Annotation information obtained from seven different databases, (B) Annotation information based on the Gene Ontology database, (C) The most common KEGG pathways.

Analysis of differential gene expression and KEGG enrichment of LC compared with the CK.(A) Differentially regulated genes under LC compared with CK. (B) KEGG enrichment analysis of the DEGs under LC and CK. (C) Heatmap analysis of differentially expressed genes between the LC and CK groups.

Analysis of differential gene expression and KEGG enrichment of HC compared with the CK.(A) Differentially regulated genes under HC compared with CK. (B) KEGG enrichment analysis of the DEGs under HC and CK. (C) Heatmap analysis of differentially expressed genes between the HC and CK groups.

Analysis of differential gene expression and KEGG enrichment of LC compared with the HC.(A) Differentially regulated genes under LC compared with HC. (B) KEGG enrichment analysis of the DEGs under LC and HC. (C) Heatmap analysis of differentially expressed genes between the LC and HC groups.

Effects of nitrite stress on miRNA levels of genes.(A) Differentially expressed miRNAs in M. rosenbergii (control vs high concentration), (B) Differentially expressed miRNAs in M. rosenbergii (control vs low concentration), (C) Differentially expressed miRNAs in M. rosenbergii (low concentration vs high concentration), (D) Related genes of linear change under Nitrous acid stress in M. rosenbergii. EF1-α gene was selected as an internal reference gene.

Go enrichment analysis of the miRNA.(A) low concentration compared with control condition, (B) high concentration compared with control condition, (C) low concentration compared with high concentration.

Top of 20 KEGG enrichment analysis of the miRNA.(A) low concentration compared with control condition, (B) high concentration compared with control condition, (C) low concentration compared with high concentration.

Full text

Introduction

Macrobrachium rosenbergii, as a freshwater shrimp species with strong adaptability and high economic value, has enormous commercial value worldwide and is an important species in freshwater aquaculture. Generally, it grows to a length of 150–200 mm within one year, with the maximum body length of male individuals reaching up to 400 mm [1]. This species is characterized by its ease of cultivation, substantial size, and high consumer preference. Due to its strong adaptability and fast growth rate, M. rosenbergii has become one of the important freshwater cultured shrimp in Southeast Asia.

The commercial value of M. rosenbergii has driven a rapid increase in market demand and industrial growth. However, high-density farming and improper feeding practices have caused significant elevations in nitrite and ammonia nitrogen levels. Consequently, these environmental stressors severely impair the health of the prawns [2]. At the same time, in response to the policy of water conservation and environmental protection, part of the aquaculture industry will use recycled water, which may lead to a sharp increase in the level of ammonia accumulation, and ammonia nitrogen can be rapidly oxidized to nitrite, thus posing a huge threat to the survival of aquatic animals. Some studies have shown that excessive nitrite in water will reduce the pH of hemolymph after it enters the body of fish and shrimp, thus affecting the transport of oxygen in the body, disrupting nitrogen excretion and damaging organs, causing changes in physiological and biochemical factors and tissue structure of shrimp and crabs [3–5]. In the stress response of Macrobrachium rosenbergii, an environmental concentration of 0.15 mg/L is relatively safe for larvae weighing 0.07 g. However, the combined presence of 2.5 mg/L total nitrogen in the water exhibits an antagonistic effect on the larvae [6]. Another study demonstrated that nitrite stress delays larval development, reduces growth rates, and increases mortality in Macrobrachium rosenbergii [7]. Nevertheless, the regulatory mechanisms underlying the response of Macrobrachium rosenbergii to nitrite stress remain unclear.

The hepatopancreas in shrimp serves as a crucial detoxification organ, similar to the liver tissue found in higher animals, and plays a significant role in energy metabolism and immunity [8]. As a result, it is often used as a model organ to explore the response of crustaceans to various environmental stresses. Recently, high-throughput sequencing technology has rapidly developed and become popular in animal research. Transcriptome technology, in particular, has been widely used to study nitrite stress in Litopenaeus vannamei [9–11]. A total of 43,405 transcripts were obtained through high-throughput sequencing for analysis of immune pathway genes and their differential expression in M. rosenbergii infected with spirochetes [12]. Comparative transcriptome analyses were conducted on the gill and hepatopancreas of M. rosenbergii under heavy metal Cd^2+^ stress, as well as on the response of larvae of M. rosenbergii to metamorphosis and salinity exposure [13,14]. Additionally, the response of juvenile M. rosenbergii shrimps infected with Vibrio harveyi was analyzed by transcriptome technique [15].

Animal miRNAs rapidly respond to diverse environmental signals and inhibit target gene expression by degrading or repressing mRNAs. Small RNA sequencing, a high-throughput method, is commonly employed to investigate stress-responsive regulatory networks involving miRNAs. Numerous studies have demonstrated miRNA involvement in regulating key biological processes in animals, including reproductive development and immunity [16–18]. For instance, miRNAs have been specifically identified and analyzed during gonad development in Macrobrachium rosenbergii [19]. Furthermore, transcriptome sequencing has been used to analyze differences in mRNA and miRNA expression related to the immune system in M. rosenbergii infected with Spirochetes [12].

Previous studies have reported changs in mRNA and miRNA expression unhder nitrite stress in plants and many model animals [20,21]. However, systematic research on the molecular response mechanisms of Macrobrachium rosenbergii—particularly the synergistic regulatory network between mRNA and miRNA—is currently lacking. Therefore, this study aims to systematically elucidate the key response pathways and their regulatory networks in Macrobrachium rosenbergii under nitrite stress through integrated high-throughput sequencing analysis of mRNA and miRNA.The objectives of this study are to identify differentially expressed mRNAs and miRNAs involved in nitrite stress responses, construct and analyze mRNA-miRNA interaction regulatory networks, identify core regulatory hubs, and simultaneously reveal key biological pathways and molecular mechanisms. The results of this study will not only deepen the understanding of the molecular mechanisms underlying Macrobrachium rosenbergii’s response to nitrite stress but also provide important theoretical references for elucidating the universal strategies of crustaceans in adapting to environmental stress.

Materials and methods

Chemicals and experimental design

Sodium nitrite (CAS: 7632-00-0, purity 99.9%) was procured from Aladdin Co., Ltd, Shanghai, China. A stock solution was prepared by dissolving the sodium nitrite powder in ultra-pure water.

Healthy M. rosenbergii were obtained from a local farm in Gaoyou (Jiangsu, China). Prior to experimentation, prawns were acclimated for 14 days in an automatic circulating water system under controlled conditions. After adaptation, A total of 270 prawns weight: 5.34 ± 1.22 g; length: 6.13 ± 0.42 cm) were randomly distributed into nine circulating water tanks (85 cm L × 65 cm W × 60 cm H; water volume: 50 L per tank), with 30 prawns per tank. Three replicate tanks were assigned to each of the following experimental groups: CK (Control): exposed to 0 mg/L sodium nitrite. LC (Low Concentration Stress): exposed to 40 mg/L sodium nitrite. HC (High Concentration Stress): exposed to 87.25 mg/L sodium nitrite. The prawns were subjected to acute nitrite exposure for a duration of 48 hours. Prawns were not fed during this period. To maintain water quality, 50% of the water in each tank was replaced daily with pre-prepared water of the corresponding nitrite concentration. Throughout the experiment, water conditions were maintained at 24 ± 0.8°C, dissolved oxygen > 6 mg/L, and a photoperiod of 12 hours light/12 hours dark (L12/D12). At the end of the 48-hour exposure, prawns were anesthetized by immersion in 10 mg/L MS-222. Hepatopancreas tissues were rapidly dissected, rinsed with ice-cold (4°C) 0.8 g/L normal saline solution, immediately snap-frozen in liquid nitrogen, and stored at −80°C for subsequent RNA extraction.

To facilitate a clearer understanding of this study, an experimental methodology flowchart has been created, which specifically illustrates the grouping parameters of the Macrobrachium rosenbergii experiment, as well as the exposure concentration and duration of sodium nitrite, sampling and processing procedures, etc (Fig 1).

Schematic diagram of experimental process.

Transcriptome sequencing and quality control

Total RNA was isolated from hepatopancreas samples using the Trizol reagent kit (Invitrogen, USA) according to the manufacturer’s protocol. RNA integrity and purity were rigorously assessed using agarose gel electrophoresis and the Agilent 2100 Bioanalyzer (Agilent Technologies, USA), ensuring the absence of RNase contamination and degradation. Next, eukaryotic mRNA was enriched from total RNA using oligo (dT) magnetic beads. Ribosomal RNA (rRNA) depletion was performed using the Ribo-Zero™ Magnetic Kit (Epicentre, Madison, WI, USA) for any potential prokaryotic contamination. The purified mRNA was fragmented using a fragmentation buffer. First-strand cDNA synthesis was performed using random hexamer primers, followed by second-strand cDNA synthesis using DNA Polymerase I, RNase H, and dNTPs. The resulting double-stranded cDNA was purified using the QiaQuick PCR extraction kit (Qiagen, The Netherlands). Library preparation involved end repair, poly(A) tailing, and ligation of Illumina sequencing adapters. The final libraries were amplified by PCR. And then,sequencing was performed on the Illumina HiSeq™ 4000 platform (Gene Denovo Biotechnology, China) to generate paired-end reads. Downstream analyses utilized high-quality clean data obtained after stringent quality filtering.

Transcriptome data analysis

Clean reads from all nine libraries (three replicates per group: CK, LC, HC) were de novo assembled into unigenes using Trinity (v2.8.5) with default parameters. Contigs were clustered using TGICL (v2.1) to form unigenes. Assembled unigenes were comprehensively annotated against multiple public databases using BLASTx (E-value threshold ≤ 1e-5). Databases included: NCBI non-redundant protein sequences (Nr), Swiss-Prot, Protein families (Pfam), Clusters of Orthologous Groups (COG/KOG), Gene Ontology (GO), and Kyoto Encyclopedia of Genes and Genomes (KEGG). The best alignment results were used for functional annotation. Annotations from different databases were integrated for each unigene. Raw read counts for each gene in each sample were generated by mapping clean reads back to the assembled transcriptome using Bowtie2 (v2.3.4.3) and RSEM (v1.3.1). Differentially Expressed Genes (DEGs) between pairwise groups (CK vs LC, CK vs HC, LC vs HC) were identified using both DESeq2 (v1.26.0) and EdgeR (v3.28.1) software packages. Genes exhibiting an absolute fold change ≥ 2.0 and a false discovery rate (FDR) < 0.05 were considered significantly differentially expressed. The results from both packages were compared, and genes identified as significant by at least one method were retained as DEGs. Putative transcription factors (TFs) were predicted by aligning protein-coding sequences against the Animal TFdb database using BLASTp (E-value ≤ 1e-5). Predicted TFs were further categorized into families based on their conserved domains.

Small RNA sequencing and quality control

For each of the 9 samples, total RNA was size-fractionated to isolate small RNA molecules (18–30 nucleotides) using 15% denaturing PAGE. Small RNA bands were excised from the gel and eluted. 3’ and 5’ adapters (Illumina TruSeq Small RNA Kit adapters) were sequentially ligated to the isolated small RNAs using T4 RNA Ligase 1 and 2 (Truncated), respectively, following the manufacturer’s protocol. The ligation products were reverse transcribed using SuperScript II Reverse Transcriptase (Invitrogen) and amplified by 15 cycles of PCR using Illumina index primers to generate cDNA libraries (140–160 bp products). Then, Small RNA libraries were sequenced single-end on the Illumina HiSeq™ 2500 platform (Gene Denovo Biotechnology, China) to generate 50 bp reads.

Small RNA sequencing data analysis

To ensure the accuracy of subsequent assembly and analysis, it is crucial to eliminate erroneous reads in the obtained sequencing engine results. These reads often contain poor quality adapters or bases, which can have a detrimental effects. Therefore, to ensure the attainment of pristine labels, the raw data underwent further filtration based on the following criteria: 1) Remove low-quality reads containing multiple bases (with a q value ≤ 20) or unknown nucleotides (denoted as N). 2) Elimination of reads lacking 3’ adaptors. 3) Elimination of reads with 5’ adaptors. 4) Delete reads that contain 3’ and 5’ adapters without a short RNA fragment in between. 5) Remove reads that contain poly A in small RNA fragments. 6) Remove reads smaller than 18nt (excluding the adaptor). The clean tags were compared with the small RNAs using the Rfam database (version 11.0) simultaneously. This comparison allowed for the identification and subsequent removal of rRNA, snRNA, sonRNA, tRNA, and scRNA. After all the tags were cleared, they were compared with the miRbase database to place existing miRNAs. In addition to existing miRNAs, the total number of miRNAs also includes known and novel miRNAs. The miRNA expression levels were determined and then normalized to TPM (transcripts per million). TargetScan (version 7.0) was employed to calculate potential targets for the miRNAs. To calculate the TPM, the actual miRNA counts were divided by the total counts of clean tags and then multiplied by 10^6^.

Validation of sequencing results by qPCR

The extraction of total RNA was conducted by the High Purity RNA Rapid Extraction Reagent (BioTeke, Beijing, China) with the manufacturer’s instructions. The cDNA was obtained by using miRNA first standard cDNA synthesis kit (Vazyme Biotech Co., Ltd.). The miRNA primers were designed based on reference sequences using Vazyme miRNA Designer software, amplified by the qPCR, and calculated by 2^-∆∆CT^ (Table 1). qPCR reactions were performed in a 20 µl reaction volume using ChamQ Universal SYBR qPCR Master Mix (Vazyme) on a QuantStudio 5 Real-Time PCR System (Applied Biosystems) in triplicate. The amplification procedure was: 95°C for 3 min and 40 cycles of 95°C for 10 s and 60°C for 30 s. Four distinct enriched genes were subjected to qPCR analysis, including P53, SLC25a28, RIO1 and HORMA. Additionally, EF1-α was utilized as an internal reference gene. Each amplification reaction was performed in triplicate. Data were analyzed using SPSS software (version 22.0) and presented as means ± standard errors of the means (SEM). Statistical significance (p < 0.05) between groups was determined by one-way analysis of variance (ANOVA) followed by Tukey’s post-hoc test for multiple comparisons.

Table 1: The primers used in this study.

Statistical analysis

Statistical analysis for this study was performed using SPSS 22.0 software. The means ± standard errors of the means were used to present the results. Statistical significance was determined by a p-value of less than 0.05. To analyze the data, one-way analysis of variance (ANOVA) along with the Tukey test was employed.

Ethical approval

All animal experiments were performed according to protocols and guidelines approved by the Institutional Animal Care and Use Committee of Zhejiang Institute of Freshwater Fisheries (Zhejiang Freshwater Fishery Environmental Monitoring Station), Huzhou, China (Animal Ethics no. 1067, March 6, 2019).

Result

Transcriptome sequencing, assembly and basic annotation

After excluding the low-quality reads and linkers, we prepared and analyzed clean reads (20.45–22.59 million reads, totaling 6.12–6.76 Gb) obtained from nine RNA libraries across the three experimental groups, with GC content ranging from 48.36%−49.97% and Q30 ranging from 94.24%−94.82% S1 Table). A total of 70054 unigenes (75.0 Mb in size) were generated by de novo assembly of the clean reads, with an average length of 791 bp, a maximum length of 35,540, an N50 length of 1284 bp, and a minimum length of 201 (S2 Table). The annotation of unigenes was conducted through BLAST alignment with six publicly available databases (COG, GO, KEGG, Nr, Pfam, and, Swiss-Prot). 22,5187 (32.14% of the total) unigenes were annotated and matched by at least one database. A total of 12,675 (18.09%), 36,146 (51.59%), 17,824 (25.44%), 24,995 (35.67%), 14,772 (21.08%), 31,185 (46.91%), and 21,902 (31.26%) unigenes were annotated in each database, respectively (Fig 2A).

Unigene annotation results.(A) Annotation information obtained from seven different databases, (B) Annotation information based on the Gene Ontology database, (C) The most common KEGG pathways.

Based on sequence homology, we identified the GO terms associated with each individual gene. A total of 36,146 individual genes were grouped into 58 gene ontology terms, which were categorized into BP, CC, and MF (biological process, cellular component, and molecular function). Within the BP category, the unigenes were further subdivided into 25 terms, including “cellular processes” (3408), “metabolic processes” (3290) and “single organism processes” (2988), which were the most important GO terms. Within the CC category, where “cell” (3065), “cell part” (3064), and “organelle” (2517) represent the main GO terms. Within the MF category, the unigenes were further subdivided into 12 terms. Most of the unigenes (3017 and 2287) were labeled with the terms “binding” and “catalytic activity” (Fig 2B). Using the KEGG database, 341 important signaling pathways involved in some developmental processes in animals were identified (Fig 2C). In addition, the “carbohydrate metabolism” (719), “translation”(1119), “signal transduction” (1434), “transport and catabolism”(1089), “endocrine system” (921), and “infectious diseases” (1502) contained significantly more unigenes than others.

Identification of DEGs and pathway enrichment analysis

The investigation focused on DEGs in various nitrite stress conditions. Compared with the CK group, a total of 640 DEGs were identified in the LC group, including 204 up-regulated genes and 436 down-regulated genes (Fig 3A). Among the top 20 KEGG pathways significantly enriched by these DEGs, the “Ribosome” (ko03010, adjusted p-value < 0.001) represented the most prominent enrichment, containing 96 DEGs (Fig 3B). The heatmap shows that low concentrations of nitrite mainly upregulated genes such as Tp63 and HSPB1, while downregulated genes including Act87E, PSAP, and MT-CO2 (Fig 3C).

Analysis of differential gene expression and KEGG enrichment of LC compared with the CK.(A) Differentially regulated genes under LC compared with CK. (B) KEGG enrichment analysis of the DEGs under LC and CK. (C) Heatmap analysis of differentially expressed genes between the LC and CK groups.

Compared with CK, the HC group identified 912 DEGs, including 403 up-regulated genes and 509 down-regulated genes (Fig 4A). 331 DEGs in the HC group belonged to the top 20 KEGG-enriched pathways, among which the “Ribosome”, “PPAR signaling pathway”, and “Staphylococcus aureus infection” pathways were extremely significantly enriched (adjusted p < 0.005) (Fig 4B). Heatmap analysis indicated that genes like Pino and sic25a42 were upregulated under high-concentration nitrite stress. Similarly, genes such as Act87E, PSAP, and MT-CO2 were also downregulated (Fig 4C).

Analysis of differential gene expression and KEGG enrichment of HC compared with the CK.(A) Differentially regulated genes under HC compared with CK. (B) KEGG enrichment analysis of the DEGs under HC and CK. (C) Heatmap analysis of differentially expressed genes between the HC and CK groups.

Compared with the LC group, the HC group identified 202 DEGs, including 80 up-regulated genes and 122 down-regulated genes (Fig 5A). Among the top 20 KEGG-enriched pathways, the “proximal tubule bicarbonate reclamation”, “glycerophospholipid metabolism”, and “citrate cycle” pathways were significantly enriched (adjusted p < 0.005) (Fig 5B). Heatmap analysis revealed that the differential genes of Macrobrachium rosenbergii under nitrite stress at different concentrations mainly included G6pdx, PISD, MFSD12, TIGD1, ANKRD66, and Gpcpd1, etc (Fig 5C). Additionally, we have exported the complete differential gene data for the CK vs. LC, CK vs. HC, and LC vs. HC groups from the database (S3, S4 and S5 Tables).

Analysis of differential gene expression and KEGG enrichment of LC compared with the HC.(A) Differentially regulated genes under LC compared with HC. (B) KEGG enrichment analysis of the DEGs under LC and HC. (C) Heatmap analysis of differentially expressed genes between the LC and HC groups.

M. rosenbergii miRNAs sequencing and characteristic analysis

The hepatic-pancreatic small RNAs of three treatment groups (CK, LC, HC) were sequenced using the Illumina Novaseq platform. The raw read counts were 50,019,908 for the CK group, 49,185,394 for the LC group, and 44,827,402 for the HC group. After data filtering, the high-quality clean reads obtained were 49,953,068 (99.87%) for the CK group, 49,121,064 (99.87%) for the LC group, and 44,774,198 (99.88%) for the HC group (S1 Table). The length distribution of small RNAs in each sample mainly concentrated at 21–23 nt, with the highest proportion at 22 nt. A total of 8,447 known miRNAs and 176 novel miRNAs were identified in the CK group; 10,675 known miRNAs and 206 novel miRNAs in the LC group; and 10,557 known miRNAs and 198 novel miRNAs in the HC group (S6 Table). The expression levels of different miRNAs varied, but miR-10465-y, miR-193-y, miR-263-x, novel-m0087-3p, novel-m0053-3p, novel-m0062-5p, and novel-m0080-5p showed high expression in all groups; while miR-145-y, miR-146-x, novel-m0022-5p, miR-202-x, novel-m0062-5p, miR-193-y, miR-263-x, and novel-m0185-5p showed low expression under different stress levels. All small RNAs were classified and annotated based on the miRBase and Rfam databases.

Differential expression analysis of miRNAs

Transcripts showing at least a 2-fold difference were classified as differentially expressed for miRNA analysis (|log2(fold change)| ≥ 1, p ≤ 0.05). The differences in the number of upregulated and downregulated miRNAs among different comparison groups reflect the concentration-dependent regulatory mechanism of nitrite stress intensity. In the LC vs CK group, the upregulation of 14 miRNAs may be related to energy metabolism and immune inflammation. In contrast, the significant downregulation of 17 miRNAs in the HC vs CK group may indicate severe energy damage and homeostasis threats. In the HC vs LC group, 13 upregulated and 20 downregulated miRNAs further reflect glycolipid metabolism reprogramming and multifunctional damage. Overall, these are closely related to the stress response of Macrobrachium rosenbergii. It is worth noting that key miRNAs, such as miR-145-y, show a consistent downregulation trend across group comparisons, suggesting their role as core regulatory factors. However, miR-m0062-5p exhibits opposite trends in different groups, indicating that different concentrations of nitrite have different regulatory effects on miRNAs. Furthermore, some newly predicted miRNAs, such as novel-m0087-3p, novel-m0053-3p, novel-m0062-5p, novel-m0080-5p, novel-m0022-5p and novel-m0185-5p, showed differential expression in inter-group comparisons, and exhibited a gradual down-regulation trend from the CK to HC groups, suggesting that they may play a role in the escalation of stress intensity (Fig 6A-C). Furthermore, we re-analyzed the data using a second reference gene, 18S rRNA, and the results were consistent with those in Fig 6, which further validates the reliability of the miRNA transcriptome (S1A-C Fig).

Effects of nitrite stress on miRNA levels of genes.(A) Differentially expressed miRNAs in M. rosenbergii (control vs high concentration), (B) Differentially expressed miRNAs in M. rosenbergii (control vs low concentration), (C) Differentially expressed miRNAs in M. rosenbergii (low concentration vs high concentration), (D) Related genes of linear change under Nitrous acid stress in M. rosenbergii. EF1-α gene was selected as an internal reference gene.

Validation of Selected miRNAs by qPCR

To validate the findings obtained through high-throughput sequencing, a subset of 16 miRNAs from various groups was randomly selected for qPCR analysis based on KEGG enrichment analysis. The expression levels of P53, SLC25a28, and RIO exhibited a significant upward trend with increasing nitrite concentrations, while HORMA showd a significant downward trend, which was consistent with the results of the high-throughput data(Fig 6D). Based on this concordance, it can be concluded that the sequencing results can be considered trustworthy and can serve as a basis for further analysis. Similarly, we performed qPCR analysis on genes using 18S rRNA as the reference gene. The results are consistent with those mentioned above (S1D Fig).

Target gene prediction and functional annotation of differentially expressed miRNAs

TargetScan (version 7.0) was used for target prediction. The software predicted differential expression of miRNA target genes in M. rosenbergii under nitrite stress, and the number of corresponding miRNA target genes was 13939. We then performed GO enrichment analysis. GO enrichment showed the similar results in the comparison of CK vs. LC (Fig 7A), CL vs. HC (Fig 7B) and LC vs. HC (Fig 7C). 51 GO-enriched miRNA target genes were identified from three aspects of CC, MF and BP. The BP classification of GO terms identified 22 groups for the target genes. The highest GO content was found in cellular process, single-organism process, and metabolic process. The classification of cellular components categorized the target genes into 19 groups. The most abundant GO terms in this classification were organelle part, macromolecular complex, and membrane part. In terms of molecular function, enhanced catalyst activity and signaling transmitter activity were the most prominent categories. According to the KEGG results, CK vs. LC (Fig 8A), CK vs. HC (Fig 8B) and LC vs. HC (Fig 8C) also observed the same trend, MAPK, apoptosis, autophagy, wnt, insulin and other signaling pathways are enriched.

Go enrichment analysis of the miRNA.(A) low concentration compared with control condition, (B) high concentration compared with control condition, (C) low concentration compared with high concentration.

Top of 20 KEGG enrichment analysis of the miRNA.(A) low concentration compared with control condition, (B) high concentration compared with control condition, (C) low concentration compared with high concentration.

Discussion

Nitrite stress represents a major environmental perturbation affecting the growth and development of fish worldwide [22] and is similarly prevalent in shrimp and crustaceans [23,24]. Aquatic organisms employ diverse defense strategies against nitrite-induced environmental stress, encompassing molecular, biochemical, and physiological alterations [20,25]. Typically, the primary metabolic shifts triggered by nitrite stress involve lipid, carbohydrate, and amino acid metabolism [26]. In this study, numerous stress-related genes and miRNAs were identified as significantly upregulated, suggesting a potential enhancemen of nitrite tolerance through the modulation of antioxidant capacity. We utilized integrated mRNA and miRNA transcriptomics via next-generation sequencing (NGS) to investigate the molecular changes in M. rosenbergii hepatopancreas under both mild (LC) and severe (HC) nitrite stress.

Our integrated analysis identified profound metabolic shifts as a core response. Compared to CK group, the DEGs and those related to energy and signal transduction in the LC group were primarily enriched in ribosome biogenesis (96 genes), protein digestion and absorption (14 genes), and the PPAR signaling pathway (11 genes). In this experiment, ribosomal protein RP-L10, RP-L13, RP-L17 and other genes were enriched and highly expressed. Ribosomes, as fundamental cellular machinery composed of RNA (rRNA) and proteins [27,28], function as the molecular factories responsible for protein synthesis, translating mRNA sequences into amino acid chains to build polypeptides [29,30]. For example, the empirical development of zebrafish is impacted by ribosomal protein L11 deficiency, leading to p53-dependent adaptive responses [31], and ribosomal protein can be utilized as an indicator of chronic bacterial heat stress in fish [32]. Furthermore, under acute environmental stress such as nitrite exposure, the significant fluctuations in the expression of ribosomal-related genes are more likely to reflect translational reprogramming or translational arrest in the organism. This serves as a defensive strategy through which the organism conserves energy expenditure while prioritizing the synthesis of critical stress-response proteins. Notably, DEGs from both the LC and HC treatment groups showed significantly enrichment in the PPAR pathway, with more metabolic pathways enriched in the HC group compared to LC. PPAR isoforms function as key transcription factors, crucially promoting lipid synthesis while inhibiting β-oxidation [33]. Zebrafish lacking the PPAR gene exhibit a severe loss of lipid synthesis function [34]. Another study has also demonstrated that the deletion of the PPARγ gene in zebrafish induces dysregulation of lipid metabolism-related protein kinases [35]. In the LC group, the miRNA-mRNA regulatory network tends to co-regulate pathways related to antioxidant enzymes and heat shock proteins, while moderately modulating metabolic rates. This reflects the organism’s active adaptive adjustments at the molecular level. In contrast, in the HC group, this regulatory balance is disrupted by severe molecular fluctuations. miRNA-mediated synergistic inhibition leads to a significant downregulation of mRNAs associated with the digestive and metabolic systems, such as those involved in amino acid metabolism and glycolysis, indicating functional impairment of the hepatopancreas under extreme toxic stress.

Nitrite stress can trigger emergency response in aquatic organisms, resulting in cell damage and death, and may also lead to growth abnormalities [36]. Different from the terrestrial environment, nitrite in the aquatic environment is more vulnerable to the impact of external factors and changes dramatically, posing a major threat to aquatic organisms. Upon stress exposure, the organism rapidly activates stress responses, altering the expression of critical signaling pathways to adapt. However, when stress exceeds tolerance limits, it can induce programmed cell death [37]. In the LC treatment group, miRNAs were primarily enriched in pathways governing dorsoventral axis formation, MAPK signaling, and Wnt signaling (p < 0.05), highlighting their regulatory importance under nitrite stress. Nitric acid stress is serious for the growth of aquatic organisms, and in severe cases, it will cause growth abnormalities and even death. Studies have shown that abnormal transduction of the Wnt/β-catenin signaling pathway leads to growth retardation in zebrafish [38], indicating the crucial role of the Wnt pathway in development. In this study, the Wnt signaling pathway was also significantly enriched (e.g., K02157, AXIN1; K02085, APC; K03083, GSK3), suggesting that it plays an important regulatory role in the growth and development of M. rosenbergii. The MAPK pathway is also a conserved signal transduction pathway, which consists of three major subtypes: ERK, JNK, and p38MAPK [39]. Primarily regulated by growth factors, the ERK and JNK pathways modulate cell cycle progression and growth [40,41]. In contrast, the p38MAPK pathway is chiefly activated by environmental stress, its high expression can mediate downstream apoptotic factors, promoting cell death [42,43]. Significant enrichment of ERK (K04371) within the MAPK pathway under LC stress suggests the organism possesses environmental adaptability at low nitrite levels, potentially achieved by modulating antioxidant capacity.

Although enrichment of dorsoventral axis formation, MAPK, and Wnt pathways was still observed in the HC treatment group, both the statistical significance and the enrichment factors were lower than in the LC group. This suggests that under high-concentration nitrite stress, the organism may redirect energy resources away from developmental processes like dorsoventral axis formation towards homeostasis maintenance and inflammation mitigation, as evidenced by the enriched lysozyme pathway [44]. Meanwhile, the enrichment of pathways including lysozyme [45], apoptosis [46], and autophagy [47] was enhanced in the HC group. These pathways have specific functions in repairing somatic cell damage: lysozyme is involved in immune defense; apoptosis is a normal physiological process for clearing damaged cells; and autophagy, which is closely related to apoptosis, can facilitate energy recovery by removing apoptotic cell debris, among other ways [48,49]. Collectively, these shifts in pathway enrichment patterns indicate adaptation to nitrite stress through energy reallocation and activation of specific damage repair mechanisms.

In general, comparing transcriptome and miRNA profiles under varying nitrite stress intensities reveals that Macrobrachium rosenbergii primarily counters nitrite stress through energy metabolism adjustments. It adapts to environmental stress-induced peripheral stimuli by altering metabolic patterns, and employs distinct signaling and metabolic pathways at different nitrite concentrations.

Conclusions

In summary, the physiological processes involved in the response of M. rosenbergii to nitrite stress are complex and diverse. This study delineates how M. rosenbergii modulates metabolic and signaling pathways to counteract nitrite toxicity. The discovery of novel miRNAs and stress-specific gene clusters advances our understanding of crustacean stress adaptation, offering biomarkers for selective breeding and sustainable aquaculture practices. At the same time, this study may provide a new data base to further elucidate the adaptation mechanisms of M. rosenbergii to nitrite stress.

Supporting information

S1 FigEffects of nitrite stress on miRNA levels of genes.(A) Differentially expressed miRNAs in M. rosenbergii (control vs high concentration), (B) Differentially expressed miRNAs in M. rosenbergii (control vs low concentration), (C) Differentially expressed miRNAs in M. rosenbergii (low concentration vs high concentration), (D) Related genes of linear change under Nitrous acid stress in M. rosenbergii.18s rRNA gene was selected as an internal reference gene.(DOCX)

S1 TableAssembly and annotation results of transcriptome.(DOCX)

S2 TableTranscriptome assembly quality statistic.(DOCX)

S3 TableDifferential gene data between CK and LC groups.(XLSX)

S4 TableDifferential gene data between CK and HC groups.(XLSX)

S5 TableDifferential gene data between LC and HC groups.(XLSX)

S6 TableSummary of miRNA sequencing results after fltering.(DOCX)