Identification, characterization and target gene prediction of miRNAs
Most of the plant miRNAs are evolutionarily conserved from species to species [22, 23] and this indicates the powerful strategy for the identification of new miRNAs by using the already known miRNAs . Many conserved miRNAs have been identified from the expressed sequence tag (EST) [25, 26] and genome survey sequence (GSS)  by using this homology search approach. For mango, there is no GSS data and the available EST data for mango is only 1709 and it was not sufficient for identification of miRNA. Hence, unigene sequences (107,744) were used for the identification of miRNAs in this study. Unigene is a unique transcript that is transcribed from a genome and many miRNAs have been identified from the unigenes of many plant species such as Artemisia annua , coconut , Litchi fruit  and black pepper .
The potential 104 pre-miRNAs of mango were predicted based on the parameters of Zhang and the MFEI values were also calculated as the MFEI gave the best prediction of miRNAs . Although the length of the predicted mature miRNAs was in the range of 18-22 nt, the length of precursor miRNAs varied significantly from 67 to144 nt with an average length of 94 nt. The predicted 104 mango miRNAs belong to 86 different families. Among them, over 70% of the miRNA families have only one family member. The highest 5 family members were found in the miR2673 family followed by miR159, with 4 family members. The remaining miRNAs have the family member of 2 or 3. Therefore, we can see that the mango miRNA distribution across various families is highly heterogeneous.
The previous studies have already proved that the plant miRNAs bind to their targets in a perfect or nearly perfect complementarity and thus the psRNATarget server was used to search the target gene of mango miRNAs in this study. Both mRNAs collected from NCBI and mRNA identified in this study were used as the target candidates of miRNAs due to the absence of Mangifera indica target candidates in psRNATarget server. Some previous studies indicated that miR156 was a master regulator of the juvenile phase in plants and it targeted the squamosa promoter binding protein-Like (SPL) gene family to regulate the transition from vegetative phase to floral phase in Arabidopsis, maize and rice [33-39]. In mango, MmiR105772, a family of miR156, also bound to its target SPL6 and thus the predicted targets of mango miRNAs were in the agreement with the previously published papers in other plants. The resulting data from psRNATarget also showed that only one miRNA (MmiR1653) had the single target gene which was the member of miR482 family and bound to monodehydroascorbate reductase 4 enzyme, the important gene related to the nutritional quality of mango fruit . All other miRNAs could target to multiple genes and some miRNAs had over two hundred target genes. For example, MmiR73030 had 230 target genes and these target genes involved in 16 KEGG pathways such as biosynthesis of antibiotics, purine metabolism, sulfur metabolism, glycerophospholipid metabolism, T cell receptor signalling pathway, steroid degradation and so on.
Sivankalyani, Sela et al. published that the mango stress-response pathways were activated by cyclic nucleotide-gated channel (CNGC) and leucine-rich repeat receptor (Lrr) . In this study, we found that MmiR90392 targeted to CNGC1, and MmiR68471 and MmiR68478 targeted to Lrr2. Moreover, MmiR10167 and MmiR15558 bound to the stress WRKY transcription factor 44 which play a major role in plant defence to biotic and abiotic stresses. MmiR78769 and MmiR101928 also bound to their target genes of phospholipase A and phospholipase D which were key factors in plant responses to biotic and abiotic stresses . The ethylene response could improve the tolerance of mango fruit to chilling stress  and ten mango miRNAs identified in this study had six ethylene responsive target genes such as ethylene-insensitive protein and ethylene-responsive transcription factor. So, these newly identified mango miRNAs have potential roles in chilling stress responsive process of mango.
Two mango miRNAs (MmiR23777 and MmiR36814) also targeted to the auxin efflux carrier which had the potential role in mango plant organ development . A total of 17 miRNAs interacted with auxin-related genes. MmiR51876 was a miRNA that targeted to auxin responsive protein. The pentose and glucoronate interconversions pathway, phenlypropanoid biosynthesis pathway and alpha-linolenic acid metabolism pathway were KEGG pathway involved in the adventitious root formation of mango cotyledon segments . In this study, 9 miRNAs bound to 8 target genes involved in these three pathways for mango root formation. MmiR10167 bound to target genes that involved in phenlypropanoid biosynthesis pathway and MmiR7519 bound to target genes involved in alpha-linolenic acid metabolism pathway. From these findings, it could be observed that these five mango miRNAs (MmiR23777, MmiR36814, MmiR51876, MmiR10167 and MmiR7519) involved in the developmental process of mango (Figure 2B).
Identification, characterization and target gene prediction of lncRNAs
As the genome sequence of mango is not available till now, the de novo assembled transcriptome sequences were used for the identification of lncRNAs in this study. A total of 277,071 RNA transcripts from Zill, Shelly and Keitt mango cultivars studied by the former researchers were used and a total of 31,226 candidate lncRNAs were predicted in this study. Among them, 24 lncRNAs were significantly expressed to heat stress and 7586 lncRNAs to cold stress. The most significantly expressed down-regulated heat responsive lncRNA was HRlnc25944 with the fold change value of -6.22. HRlnc11351 and HRlnc27371 were the mostly expressed up-regulated lncRNAs with fold change value greater than 7. For the cold responsive lncRNAs, CRlnc10871 was the mostly expressed down-regulated lncRNA (FC value -11.19), and CRlnc26299, CRlnc30496 and CRlnc36473 were the most significantly expressed lncRNAs with fold change value greater than 11.
No heat responsive lncRNAs was conserved but 0.29% of cold responsive lncRNAs were conserved with 12 different plant species. Among them, the highest conserved lncRNAs were CRlnc32663 and CRlnc47883. Each of which was conserved with 4 different lncRNAs of other plants. CRlnc32663 conserved with 4 different lncRNAs of 3 three different plant species such as Manihot esculenta, Malus domestica and Populus trichocarpa. CRlnc42883 also conserved with 4 lncRNAs of Oryza rufipogon, Oryza barthii and Solanum lycopersicum (Figure 3B).
For heat responsive lncRNAs, 8 bound to 115 target genes involved in the plant development and stress response. HRlnc11351 was the most significantly expressed up-regulated lncRNAs with fold change value of 7.55 and bound to six heat shock proteins. In cold responsive lncRNAs, CRlnc26299 was one of the most significantly expressed up-regulated lncRNAs and bound to RC12B (JK513200_1) which is the low temperature and salt responsive protein found in Arabidopsis thaliana . The WRKY proteins are a large family of transcriptional regulators in higher plant and are exhibited the variable expression patterns in response to chilling stress in cucumber, mango and rice [40, 46, 47]. In this study, 64 cold responsive lncRNAs have interaction with WRKY gene family. So, we can observe that the cold responsive lncRNAs of mango have the interaction with the target genes that are expressed at the low temperature stress.
GO enrichment analysis and KEGG pathway analysis were performed for the better understanding of the target genes of newly identified lncRNAs. From GO enrichment analysis result, we could see that both types of heat responsive lncRNAs and cold responsive lncRNAs had interaction with the target genes involved in all three broad categories such as biological process, cellular component and molecular function. For the biological processes, both heat responsive and cold responsive lncRNAs were highly enriched in metabolic processes and cellular processes. In the cellular component analysis, GOs related to membranes, intracellular and cytoplasm were highly enriched for both type of lncRNAs. Also for molecular function analysis, most of the enriched GO terms in both type of lncRNAs were related to catalytic activity and binding. Therefore, we could see that, the GO terms highly enriched in both heat responsive and cold responsive lncRNAs were not quite different.
Among 17 KEGG pathways of the target genes of the heat responsive lncRNAs, amino sugar and nucleotide sugar metabolism was the most significant pathway and 8 target genes involved in this pathway. As mentioned above, HRlnc11351 was the most significantly expressed up-regulated lncRNAs and its target gene, JK513625_1 is 3-ketoacyl-CoA thiolase 2 (KAT2, EC:126.96.36.199) which could be mapped to 9 different pathways such as benzoate degradation, fatty acid elongation, biosynthesis of unsaturated fatty acids, alpha-linolenic acid metabolism, fatty acid degradation, valine, leucine and isoleucine degradation, biosynthesis of antibiotics, geraniol degradation and ethylbenzene degradation according to the result of KEGG pathway analysis. In Arabidopsis, KAT2 is an enzyme that catalyses the β-oxidation of fatty acid and involves in abscisic acid (ABA) signal transduction . The phytohormone ABA plays an important role in plant development and adaptation to diverse environmental stresses. Therefore, HRlnc11351 may involve in the important role of mango development and stress response by targeting to KAT2. For cold responsive lncRNAs, 209 target genes had mapped to 86 KEGG pathways. JK513026_1, alcohol dehydrogenase 1 (ADH1, EC:188.8.131.52) was the most enriched target gene and involved in 12 different pathways including glycolysis/gluconeogenesis, metabolism of xenobiotics by cytochrome P450, glycine, serine and threonine metabolism, methane metabolism, fatty acid degradation and so on. In plants, ADH genes are involved in mediating stress responses and developments. In mango, ADH1 has important role in the fruit ripening  and thus, cold responsive lncRNAs that target to ADH1 gene may play important role in mango fruit ripening process. According to the KEGG pathway analysis results, purine metabolism and biosynthesis of antibiotics were the highly enriched pathways among 86 pathways and more than 50 target genes were enriched in each pathway.
Interaction between lncRNAs and miRNAs
The interaction between the miRNAs and lncRNAs showed that the most of the miRNAs had targeted to more than one lncRNAs and only 8 miRNAs had single target lncRNAs. The number of lncRNAs targeted by a single miRNA was in the range of 1 to 90. A total of 90 target lncRNAs were found for MmiR73030 which also targeted 230 mRNAs. This miRNA had the highest target numbers in both lncRNAs and mRNAs.
LncRNAs not only can target to miRNAs to reduce the stability of lncRNAs but also can function as molecular decoys or sponges of miRNAs [[10, 50]]. So the miRNA target mimicry search was performed by using TAPIR, which is a web server for the prediction of plant miRNA targets including target mimics. Although no heat responsive lncRNA act as the target mimic of miRNAs, 20 cold responsive lncRNAs were predicted as the target mimics of 20 miRNAs. CRlnc31221 was the target mimic of MmiR5408 which targeted to 8 cold responsive lncRNAs and 47 target genes. These target genes were involved in starch and sucrose metabolism, inositol phosphate metabolism and phenylpropanoid biosynthesis pathways which pathways are important for plant growth and development, and plant response towards biotic and abiotic stress. During target mimicry, the interactions between miRNAs and their authentic targets were blocked by binding of decoy RNA to miRNAs via partially complementary sequences . So, the target mimicry of CRlnc31221 had the potential regulation effect to the interaction between the target genes and MmiR5408.