We promise that the present research complies with the applicable international and institutional policies relating to animal experiments. Meanwhile, all experimental methods involved in this study were approved by the Institutional Animal Care and Use Committee of Yantai University and performed in accordance with relevant guidelines and regulations. Additionally, all L. lota used in the present study were quickly transported to the laboratory and then were anesthetized with tricaine methanesulfonate (100 mg/L) before dissection.

L. lota (28.2 ± 1.03 g, 15.6 ± 0.90 cm) were collected from the aquafarm of the Irtysh River, Burqin, Xijiang Province (China) on 29 October 2019. Before poly (I:C) stimulation, all L. lota were placed temporarily in the same polyethylene aquariums (volume: 40 L) for a week to acclimatize, and the temperature, pH, ammonia nitrogen, nitrite, nitrate, and dissolved oxygen of water were maintained at 17°C, 7.9, 0.27 mg/L, 0.001 mg/L, 0.05 mg/L, and 5.0 mg/L, respectively. Normal feeding (3% of body weight; at 08:00 and 14:00) was performed throughout the experimental process. Meanwhile, we change the water every morning to ensure the stability of all water factors. After acclimation, all L. lota were randomly assigned to two separate aquariums (length × width × height: 72 cm × 53 cm × 44 cm), each of which was assigned 10 individuals (6 were used for subsequent experiments and the remaining 4 were used as backup), and the breeding density of L. lota was approximately 1,679.57 g/m³. We intraperitoneally injected 200 μL of PBS reagent (Beyotime Biotechnology, Shanghai, China) into individual L. lota in an aquarium as the control group, and injected 200 μL of poly (I:C) (Beyotime Biotechnology, Shanghai, China) at a concentration of 0.5 mg/mL (dissolved using PBS reagent) into individual L. lota in another aquarium as the treatment group. At 12 h and 48 h after injection, L. lota in the two groups were sacrificed and dissected to obtain the liver and spleen, and three individuals were randomly collected as biological replicates at each time point in each group. The obtained livers and spleens were immediately frozen with liquid nitrogen and then stored at −80°C. Finally, 24 tissue samples (2 tissues × 2 groups × 2 time nodes × 3 biological replicates) were obtained for subsequent RNA extraction. All experimental groups were named L_12PBS, S_12PBS, L_12po, S_12po, L_48PBS, S_48PBS, L48_po, and S_48po, respectively.

We used the Trizol Reagent Kit (Vazyme, Nanjing, China) to extract the total RNA of each tissue sample following the manufacturer’s protocol. The degradation or contamination degree, purity (OD260/280), concentration, and integrity of RNAs was detected using the agarose gel electrophoresis, Nanodrop (Implen, MUC, Germany), Qubit 2.0 (Thermo Fisher Scientific, MA, USA), and Agilent 2100 (Agilent Technologies, CA, USA), respectively.

The qualified total RNAs were used for library construction. We removed rRNA from the 1 μg of total RNA per sample and purified the remaining mRNA using RNA Purification Beads (Illumina, San Diego, CA, USA), and then cleaned all mRNAs three times using the Beads Binding Buffer (Illumina, San Diego, CA, USA). We used the NEBNext® UltraTM RNA Library Prep Kit for Illumina ® (NEB, USA) to construct 24 sequencing libraries after mRNA fragmentation. Specifically, we first reversely transcribed the mRNA into first-strand cDNA using random hexamers, and then added buffer, dNTPs, and DNA polymerase I to synthesize the second-strand cDNA. The synthesized double-stranded cDNAs were purified using AMPure XP beads (Illumina, San Diego, CA, USA), and terminal repair, A-tail addition, and sequencing connector were then performed. Subsequently, we added A-tails and adapters to the double-stranded cDNAs, respectively. Then, AMPure XP Beads were applied to select the fragment size, and these fragments with appropriate size were amplified by PCR to obtain the cDNA libraries. The constructed cDNA libraries were quantified using Qubit 2.0 (Thermo Fisher Scientific, MA, USA) and then diluted to 1 ng/µL. Finally, the cDNA libraries qualified by the Agilent 2100 (Agilent Technologies, CA, USA) were sequenced on an Illumina HiSeq 2000 platform across one lane with 150 bp paired-end.

Trimmomatic software (version 0.36, 27) was used to eliminate low-quality raw RNA-seq reads, mainly referring to these reads that contained sequencing adapters, unknown nucleotides ratio > 10%, and quality scores < 20. Hisat software (version 2.0.4; 28) was used to compare the high-quality clean RNA-seq reads to the published L. lota whole genome sequences (29). The regional distribution and density distribution of clean RNA-seq reads on the genome sequences were then analyzed. Meanwhile, we use rMATS (version 3.2.5; 30) software to predict the classification, number, and structure of alternative splicing (AS) events.

To identify the transcriptional responses of liver and spleen in L. lota to the poly (I:C) stimulation, we first quantified the gene expression levels of each sample using the HTSeq (version 2.0; 31) software and described them using FPKM (expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced). Then, we identified the differentially expressed genes (DEGs) among different experimental groups using DESeq software (version 2.0, 32), and false discovery rate (FDR) adjusted p-value < 0.005 and fold change (FC) ≥ 2 (which corresponds to |log₂FC| ≥ 1) were used as the filtering thresholds. A Venn diagram was applied to visualize the number of DEGs. Furthermore, we determined the clustering relationship of DEG expression level in four experimental groups according to the FPKMs. Finally, we attempted to analyze the biological functions of these DEGs and their products based on the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) annotation. The hyper-geometric distribution corrected p-value < 0.05 was taken as the criterion of significant enrichment of the GO term and KEGG pathway.

In the present study, quantitative reverse transcription PCR (qRT-PCR) was applied to verify the reliability of transcriptome data. According to the functional annotation information, a total of eight immune-related DEGs (four in the liver and four in the spleen) were randomly selected for qRT-PCR experiments. Meanwhile, Ribosomal protein S29 (rps29) and Ribosomal protein L26 (rpl26) were used as reference genes for standardization (22). Primer Premier 6.0 was applied to design the reference gene-specific and DEG-specific primers ( Table 1 ). Twenty-four cDNA samples from 12 experimental individuals were diluted 25 times using nuclease-free water according to standard curves and then used as templates for qRT-PCR. The qRT-PCR was conducted on the StepOne Plus Real-Time PCR system (ABI, USA) according to the manufacturer’s instructions of the TaKaRa TB Green Premix Ex Taq (Til RAaseH Plus, RR420A). A 20-μL reaction system was formulated, including 2 μL of diluted cDNA template, 0.4 μL of forward primer (10 μM), 0.4 μL of reverse primer (10 μM), 0.4 μL of ROX Reference Dye (50×), and 6.8 μL of RNase-free water. The qRT-PCR cycling conditions were as follows: one cycle of 95 °C for 30 s, followed by 40 cycles of 5 s at 95°C, 30 s at 60°C, and then enter the dissociation stage. Three experimental triplicates were performed for each qRT-PCR cycling condition to ensure the accuracy of qRT-PCR results. The relative expression levels of eight DEGs were calculated by the 2^−ΔΔCT method (ΔCT = CT_DEG − CT_{internal gene}, ΔΔCT = ΔCT_{treatment group} – ΔCT_{control group}).

All raw RNA-seq reads have been submitted to the NCBI sequence read archive under the BioProject PRJNA874751. A total of 228.08 Gb of clean RNA-seq reads were captured after removing the low-quality raw RNA-seq reads, and the summary of RNA-seq reads is shown in Supplementary File 1 . All clean RNA-seq reads from 24 samples were compared to the reference genome ( Supplementary File 2 ) and results showed that 90.97% of the clean RNA-seq reads can be located on the reference genome. Meanwhile, the percentages of clean RNA-seq reads with multiple mapped locations and uniquely mapped locations on the reference genome are 7.14% and 83.83%, respectively. Furthermore, although many clean RNA-seq reads are mapped to the exon region of reference genome, a small number of reads are mapped to the intergenic region, and a minimal proportion of clean RNA-seq reads were compared to intron regions ( Supplementary File 3 ). Meanwhile, there was a positive correlation between the distribution density of clean RNA-seq reads and chromosome length ( Supplementary File 4 ). Additionally, a total of 20,646 predicted genes were successfully annotated on at least one protein database. Of all annotated predicted genes, 3,432, 14,369, 13,246, 17,187, 20,576, 18,218, 13,388, and 11,282 predicted genes had significantly matched with the sequences in the AnimalTFDB, GO, KEGG, KOG, NR, Pfam, Swiss-Prot, and TrEMBL databases, respectively ( Figure 1 ).

A total of 3,347 and 3,361 AS events were predicted based on “junction count only-” and “reads on target and junction counts-” quantitative methods, respectively. Meanwhile, a total of 198 and 202 different AS events were predicted based on “junction count only-” and “reads on target and junction counts-” quantitative methods, respectively. Among the five common AS events (namely, skipped exon [SE], alternative 5′ splice site [A5SS], alternative 3′ splice site [A3SS], mutually exclusive exons [MXE], and retained intron [RI]), SE and MXE were found to be the most abundant AS events ( Table 2 ).

In the present study, we have identified the poly (I:C)-stimulated DEGs in liver (L_12po-vs-L_12PBS, L_48po-vs-L_48PBS) and spleen (S_12po-vs-S_12PBS, and S_48po-vs-S_48PBS) of L. lota ( Supplementary File 5 ). Results showed that the number of DEGs increased with the increase in infiltration time of liver (2,186 at 12 h and 3,123 at 48 h; Figure 2A ) and spleen (1,542 at 12 h and 2,622 at 48 h; Figure 2B ) in poly (I:C). Compared with the PBS treatments, a total of 1,231 (56.31%) and 1,776 (56.87%) upregulated genes were identified from the liver exposed to poly (I:C) at 12 h and 48 h, respectively. Meanwhile, a total of 769 (49.87%) and 1,502 (57.28%) upregulated genes were identified from the spleen exposed to poly (I:C) at 12 h and 48 h, respectively. Additionally, cluster analysis showed that the expression levels of DEGs from the same tissue and the same treatment were similar, and the similarity of DEG expression levels in the same tissue was higher than that between tissues ( Figure 3 ).

We found that a substantial overlap of 300 DEGs existed in four experimental pairs (L_12po-vs-L_12PBS, L_48po-vs-L_48PBS, S_12po-vs-S_12PBS, and S_48po-vs-S_48PBS) ( Figure 4 ). Considering that these genes contributed more to the responses of L. lota to poly (I:C), we obtained annotation information of 300 DEGs ( Supplementary File 6 ). Not unexpectedly, there are abundant immune-related genes (such as Nuclear valosin-containing protein [NVLp], Tumor necrosis factor [TNF], Zinc-binding protein [ZBP], Antigen peptide transporter (APT), Calcium/calmodulin-dependent protein kinase type 1G [Camk1g], and Interferon regulatory factor 4 [IRF-4]) in L. lota that were ultimately predicted (the yellow marked part of Supplementary File 6 ).

Unsurprisingly, poly (I:C) injection significantly affected the biological processes of L. lota, and this phenomenon had time and tissue specificity. Firstly, we matched the up- and downreg ulated genes in the liver and spleen into the GO database and then obtained the GO annotation information ( Supplementary File 7 ; Figure 5 ). At 12 h after poly (I:C) injection, upregulated genes in the liver and spleen were significantly enriched in 79 and 156 GO terms (corrected p-value < 0.05), respectively. Meanwhile, downregulated genes in the liver and spleen were significantly enriched in 78 and 74 GO terms (corrected p-value < 0.05), respectively. At 48 h after poly (I:C) injection, upregulated genes in the liver and spleen were significantly enriched in 162 and 73 GO terms (corrected p-value < 0.05), respectively. Meanwhile, downregulation in the liver and spleen were significantly enriched in 12 and 24 GO terms (corrected p-value < 0.05), respectively. The terms of up- or down regulated genes are related to protein regulation, amino acid regulation, energy regulation, cellular homeostasis, immune cell integrity and ion transport.

We recorded the networks of molecular interactions in the cells and the variants specific to L. lota by mapping the up- and downreg ulated genes in the liver and spleen to the KEGG database, respectively ( Supplementary File 8 ; Figure 6 ). At 12 h after poly (I:C) injection, upregulated genes in the liver and spleen were significantly enriched in the 14 (namely, Herpes simplex infection pathway, Proteasome pathway, Apoptosis pathway, Jak-STAT signaling pathway, p53 signaling pathway, RIG-I-like receptor signaling pathway, Toll-like receptor signaling pathway, Cytosolic DNA-sensing pathway, Glycerophospholipid metabolism, Ubiquitin-mediated proteolysis, Ether lipid metabolism pathway, Cytokine–cytokine receptor interaction pathway, alpha-Linolenic acid metabolism pathway, and NOD-like receptor signaling pathway) and 11 (namely, RIG-I-like receptor signaling pathway, Herpes simplex infection pathway, Apoptosis pathway, Toll-like receptor signaling pathway, p53 signaling pathway, NOD-like receptor signaling pathway, Jak-STAT signaling pathway, Cytosolic DNA-sensing pathway, Cytokine–cytokine receptor interaction pathway, Proteasome pathway, and Drug metabolism - other enzymes pathway) metabolic pathways, respectively (corrected p-value < 0.05). Meanwhile, downregulated genes in the liver and spleen were significantly enriched in the seven (namely, Phosphatidylinositol signaling system pathway, Adherens junction pathway, Focal adhesion pathway, Notch signaling pathway, Glycerolipid metabolism pathway, Inositol phosphate metabolism pathway, and TGF-beta signaling pathway) and five (namely, Focal adhesion pathway, ECM–receptor interaction pathway, Vascular smooth muscle contraction pathway, Regulation of actin cytoskeleton pathway, and Cytokine–cytokine receptor interaction pathway) metabolic pathways, respectively (corrected p-value < 0.05). However, the number and type of metabolic pathways significantly enriched at 48 h of poly (I:C) stimulation varied, which was present in both tissues. Results showed that upregulated genes in the liver and spleen were significantly enriched in the 13 (namely, Ribosome pathway, Proteasome pathway, Spliceosome pathway, Ribosome biogenesis in eukaryotes pathway, Pyrimidine metabolism pathway, RNA transport pathway, DNA replication pathway, Herpes simplex infection pathway, RNA polymerase pathway, Purine metabolism pathway, Cytosolic DNA-sensing pathway, RIG-I-like receptor signaling pathway, and Protein processing in endoplasmic reticulum pathway) and 5 (namely, Proteasome pathway, Ribosome biogenesis in eukaryotes pathway, Herpes simplex infection pathway, RIG-I-like receptor signaling pathway, and Pyrimidine metabolism pathway) metabolic pathways, respectively (corrected p-value < 0.05). Meanwhile, downregulated genes in the liver and spleen were significantly enriched in the two (namely, ECM–receptor interaction pathway and ABC transporters pathway) and six (included Focal adhesion pathway, Phosphatidylinositol signaling system pathway, ECM–receptor interaction pathway, Pentose phosphate pathway, Vascular smooth muscle contraction pathway, and Regulation of actin cytoskeleton pathway) metabolic pathways, respectively (corrected p-value < 0.05).

The expression trend of eight genes was concordant based on the qRT-PCR data and the RNA-seq data, which meant that the sequencing results were reliable. The different genes showed different expression trends after the poly (I:C) challenge ( Figure 7 ). Detailed casp3, hsp90a1, and lgp2 in the liver were upregulated at 12 h and 48 h after poly (I:C) stimulation, while met was downregulated. In the spleen, csf1r, ccnd2, and flt3 were downregulated at 12 h and 48 h after poly (I:C) stimulation. Meanwhile, pmp22 in the spleen was downregulated at 12 h after poly (I:C) stimulation, but upregulated at 48 h after poly (I:C) stimulation.

The present clean RNA-seq reads covers more than 90% of the L. lota genome, and the percentage of clean RNA-seq reads with multiple locations on the genome were less than 10%. This means that our sequenced reads were uncontaminated. Meanwhile, a large proportion of reads could be compared to the exon region and a small proportion to the intergenic region, which provided evidence for the integrity of the L. lota genome annotation information. Residues from pre-mRNA and intron retention events during AS processes resulted in a small number of reads being compared to the intron region. Additionally, many protein-coding genes (20,646/21,664; 95.30%; 29) can be functionally annotated in at least one protein database. In conclusion, our sequencing reads can provide a basic resource for studying the regulatory mechanisms of L. lota to poly (I:C).