For the study, a total of 468 raw RNA-Seq datasets were obtained from the National Center for Biotechnology Information (NCBI) database (https://www.ncbi.nlm.nih.gov/) for common carp, covering >9.7 billion transcript reads. The dataset encompassed 42 bioprojects from 23 institutions across 7 countries. These included data across 28 distinct tissues from different varieties of common carp (viz. koi, haematopterus, specularis, color, huanghe, wuyuanesis, singuonensis, jian), the details of which are provided in Supplementary Table S1. All the bioinformatics analysis were tissue-specific, i.e., performed separately for each of the 28 tissues.

The raw reads obtained from the NCBI were first visually assessed for the quality using FastQC tool ver. 0.11.8 (Schmieder and Edwards, 2011) followed by elimination of the adaptor sequences and low-quality reads using Trimmomatic software ver. 0.39 (Bolger et al., 2014). To ensure non-interference of polyA tail, reads were trimmed off with primer, followed by automatic transcriptome assembly. Moreover, lncRNAs are identified based on alignment with reference sequence which further ensures exclusion of such tails as they are not genetically coded. They are added post-transcriptionally for protection/shelf life regulating translational efficiency. Common carp reference genome and annotation files were downloaded from NCBI (https://www.ncbi.nlm.nih.gov/assembly/GCF_000951615.1/) (Xu et al., 2014). Index files of the reference genome were generated using the HISAT2-build function of HISAT2 version 2.2.0 (Kim et al., 2015; Pertea et al., 2016). Sam files were aligned and converted to binary bam files using Samtools software version 1.9 (Li et al., 2009). StringTie software version 2.1.4 was used for transcriptome assembly of the individual bam files and generate Gene Transfer Format (gtf) files for each transcriptome reads (Pertea et al., 2015; Pertea et al., 2016). Finally, the StringTie-merge feature was used to combine tissue-specific files and generate a single gtf file per tissue (Banerjee et al., 2021).

To identify potential lncRNA transcripts, the FASTA sequences that corresponded to every transcript within the combined assembly file were obtained by the gffread program version 0.12.3 using the respective reference genome (Pertea and Pertea, 2020). Owing to the longer size of lncRNAs, transcripts <200 base pairs were removed using perl scripts. ORFs were predicted using ORFPredictor, and those exceeding 300 nucleotides were eliminated (Min et al., 2005; Wang et al., 2021). The coding potential of the transcripts were assessed using CPC2 ver 1.0.1 (Kang et al., 2017), and PLEK ver 1.2 (Li et al., 2014), eliminating the coding RNAs. Non-coding RNAs were discovered through a BlastN search against RNACentral (https://rnacentral.org/) and transcripts with at least 95% identity were excluded (Shumayla et al., 2017). The remaining transcripts might contain small classes of non-coding RNAs, such as mRNA, tRNA, rRNA, miRNA, and snRNA. These underwent Blastp search (Mount, 2007) against the Pfam (http://pfam.xfam.org/) database (Altschul et al., 1990) and non-redundant database (https://www.ncbi.nlm.nih.gov/protein/) to eliminate recognized protein-coding RNAs according to the C. carpio annotation. The remaining transcripts were considered to be the potential lncRNAs for further analysis. The FPKM values for these transcripts were also calculated. The pipeline for identification of lncRNA in C. carpio is delineated in Figure 1.

The gffcompare software was used to categorize lncRNAs into different groups based on their position relative to protein-coding genes, namely, lincRNA (intergenic lncRNAs), labeled as “u”; exonic lncRNAs (intersecting with protein-coding exons), labeled as “x”; and intronic lncRNAs (existing in introns without sharing sequences with exons), labeled as “i” (Leiva et al., 2020). Further, the protein-coding genes were eliminated through an evaluation and elimination process.

Compared to protein-coding mRNAs in different species, lncRNAs have lower conservation (Derrien et al., 2012). The conservation of lncRNAs in C. carpio with zebrafish, rainbow trout, and large yellow croaker was assessed using BLAST, with a cutoff E-value < 1e⁻⁶ (Zhang et al., 2022). Zebrafish data were obtained from the ZFLNC database (https://www.biochen.org/zflnc/) (Hu et al., 2018), while information on rainbow trout and large yellow croaker came from previous lncRNA profiling studies (Al-Tobasei et al., 2016; Zhang et al., 2020).

A total of 9,775,580,802 raw reads generated by Illumina HiSeq platform from 28 different tissues was collected from 468 RNA-seq datasets from the NCBI SRA database. After discarding adaptor sequences and low-quality reads, 9,385,180,448 clean reads (94.88%) were obtained. Approximately 80% of the clean reads from the 468 samples were aligned to the reference genome of the common carp using HISAT2 (Supplementary Table S2). A total of 268,353 transcripts for tissue-specific lncRNAs were generated using StringTie-merge module.

On comparing the four GTF files to the existing C. carpio annotation file using GffCompare, allowing transcript annotation based on genomic location relative to known genes, 300,218 transcripts were found to lie in the u (unknown intergenic), x (genic antisense) and i (intronic) classes (Figure 3A). A total of 33,990 putative lncRNAs were discovered in C. carpio after applying filtering criteria, namely, eliminating nucleotide sequences <200 nucleotides, ORF length >300 nucleotides, CPC2 score >0.5, housekeeping RNAs with >95% identity and transcripts similar to protein families or genes. Majority of the sequences showed GC% between 20 and 40 (Figure 3B), while the most abundant sequence length was 200–400 bp, followed by 400–600 bp (Figure 3C). Almost 56.13% of sequences had one exon, followed by 40.36% having two exons (Figure 3D). The distribution of these lncRNAs across different tissues shows kidney tissues to exhibit the highest number of identified lncRNAs, i.e., 8,003 (23.54%) (Figure 3E). We found three lncRNAs, TCONS_00121934, TCONS_00177318, and TCONS_00328247, abundantly expressed in all 28 tissue types we studied. Figure 3E shows the chromosome-wise distribution of lncRNAs in common carp. The calculate FPKM value across all 28 tissues and distinct levels of expression were discerned, as revealed by the comprehensive analysis of the bioprojects. The study unveiled a diverse range of average expression values for lncRNAs across multiple tissues as provided in Supplementary Table S3.

Out of the total 33,990 lncRNAs, the majority (65.30%) belonged to intergenic class of lncRNAs (u), followed by intronic lncRNAs (i) and exonic lncRNAs (x) which were 27.28% and 7.39%, respectively (Figure 3A). These findings suggest that the majority of lncRNAs (i.e., intergenic or lincRNAs), do not overlap with protein-coding genes. To investigate the conservation of lncRNAs among various species, these putative lncRNAs discovered in C. carpio were compared to those of zebrafish, rainbow trout, and large yellow croaker through a blast analysis. The results showed that only a small number of lncRNAs, i.e., 3,484 (10.25%) in zebrafish (cyprinid), 138 (0.4%) in rainbow trout (salmonid), and 43 (0.12%) in large yellow croaker (sciaenidae) were conserved (Supplementary Table S4). This low level of similarity may be due to poor conservation of lncRNAs across species and tissues.

It was observed from the analysis that the most common carp lncRNAs had a GC content ranging from 0% to 66.4%, but the majority (99.70%) were within 20%–60%, with an average GC content of 37.88% (Figure 3B). The length of lncRNAs in common carp varied from 200 to 15,799 nucleotides (Figure 3C), with 72.17% being between 0 and 999 nucleotides and 27.82% between 1,000 and 1,999 nucleotides. The average length of the lncRNAs was found to be 837 nucleotides, which is shorter than the length of protein-coding genes (2.914 kb). Common carp lncRNAs had 1–7 exons on average (1.48 exons). The majority (96.49%) were single or double exon types, with an average length of 849.82 nucleotides (Figure 3D). The chromosomal distribution of the identified lncRNA in common carp is shown in (Figure 3F), with chromosome 38 having the highest count of 1,105 lncRNAs and chromosome 19 having the lowest count of 252. This distribution was visualized using Circos software (Figure 4).

The identified lncRNAs were analyzed for their potential target genes based on their proximity to protein-coding genes. Additionally, the sequences of the lncRNAs were compared to miRNAs from miRBase using the BlastN program to determine if they function as miRNA precursors. Results showed that 9 pre-miRNAs matched with 19 distinct lncRNAs in common carp, indicating the potential for these lncRNAs to produce mature miRNAs at matching accuracy cut-off of ≥90% (Supplementary Table S5). Few lncRNAs were discovered to have a stable hairpin structure, indicating the presence of a miRNA precursor. In our study, TCONS_00417363 lncRNA (dark green), which contained the miRNA precursor ami-mir-133a-2 (red) was identified by Vienna RNA package and RNAfold program (Supplementary Figure S1A).

LncRNAs play a crucial role in regulating gene expression by acting as competitive endogenous RNAs (ceRNAs), which capture miRNAs and prevent them from binding to their target mRNAs. This interaction between lncRNAs and miRNAs significantly affects gene expression and various biological processes. The miRNAs also regulate gene expression by binding to the 3′ untranslated region (UTR) of target mRNAs, leading to mRNA degradation or translational repression. These miRNA-mediated gene expression regulations are essential for the normal development and function of fish organs and tissues. Using the psRNAtarget server, 18,977 interactions (lncRNA-miRNA) involving 2,837 distinct miRNAs and 4,782 distinct lncRNAs in common carp were identified. Using the RNAfold program with the Vienna RNA package, a visualization of the secondary structure of lncRNA TCONS_00239471 (displayed in dark green), along with the locations of its binding sites for miRNAs ccr-miR-338 (depicted in red) and mir-338-y (depicted in blue) were generated (Supplementary Figure S1B). In the case of miRNA-mRNA, a total of 8,079 interactions, which involved 1,900 unique miRNAs and 3,718 unique mRNAs were revealed. The Cytoscape software was utilized to construct individual lncRNA-miRNA and miRNA-mRNA interaction networks after that merged both the network in a single network and visualizing the whole lncRNA-miRNA-mRNA interaction network. From this interaction analysis, we observe that mir-8499-y target 2 lncRNAs of TCONS_00413548, TCONS_00007855 and six mRNAs these are lcl|LHQP01046800.1_mrna_46126, lcl|LHQP01016561.1_mrna_24063, lcl|LHQP01011828.1_mrna_17740, lcl|LHQP01021498.1_mrna_29844, lcl|LHQP01064132.1_mrna_49370 and lcl|LHQP01015883.1_mrna_23319 (Supplementary Figure S2).

The identification of circRNAs in common carp involved four main step: i) quality check of 9,775,580,802 raw RNA-seq reads from NCBI ii) final clean reads aligned to the common carp reference genome using BWA-mem -T 19 ii) generation of 468 SAM files for each sample, followed by merging for each of the 28 tissues. iv) Submission of these merged files to the widely used command-line tool, CIRI2. This resulted into 22,854 distinct potential circRNAs (Supplementary Table S6). In this study, we have not only identified tissue-specific common carp circRNAs but have also uncovered striking differences in their abundance across tissues. A total of 11,621 circRNAs (50.84%) were derived from intergenic regions, while only 2,575 (11.26%) were generated by introns (Figure 5A). Out of the identified circRNAs, 8,655 (37.87%) originated from exons of protein-coding genes, indicating they were exonic circRNAs. These exonic circRNAs had both back-splice sites aligned with known exonic boundaries. A visual representation of the length distribution of circRNAs reflected that majority of circRNAs were over 1,600 nucleotides long (Figure 5B). Specifically, embryonic tissue showed the abundance of circRNAs (8635) (Figure 5C). Chromosome 38 was seen to have highest number of circRNA (745) (Figure 5D). The chromosome-wise circos-map of circular RNAs in common carp is summarized in Figure 4 (outermost dark purple circle).

CircRNA-miRNA-mRNA interactions have been identified as an important regulatory mechanism in gene expression which is also called the ceRNAs network, where circRNAs act as a sponge for miRNAs, thereby inhibiting the degradation of target mRNAs by miRNAs. By using TargetFinder software we get a total of 970,159 interactions (circRNA-miRNA) where 5,906 unique miRNAs were found targeting 15,731 circRNAs in the tissues under study, except the pituitary tissue. As a result, 10,484 interaction (miRNA-mRNA) was found from the psRNATarget server analysis where 4,524 distinct mRNAs were targeted by 2,871 unique miRNAs of common carp. An entire circRNA–miRNA-mRNA interaction network was delineated by Cytoscape (Supplementary Figure S3). The figures show that miRNA, mir-6627-y targeted a total of seven mRNAs, miRNA mir-6651-x targeted nine circRNAs and two mRNAs, while another miRNA, mir-7371-x targeted a total of 13 circRNAs.

In this study, we utilized GO and KEGG annotations to gain understanding of the functions of lncRNAs/circRNAs, based on the hypothesis that their functions may be linked to those of their parent genes. An analysis of GO categories and KEGG pathways was done on 3,718 host genes of 33,990 lncRNAs and 4,524 host genes of 22,854 circRNAs to investigate their possible roles in common carp. For the lncRNAs in cellular component category, the top three largest groups were the nucleus (7.31%), plasma membrane (5.85%), and membrane (2.62%); for the biological process: anatomical structure development (10.02%), signaling (8.02%), and protein modification process (5.69%) and for molecular function: transferase activity (7.29%), catalytic activity (5.38%), and hydrolase activity (3.78%). The significantly enriched KEGG pathways were mainly necroptosis, NOD-like receptor signaling pathway, hypertrophic cardiomyopathy, small cell lung cancer, MAPK signaling pathway, pathways of neurodegeneration, and axon guidance (Supplementary Tables S7, S8). The analyses of GO and KEGG indicate a strong potential for lncRNAs and circRNAs to play various roles in biological processes within the common carp fish.

A web-based database, named CCncRNAdb for common carp ncRNAs is accessible at http://backlin.cabgrid.res.in/ccncrnadb/index.php, which contains information on 33,990 lncRNAs and 22,854 circRNAs, including their characterization, and interactions with miRNAs and mRNAs. CCncRNAdb has six tabs, namely, Home, lncRNA, Interactions, circRNA, Download, and Teams (Figure 6). The main features of this database are: CCncRNAdb offers extensive information on tissue-specific lncRNAs and circRNAs, such as their chromosome locations, sequence length, and coding potential. To obtain details regarding the lncRNA/circRNA-miRNA-mRNA interaction network, users can visit the “Interaction” tab. The results can be downloaded directly for all the tissues.