With the aim to perform a comprehensive study, a total of 53 raw RNA-Seq datasets (>1.2 billion reads) of black pepper were downloaded from the NCBI database. The dataset covered eight BioProjects, seven institutes across four countries, and six tissues (flower, fruit, leaf, panicle, root, and stem) from six black pepper cultivars viz. Reyin-1, Bragantina, Thottumuriyan, IPN No. LK-0-WU-0014181, panniyur-1, and Genotype 4226 ( Supplementary Table 1 ).

Raw reads obtained from NCBI were subjected to quality check using FastQC tool ver. 0.11.8 (Andrews, 2010), which helps us visualize various quality parameters like per base sequence quality, per base sequence content, presence of adaptor sequences, etc. The Trimmomatic ver. 0.39 (Bolger et al., 2014) software was then used to trim out the probable contaminants like adaptor sequences and low-quality reads with a Phred score of less than 30. The reference genome and annotation file of black pepper were downloaded from the GCGI (Group of Cotton Genetic Improvement, http://cotton.hszau.edu.cn/EN/index.php) website (Hu et al., 2019). The HISAT2-build command from HISAT2 ver. 2.2.0 (Kim et al., 2015; Pertea et al., 2015; Pertea et al., 2016) software was used to index the reference genome with splice sites and exons information retrieved from the annotation file. Index files were then used for aligning individual clean reads. Sam files obtained after the reads alignment were converted into binary bam files using Samtools software ver. 1.9 (Li et al., 2009), then the transcriptome assembly of the individual bam files was performed using StringTie ver 2.1.4 (Pertea et al., 2015) software to give out gtf files for each of the reads. Individual files corresponding to each tissue were then merged using StringTie–merge to get a single gtf (Gene Transfer Format) file for each tissue.

Identification of lncRNA candidate transcripts from the assembled transcripts involved various steps as shown in Figure 1 . First the fasta sequences corresponding to each transcript in the merged assembly file were extracted from their respective reference genome fasta files using the gffread program ver 0.12.3. As lncRNAs are RNA transcripts longer than 200 base pairs, using in-house Perl scripts, transcripts shorter than 200 bp were removed. Studies have shown that lncRNAs in general have a lower quality shorter ORFs than the protein-coding mRNAs. ORFPredictor (Min et al., 2005) was used to predict the ORFs in each transcript and those having ORFs longer than 300 nucleotides were removed. Binary classifiers like CPC2 ver 1.0.1 (Kang et al., 2017) and PLEK (Li et al., 2014) were used to classify the remaining sequences into coding or non-coding. Transcripts showing coding labels in all of them were discarded and kept as coding transcripts. The remaining transcripts may contain some housekeeping RNAs like rRNAs, tRNAs, snoRNAs, and other ncRNAs. So to get novel ncRNAs we did a BlastN search against the databases like SILVA database (https://www.arb-silva.de/download/arb-files/release138.1), GtRNAdb (http://gtrnadb.ucsc.edu/release18.1), and RNACentral (https://rnacentral.org/release16.0) and removed the transcripts showing >=95% identity (Sharma et al., 2017). Transcripts matching any reported proteins or protein families were identified and removed using Blast against NCBI-nr protein and Pfam (http://pfam.xfam.org/release33.1) databases (Altschul et al., 1990).

LncRNAs are poorly conserved compared to protein-coding mRNAs across species. Yet to check the conservativeness of black pepper’s lncRNAs, reported lncRNAs of 38 plant species available at the CANTATAdb ver. 2.0 (http://cantata.amu.edu.pl/) (Szcześniak et al., 2019) database were downloaded and a local BlastN search was performed against our identified black pepper lncRNAs as a query with 10^-20 e-value as the cutoff.

A web-resource, a Black pepper ncRNA database BPncRDB was created using the three-tier architecture, viz., client, middle, and database tiers that house all the results of this study related to the lncRNAs, circRNAs and their interactions with the miRNAs. The database was developed in MySQL database (https://www.mysql.com/) while the web-interface was prepared in PHP (https://www.php.net/) and HTML while designing was done using CSS and made dynamic using JavaScript. It was hosted on Apache server (https://httpd.apache.org/). XAMPP framework was used to design and deploy the webpage. The user can retrieve data as: (a) request from user to webserver, (b) query sent to MySQL database, (c) response generated by database and sent to web-interface, and (d) response of web-server to user. BPncRDB includes information of tissue-specific lncRNAs and circRNAs and their relationship with miRNAs, interaction between lncRNA-miRNA and miRNA-circRNA, etc.

The library was sequenced using the Illumina HiSeq X platform and 1,278,079,048 raw reads in 16 samples were collected. After discarding the Illumina platform’s adaptor sequences and other low-quality reads using the Trimmomatic software, we obtained 1,251,005,293 (97.88%) clean reads ( Supplementary Table 2 ). The trimmed clean reads were aligned to black pepper’s reference genome using HISAT2 and approximately 70%–94% of reads were mapped across the 53 samples ( Supplementary Table 2 ). StringTie software was then used to assemble the mapped reads of individual samples with respect to the reference annotation file of black pepper. Individual assembly files for each tissue were merged using the StringTie-merge module in order to get tissue-wise lncRNAs, and 294,777 transcripts were obtained.

A stringent filtering pipeline for the assembled transcripts was developed and used for the identification of those transcripts which fulfill the criteria of lncRNAs ( Figure 1 ). GffCompare program was used to compare the four GTF files with an annotation file of black pepper and annotate the transcripts corresponding to their location on the genome with respect to the known genes (Pertea and Pertea, 2020). Out of the 15 class codes, i (intronic), u (intergenic), and x (natural antisense transcripts) represent the most probable non-coding transcripts and were selected for the downstream analysis, and the rest were discarded. Gffread was used to extract fasta sequences corresponding to the class codes and a total of 41090 sequences were found. In-house developed Perl scripts were used to remove sequences smaller than 200 nucleotides. To remove the potentially coding transcripts, ORFpredictor was used to find out the ORFs in each transcript and those having an ORF length >300 nucleotides were discarded. The coding probability of the remaining transcripts was calculated using CPC2 and PLEK software and taking the intersection of the results. Transcripts with CPC2 score >0.5 and predicted coding by PLEK were considered coding and discarded. Housekeeping RNAs were removed by BlastN against ncRNA databases, viz., Silvadb, gtRNAdb, and RNACentral and having percent identity >95%. Transcripts having similarity matching with any of the protein families in Pfam or genes in the NCBI-nr database were removed using BlastX (e-value 10^-3). Finally, we identified 6406 novel black pepper long non-coding RNAs in black pepper, out of which 1115, 2621, 2727, 828, 1214, and 1003 were expressed in flower, fruit, leaf, panicle, root, and stem tissues, respectively ( Figure 3A ; Supplementary Table 3 ).

The study presents the first ever identification and characterization of black pepper lncRNAs to understand the functional importance of lncRNAs via various mechanisms affecting the crop’s yield and stress responses. In consistency with the previous lncRNA studies most (~80%) of the identified black pepper lncRNAs are intergenic in nature ( Figure 3B ). The length of the identified lncRNAs was distributed in the range of 200–10667 nucleotides which is much shorter than mRNAs which ranges up to 16.6 mega bases in length. More than 80% of the lncRNAs were shorter than 1500 nucleotides while only 10% were longer than 2000 nucleotides ( Figure 3C ). With an average of 1.88 exons per lncRNA, a major proportion (84.85%) of the lncRNAs in all six tissues were derived from 1 or 2 exons while only 1.68% were from ≥5 exons ( Figure 3D ). Chromosomal distribution of the identified lncRNA was depicted in ( Figure 3E ) and further visualized using the Circos software ( Figure 4 ).

To check on the conservation level of identified black pepper lncRNAs, BlastN against the previously known and reported lncRNAs of 38 plant species was performed with 10^-20 as e-value cutoff. Only 45 high-confidence matches were found where 42 identified black pepper lncRNAs matched with 27 database lncRNAs from 13 plant species with a maximum of five matches to Oryza barthii and four matches to Oryza rufipogon and Medicago truncatula each ( Supplementary Table 4 ). This result agrees with the literature suggesting the very poor conservation levels of lncRNAs compared to protein-coding mRNAs and are species- and tissue-specific.

lncRNAs are long RNAs present in the nucleus, nucleolus, and or cytoplasm and can function as the precursor for smaller ncRNAs such as snRNAs, snoRNAs, as well as the miRNAs (18–23 nts). To find out the black pepper lncRNAs which can possibly be precursors of known miRNAs, a BlastN search against the mirBase database was done ( Supplementary Table 5 ). A total of 14 pre-miRNAs were matched >=90% with 36 unique lncRNAs, which suggests those lncRNAs can give rise to the mature miRNAs after being acted upon by nuclease enzymes like dicer and/or drosha. Figure 5A shows lncRNA TCONS_00294108 containing the precursor and mature sequences of miRNA vvi-miR156i.

Micro RNAs (miRNAs) are small ncRNAs (18–23 nts) which have the major function of mRNA expression regulation by binding the 3′ UTR of the protein coding mRNAs. The expression is suppressed or silenced depending upon the complementarity of the miRNA-mRNA binding. A full complementarity leads to mRNA degradation thus silencing while a partial binding decreases the mRNA expression level downregulating the genes. LncRNAs sometimes interfere in the process and act as a miRNA sponge and prohibit miRNA-mRNA binding. Identified lncRNAs and known plant miRNAs available at the psRNAtarget server were taken for this analysis. Identified black pepper lncRNAs were submitted to the psRNAtarget server as target sequences against the available plant miRNAs and run with default parameters of max UPE 25 and expectation ≤2 ( Supplementary Table 6 ); 1702 interactions with 1054 unique miRNAs and 396 unique lncRNAs were found ( Figure 6 ). Figure 5B shows the RNAFold structure of lncRNA TCONS00240461 (blue) with the binding sites of miRNAs Osa-miR159e (red) and Zma-miE159c (green). Target mRNAs of the identified miRNAs were also found using psRNAtarget by submitting black pepper CDS (coding sequence) as target and identified 1054 miRNAs as small RNAs ( Supplementary Table 7 ); 3274 total interactions were found between the black pepper genes and miRNAs suggesting the possible gene regulations involved. Individual lncRNA-miRNA and miRNA-mRNA networks were merged and visualized using the Cytoscape software ( Figure 7A ), where we can see various networks involving miRNA and target mRNAs which can potentially be meddled by lncRNA thus regulating the normal gene regulation mechanism. For instance, miRNA pta-miR156b meant to be targeting mRNA Pn2.1339 is also capable of targeting four lncRNAs viz. TCONS_00067586, TCONS_00076072, TCONS_00072122, and TCONS_00076071. ( Figure 8A ).

After removing the low-quality reads from the raw sequences, the clean reads (>1.2 billion reads) were aligned to the reference genome using BWA-mem generating 53 SAM files for each sample which were then merged corresponding to each of the six tissues and submitted to a reliable, sensitive, and widely used command line tool CIRI2. A total of 4621 distinct circRNAs with ≥2 backspliced reads were identified including 3871 novel circRNA. A total of 750 circRNAs were found to be overlapping with known plant circRNAs in PlantcircBase. Tissue-wise analysis revealed 333, 2939, 1004, 156, 400, and 386 circRNAs identified in the flower, fruit, leaf, panicle, root and stem tissues, respectively ( Supplementary Table 8 ). Most the circRNAs were found to be ≤1000 nt (3077 circRNA smaller than 1000 nt) with a median length of approximately 400 nt, which is in consist with previous reports ( Figure 9A ) (Zheng et al., 2016). As described in previous studies, the abundance of circRNA was found to be lower than that of mRNA and lncRNA. Despite being widely distributed over all the chromosomes, genomic origin analysis showed most of the identified circRNA were coming from intergenic (47.98%) region followed by exonic (43.37%) regions and only few (8.66%) came from the intronic portion of the genome ( Figures 9B, C ).

Various studies suggesting circRNAs can regulate gene expression by acting as ceRNAs for the miRNAs by competing and inhibiting the miRNA binding with the mRNA molecules. The result of TargetFinder revealed 2591 interactions where 690 unique miRNAs were found targeting 432 circRNAs in all tissues combined ( Supplementary Table 9 ). CircRNA Pn4:30202222-30202678 was found to contain the putative miRNA binding site of aly-miR838-3p. mRNA targets of the identified miRNAs were identified using the psRNAtarget web server and a total of 6213 miRNA-mRNA interactions were found in which 619 unique miRNAs were targeting 744 unique black pepper mRNAs. Unique miRNA Pn3.1899 was found to be targeting 15 different mRNAs of black pepper ( Supplementary Table 10 ). The miRNA–circRNA–mRNA network relationship was visualized using Cytoscape software ( Figures 7B , 8B ).

In this study, we have focused on the trans-regulating action of lncRNAs and circRNAs where they regulate the target gene function by sequestering the miRNAs thus acting as ceRNAs in the cell. A total of 781 and 744 mRNAs were found to have an interaction in 6406 lncRNAs and 4621 circRNAs. To know the functional significance of the identified target mRNAs, GO annotation was performed. Results of the GO analysis revealed important GO subcategories involved, for instance “protein modification process”, “transmembrane transport”, “anatomical structure development”, and “signaling” annotated in the biological process category of GO. Similarly, “transferase activity”, “catalytic activity”, and “oxidoreductase activity” were annotated in the molecular function category and “nucleus”, “membrane”, and “plastid” in the cellular component category of GO. KEGG pathway analysis further revealed gene functions. Highly enriched pathways included the ras singaling pathway, diterpenoid biosynthesis, fatty acid biosynthesis, and plant–pathogen interaction. Both the GO and KEGG analyses strongly suggest the potential role of lncRNAs and circRNAs in diverse biological processes in the black pepper plant ( Supplementary Tables 11–14 ).

The black pepper non-coding RNA database, BPncRDB, is freely accessible at http://backlin.cabgrid.res.in/bpncrdb/index.php. The website features six tabs including Home, Search lncRNA, Interaction, circRNA, Download and Team ( Figure 10 ). The database catalogues the black pepper’s 6406 lncRNAs and 4621 circRNAs along with the miRNAs’ interactions involved with them. The database also contains the information regarding the mRNA/genes involved in the ceRNA pathways with the lncRNA/circRNAs. All of the analyzed data that is available in the MySQL database can be downloaded using the links on the “Download” page.