Hundreds of LncRNAs Display Circadian Rhythmicity in Zebrafish Larvae

Wild-type AB strain zebrafish were raised at the Soochow University Zebrafish Facility according to standard protocols [40]. Wild-type embryos were produced by pair mating and then raised in E3 (5 mM NaCl, 0.17 mM KCl, 0.33 mM CaCl2, and 0.33 mM MgSO4) embryo medium at 28.5 °C. To obtain larvae under constant dark (DD) or constant light (LL) conditions, embryos were first raised under the normal light/dark (LD, 14/10 h) condition for the first four days post fertilization (dpf) to activate and entrain the circadian system, and then transferred into the DD or LL environment. The samples were collected in darkness by employing a faint red flashlight under the DD condition. The larvae were anesthetized by ice. All the samples were collected within 2 min. Total RNAs were extracted from approximately 50 of the wild-type larvae every four hours from 120–140 hpf under the constant dark or constant light conditions using TRIzol reagent (Invitrogen, Carlsbad, CA, USA). All procedures were approved by the Soochow University Animal Care and Use Committee (#SUDA20211013A01) and were in accordance with regulations of the government of the People’s Republic of China.

Total RNAs extracted from six stages, 120 hpf (CT0), 124 hpf (CT4), 128 hpf (CT8), 132 hpf (CT12), 136 hpf (CT16), and 140 hpf (CT20) were examined on 1% agarose gels for integrity. Their concentrations were examined with a Nanodrop 2000 spectrophotometer (Thermo scientific, Waltham, MA, USA). A total amount of 3 µg RNA per sample was used in the sequencing library preparations, which were generated using NEBNext^®Ultra^TM RNA Library Prep Kit (Ipswich, MA, USA) following the manufacturer’s instructions. Construction and sequencing of complementary DNA (cDNA) libraries were conducted at Biomarker Technologies Corporation (Beijing, China). Clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit (Illumina, PE-401-3001) according to the manufacturer’s instructions. After clustering, the library preparations were sequenced on an Illumina Hiseq 2000 platform. We used Perl scripts for removing the adapters for clean reads, calculated the Q20, Q30, and GC content, and duplication data, and then generated the raw reads. All the analysis in this study is based on clean FPKM (fragments per kilo base per million mapped reads) data with high quality.

We generated time-course data under both DD and LL conditions from the transcriptome analysis [10] of wild-type zebrafish larvae. The two replicates were collected at six time points with a 4-h long interval. The six time point data under the constant darkness (DD) condition included WTDD120 hpf (CT0), WTDD124 hpf (CT4), WTDD128 hpf (CT8), WTDD132 hpf (CT12), WTDD136 hpf (CT16), and WTDD140 hpf (CT20), whereas the data under the constant light (LL) condition included WTLL120 hpf (CT0), WTLL124 hpf (CT4), WTLL128 hpf (CT8), WTLL132 hpf (CT12), WTLL136 hpf (CT16), and WTLL140 hpf (CT20). We compared the sequences from both replicates under the same condition, i.e., Sample 1 WT (DD) and Sample 2 WT (DD) as well as Sample 1 WT (LL) and Sample 2 WT (LL), to find the common transcripts that turned up in both samples for each condition. The sequences that matched with an E-value close to 0 were considered to be the same one in both samples in each of the DD and LL conditions. The data for each of the time points from CT0 to CT20 were determined by averaging the time-course-specific numeric expression levels from each sample. Together, we generated 3446 unannotated transcripts with six-point data under both LL and DD conditions. The expression profiles of 3446 transcripts are listed in the Supplementary Table S1 (DD) and Supplementary Table S2 (LL).

In order to identify the lncRNAs from the 3446 unannotated zebrafish larval transcripts, we compared these larval transcript sequences with the known zebrafish lncRNA sequences from the ZFLNC lncRNA database [29] and uncovered 3048 larval transcripts, out of the 3446 unannotated transcripts, which shared high similarities, with E-values lower than E-50, with ZFLNC lncRNAs. However, the remaining 398 larval transcripts could not be annotated with the ZFLNC database due to their higher E-values (>10⁻⁵⁰). We hypothesized that the lncRNAs catalogued in the ZFLNC database were far from complete. Hence, it is possible that the E-value-based sequence comparisons may not be able to identify some of the zebrafish larval lncRNAs. Several state-of-the-art computational tools [39] can be employed to identify the coding potentials of the unannotated transcripts. A previous study revealed that the Coding Potential Assessing Tool (CPAT) [41] has sensitivity and specificity of 0.96 and 0.97, respectively, to assess the coding potential of the novel transcripts. Therefore, we analyzed the coding abilities of the remaining 398 larval transcripts with the CPAT. Out of 398 unannotated transcripts, 172 transcripts were found to have coding abilities lower than 0.01. Hence, we designated these 172 transcripts as additional lncRNAs. With the ZFLNC database-annotated 3048 lncRNAs and CPAT-determined 172 lncRNAs, we obtained a total of 3220 lncRNAs from zebrafish larvae (Supplementary Table S3). The identifiers of the lncRNAs, such as Gene Bank IDs, and Ensemble IDs, were taken from the ZFLNC database as well as NCBI BLAST search results (Supplementary Table S3).

To compare zebrafish larval lncRNAs with those of two different zebrafish organs/tissues, we exploited the rhythmically expressed 586 and 165 lncRNAs from the zebrafish pineal gland and testis, respectively, from our recently published study [8].

We employed MetaCycle [42] to determine the rhythmicity of zebrafish larval lncRNAs with a statistical significance p-value of 0.05. Out of the 3220 zebrafish larval lncRNAs, the MetaCycle analysis determined 269 and 309 lncRNAs displaying circadian expression under DD and LL conditions, respectively (Supplementary Tables S4 and S5). For visualization of the expression patterns of these circadianly expressed lncRNAs, we further applied the BioDare2 system (https://biodare2.ed.ac.uk/↗, accessed on 1 October 2021) [43].

Numerous studies suggest the regulation of rhythmic expression of genes by cis regulatory motifs [8,44,45]. In particular, E-Box, D-Box and RORE motifs can regulate the expression profiles in the morning, evening, and night, respectively [46,47,48]. However, due to the lack of the promoter sequence of some lncRNAs, we investigated the promoter sequences of 115 (DD) and 112 (LL) circadianly expressed zebrafish larval ncRNAs (Supplementary Figure S1B) [49]. We BLASTed these larval lncRNAs against the NCBI’s Nucleotide collection (nt) database and selected the orthologs based on E-value-based sequence similarity and determined the Gene Bank Accession Number for each of the orthologs. The Gene Bank Accession Numbers were mapped to Ensemble IDs by bioDBnet (biological DataBase network) [50]. Supplementary Tables S6–S8 list Ensemble IDs of zebrafish larval lncRNAs under the DD condition, whereas Supplementary Tables S9–S11 show Ensemble IDs of zebrafish larval lncRNAs under the LL condition. Further, we downloaded a 5000-nucleotide 5′ upstream promoter sequence for each of the Ensemble IDs from Ensembl BioMart. We searched the E-Box, D-Box and RORE elements with a Matlab (The Mathworks, Inc., Natick MA, USA) program and the motif search tool Find Individual Motif Occurrences (FIMO) [51]. We downloaded the probability distribution for each of the E-Box, D-Box and RORE elements in FIMO from JASPAR [52] and kept the p-value threshold to 0.01 in the FIMO motif search. The promoter analysis helped establish a direct correlation between promoter elements and rhythmic expression of these zebrafish larval lncRNAs.

We conducted Gene Ontology (GO) [53] analysis, Clusters of Orthologous Groups (COGs) [54] analysis, and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis for all these circadianly expressed zebrafish larval lncRNAs. The GO analysis was performed with a Cytoscape open source Java application (https://cytoscape.org↗, accessed on 1 October 2021) [55] as well as with Blast2GO [56] software that employs the BLAST [57] algorithm to perform sequence alignment with the sequences in the nucleotide collection (nt) database in order to determine potential orthologs and the GO annotations. For all the lncRNAs orthologs with known Ensembl IDs, we also used BMKCloud (https://international.biocloud.net/↗, accessed on 1 October 2021) tools to visualize the GO annotation as well as to perform COGs analysis. A KEGG pathway enrichment analysis was performed with both BMKCloud and Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.8 (https://david.ncifcrf.gov/↗, accessed on 1 October 2021) [58].

We applied Principal Component Analysis (PCA) to compress the multidimensional time-course data into the two most important principal components [59]. The PCA assigned different scores derived from the gene expression profiles for each of the lncRNAs for all the principal components. We separately performed PCA for expression profiles of the morning, evening, and night lncRNAs. The absolute values of the PCA scores from the first principal component were used to rank the lncRNAs. The top-ranked representative lncRNAs from both the WT (DD) and WT (LL) datasets were plotted for visualization.

We employed both local and NCBI BLAST algorithms implemented in Blast2GO [56] to compare the DNA sequences. The execution of the BLAST algorithm provided expected value (E-value) output as the measure of statistical significance of the similarity for a sequence comparison. A lower E-value (≤10⁻⁵) represented a match of high statistical significance.

The amino acid sequences for the selected zebrafish larval lncRNAs and their mouse and human orthologs were translated into their corresponding peptide sequences with the ExPASy translate tool [60]. We separately applied the Clustal Omega (https://www.ebi.ac.uk/Tools/msa/clustalo/↗, accessed on 1 October 2021) tool to perform multiple sequence alignment for each of these selected zebrafish larval lncRNA-encoded peptides with the peptides encoded by testicular and pineal gland orthologs, as well as the mouse and human orthologs. The 3D models of conserved peptides were predicted with the (PS)2-v2: Protein Structure Prediction Server [61,62] and subsequently visualized with Jmol (http://www.jmol.org/↗, accessed on 1 October 2021).

The circadian clock genes display unique rhythmic expression patterns at different times over the 24-h cycle. Some genes peak during the morning, some during the evening, and others at night [46,47,48]. We investigated whether these zebrafish larval lncRNAs also exhibit peak expression during morning, evening, and night. Hence, rhythmicity analysis with MetaCycle revealed that approximately 8.35% (269 out of 3220) lncRNAs displayed circadian rhythmicity under the DD condition (Figure 1A–D; Supplementary Table S4). The correlative expression profiles of these 269 circadianly expressed larval lncRNAs are shown in a heat map (Figure 1A), while their phases are exhibited as BioDare2 plots (Figure 1B). We performed PCA of these circadianly expressed larval lncRNAs and ranked them based on the scores of the first principal component (Supplementary Figure S2A). The expression patterns and phases of the top-ranked representative lncRNAs are shown in Figure 1C,D. Subsequently, based on the expression patterns, we classified these 269 larval lncRNAs into 100 morning (CT0 and CT4) lncRNAs (E-H), 75 evening (CT8 and CT12) lncRNAs (I-L), and 94 night (CT16 CT20) lncRNAs under the DD condition (Supplementary Table S5), which were ranked according to their PCA scores (Supplementary Figure S2B,D). As visualized by expression patterns (Figure 1G,K,O) and the BioDare2 plots (Figure 1H,L,P) for representative lncRNAs, distinct expression patterns were observed for these morning lncRNAs, evening lncRNAs and night lncRNAs.

We hypothesized that the circadian expression of these 269 larval lncRNAs could be mediated via E-box, D-box or RORE promoter motifs. Therefore, we searched for the presence of these promoter motifs in the 5′ promoter sequences of each of these 100 morning lncRNAs, 75 evening lncRNAs, and 94 night lncRNAs. However, due to the lack of promoter sequences of some larval lncRNAs, we could only find the Ensembl IDs and the corresponding promoter sequences for 42 morning lncRNAs (), 33 evening lncRNAs (), and 40 night lncRNAs (). Subsequently, we searched for the presence of E-box, D-box or RORE elements in these promoter sequences with both the Matlab program and FIMO online application (). Initially, we searched the presence of the motif elements using the Matlab program (). We found that all 42 morning lncRNAs contain the E-box variable motif CANNTG, where N can be any nucleotide (). Further, seven out of thirty-three evening lncRNAs contain the D-box variable motif TTAYGTAA (), while seven out of forty night lncRNAs contain the RORE elements (). We then investigated all 114 lncRNAs’ promoter sequences with the FIMO utility. Interestingly, FIMO analysis revealed that all 42 morning lncRNAs’ promoters contain E-box motif (), all 33 evening lncRNAs’ promoters contain D-box motif (), and all 40 night lncRNAs’ promoters contain RORE elements (). Taken together, these analyses suggest circadian regulation of these circadianly expressed zebrafish larval lncRNAs. Although further experimental validation would be needed to confirm these findings, our computational analyses strongly suggest that under the DD condition, the rhythmicity of zebrafish larval morning lncRNAs are likely to be mediated via E-box, the rhythmicity of zebrafish larval night lncRNAs via D-box, and the rhythmicity of zebrafish larval night lncRNAs via the RORE elements. Supplementary Table S6 Supplementary Table S7 Supplementary Table S8 Supplementary Table S9 Supplementary Tables S9–S12 Supplementary Table S10 Supplementary Table S11 Supplementary Table S12 Supplementary Table S10 Supplementary Table S11 Supplementary Table S12

We then performed Gene Ontology (GO) analysis, COG functional annotation [54], and KEGG pathway enrichment analyses for all these 269 zebrafish larval lncRNAs displaying circadian rhythmicity under the DD condition (Figure 2A, Supplementary Figure S3A and Table S13). Out of the 269 lncRNAs, we were able to find the GO annotation, such as #GO, GO IDs, GO Names, and Enzyme codes, for 151 lncRNAs (Supplementary Table S13), and revealed their potential involvement in several fundamental biological activities, such as regulation of response to DNA damage stimulus, regulation of cellular response to stress, DNA replication, and regulation of DNA. The COG functional classification suggested numerous novel annotations, including amino acid transport and metabolism, translation and transcription (Figure 2B). Further, the KEGG pathway enrichment analysis uncovered gene numbers and corresponding q-values enriched in a diverse set of pathways, including Wnt signaling and sugar metabolism (Figure 2C, Supplementary Figure S3B and Table S14). Among the morning, evening and night lncRNAs with the known orthologs Ensemble IDs (Supplementary Tables S6–S8), we found 75 lncRNAs could be involved in different biological pathways (Supplementary Figure S3B). For instance, lncRNA Danio_rerio_newGene_6700 (FP245396.14↗) is potentially involved in Wnt signaling and melanogeneis, whereas lncRNAs Danio_rerio_newGene_21806 (ENSDARG00000032997) could be involved in protein processing in the endoplasmic reticulum. Taken together, we uncovered 269 zebrafish larval lncRNAs displaying circadian rhythmicity under the DD condition and determined their GO, COG and KEGG pathway enrichments.

Similar rhythmicity analysis revealed approximately 9.6% (309 out of 3220) larval lncRNAs showing circadian expression under the LL condition (Supplementary Table S15). The correlative expression profiles (Figure 3A) and expression phases (Figure 3B) of these 309 circadianly expressed larval lncRNAs under the LL condition are shown as heat map and BioDare2 plots. These 309 circadianly expressed larval lncRNAs were ranked based on their PCA scores (Supplementary Figure S4A), and the expression patterns and phases of the top-ranked representative larval lncRNAs are visualized in Figure 3C,D. We further categorized these circadianly expressed lncRNAs into 78 morning (CT0 and CT4) lncRNAs (Figure 3E–H), 144 evening (CT8 and CT12) lncRNAs (Figure 3I–L), and 87 night (CT16 and CT20) lncRNAs (Figure 3M–P) (Supplementary Table S16). For each of the morning, evening, and night lncRNAs, we separately performed PCA (Supplementary Figure S4B–D) and ranked the lncRNAs according to their PCA scores (Figure 3G,K,O). BioDare2 plots show the distinct expression patterns for these morning lncRNAs, evening lncRNAs, and night lncRNAs (Figure 3F,H,J,L,P).

Subsequently, the 5′ promoter sequences of all these circadianly expressed lncRNAs under the LL condition were investigated (). Out of the 78 morning lncRNAs, 144 evening lncRNAs, and 87 night lncRNAs, we were able to find the Ensembl IDs and the corresponding promoter sequences of 30 morning lncRNAs (), 58 evening lncRNAs (), and 24 night lncRNAs (). Together, we found 112 zebrafish larval lncRNAs under the LL condition with known Ensemble IDs. Similarly, we initially searched for the presence of the E-box, D-box or RORE promoter motifs in the 5′ promoter regions of these larval lncRNAs. The batch search analysis of motif sequences with the Matlab program revealed that all 30 morning lncRNAs’ promoter regions contain the E-box variable motifs (), whereas 14 evening lncRNAs’ promoter regions contain the D-box () and only one night lncRNAs’ promoter region contains RORE motifs (). We further interrogated these promoter regions with the FIMO utility. Interestingly, FIMO analysis showed that all 30 morning lncRNAs’ promoter regions contain the E-box motifs (), all 58 evening lncRNAs’ promoter regions contain D-box (), and all 24 night lncRNAs’ promoter regions contain RORE motifs (). Our findings, despite constrained by the availability of the promoter region sequences of lncRNAs, strongly suggest circadian regulation of these circadianly expressed zebrafish larval lncRNAs under LL condition. In particular, the zebrafish larval morning lncRNAs, evening lncRNAs, and night lncRNAs are likely to be mediated via the E-box, D-box or RORE promoter motifs, respectively. Supplementary Tables S17–S23 Supplementary Table S17 Supplementary Table S18 Supplementary Table S19 Supplementary Tables S20 and S21 Supplementary Table S22 Supplementary Table S23 Supplementary Tables S20 and S21 Supplementary Table S22 Supplementary Table S23

We also performed GO, COG, and KEGG pathway analyses of these zebrafish larval lncRNAs under the LL condition with known orthologs. Out of the 309 circadianly expressed lncRNAs under the LL condition, we were able to determine GO annotations for 173 lncRNAs with various genomic annotations, such as GO IDs, GO Names (Supplementary Table S24). In particular, GO analysis established the potential roles of circadianly expressed lncRNAs in a diverse set of fundamental biological processes, such as digestive system development, exocrine system development, RNA catabolic process, exocrine pancreas development, mRNA catabolic process, nuclear-transcribed mRNA catabolic process, response to inorganic substances, and response to metal ion protein-containing complex localization (Figure 4A, Supplementary Figure S5A and Table S24). The COG classification elucidated potential involvement of these zebrafish lncRNAs in important biological activities, such as cell cycle control, cell division, translation, transcription, and amino acid transport and metabolism (Figure 4B). KEGG pathway enrichment analysis determined statistics of pathway enrichment as well as direct role of lncRNAs in biological activities. For example, the KEGG statistics of pathway enrichment highlighted gene numbers with statistically significant q-values that regulated cellular activities and pathways (Figure 4C, Supplementary Figure S5B and Table S25). Further, 93 zebrafish larval lncRNAs under the LL condition are potentially involved in multiple biological pathways. For example, a zebrafish larval lncRNA Danio_rerio_newGene_5943 (ENSDARG00000104170) is potentially involved in nucleic acid binding and metal ion binding, the other lncRNA, Danio_rerio_newGene_17037 (ENSDARG00000093359), is likely to be involved in inositol phosphate metabolism, metabolic pathways, and the phosphatidylinositol signaling system (Supplementary Table S25). The percentages of genes potentially involved in biological processes, cellular components, and molecular functions are shown in Supplementary Figure S5B. Overall, we determined 309 circadianly expressed larval lncRNAs under the LL condition and their GO, COG, and KEGG pathway annotation analyses suggest these circadianly expressed lncRNAs likely contribute to various biological processes.

In order to understand the effect of light on gain or loss of rhythmicity we compared circadianly expressed lncRNAs under the DD condition with those under LL condition and found that a total of 239 lncRNAs that were circadianly expressed under the DD condition, but not the LL condition (Supplementary Table S26, Figure 5 and Supplementary Figure S6), implicating loss of rhythmicity, likely due to exposure to constant light. Interestingly, a total of 279 lncRNAs were circadianly expressed only under the LL condition (Supplementary Table S26 and Figures S7 and S8), implying gain of rhythmicity under constant light. Further, 30 lncRNAs were coexpressed under both DD and LL conditions (Supplementary Table S26). Heat map and BioDare2 plots show the correlative expression profiles (Figure 5A) and phases (Figure 5B) of these 30 coexpressing lncRNAs. The expression patterns and phases of the top-ranked lncRNAs based upon their PCA scores (Supplementary Figure S6A) are shown in Figure 5C,D. We further categorized these 30 lncRNAs into morning, evening, and night lncRNAs.

Based on the expression profiles under the DD condition, we found eleven, eleven and eight lncRNAs exhibited peak expression in the morning (Figure 5E–H), evening (Figure 5I–L), and night (Figure 5M–P), respectively (Supplementary Table S27). Further, we separately performed PCA of the expression profiles of these morning lncRNAs, evening lncRNAs, and night lncRNAs under the DD condition (Supplementary Figure S6B–D). BioDare2 plots show the distinct expression patterns for the morning lncRNAs, evening lncRNAs, and night lncRNAs (Figure 5F,H,J,H,N,P). In contrast, based on the expression profiles of 30 lncRNAs under LL conditions, there were nine, fourteen and seven lncRNAs that peaked in the morning (Supplementary Figure S7E–H), evening (Supplementary Figure S7I–L), and night (Supplementary Figure S7M–P), respectively (Supplementary Table S28). PCA of the expression profiles of these morning lncRNAs, evening lncRNAs, and night lncRNAs under the LL condition allowed for selecting top-ranked representatives (Supplementary Figure S8A–D). BioDare2 plots show the distinct expression patterns for the morning lncRNAs, evening lncRNAs, and night lncRNAs (Supplementary Figure S7F,H,J,N,P). Similarly, promoter analysis suggested that rhythmicity of these 30 circadianly expressed zebrafish larval lncRNAs under both DD (Supplementary Tables S9–S12) and LL (Supplementary Tables S20–S23) conditions are likely mediated via the E-box, D-box or RORE motifs in the morning, evening, and night, respectively.

We also compared the changes in BioDare2 phase values of the 30 lncRNAs under DD and LL conditions (). The results suggest that 12 lncRNAs gained phase due to exposure to light as their phase values increased from the DD to the LL condition, while 17 lncRNAs exhibited decreased phase values under the LL condition. Interestingly, one lncRNA Danio_rerio_newGene_3412 (ENSDARG00000029881) showed exactly the same phase value under both the DD and LL condition because of very similar expression values under both DD and LL conditions (). Supplementary Table S28 Supplementary Table S28

All 30 coexpressed lncRNAs were investigated to determine their GO, COG, and KEGG pathway annotations. GO analysis suggested the potential roles of these coexpressed circadian lncRNAs in diverse sets of biological processes, such as cellular response to metal ions, cellular response to calcium ions, and cellular response to inorganic substances (Figure 6A, Supplementary Figure S9A and Table S13). The COG classification elucidated their involvement in important biological activities, such as posttranslational modification, protein turnover, and chaperones (Figure 6B). The KEGG pathway analysis determined the statistics of pathway enrichment as well as the possible roles of lncRNAs in biological activities (Figure 6C, Supplementary Figure S9B and Table S25). Numerous circadianly expressed zebrafish larval lncRNAs are potentially involved in multiple biological pathways. For example, a zebrafish larval lncRNA Danio_rerio_newGene_11483 (ENSDARG00000098252) is likely to be involved in calcium-dependent membrane targeting, whereas the other lncRNA Danio_rerio_newGene_6700 (ENSDARG00000063437) is likely to be involved in Wnt signaling and melanogenesis. Taken together, we identified 30 zebrafish larval lncRNAs coexpressed under both DD and LL conditions and predicted their possible functional roles.

Our recent study [8] identified 586 and 165 rhythmically expressed lncRNAs in zebrafish pineal gland and zebrafish testis, respectively, and demonstrated that 26 rhythmically expressed lncRNAs are coexpressed in both the zebrafish pineal gland and testis [8]. We hypothesized that a certain number of circadianly expressed larval lncRNAs may be rhythmically/circadianly expressed in zebrafish pineal gland and testis, even though larvae and adult organs exhibit distinct developmental programs [63,64]. Therefore, we compared these circadianly expressed zebrafish larval lncRNAs under DD and LL conditions with those of the pineal gland testis lncRNAs (Figure 7, Supplementary Tables S29–S34) and revealed that under the DD condition, 60 circadianly expressed larval lncRNAs were coexpressed in the pineal gland (Supplementary Table S29), while 39 circadianly expressed larval lncRNAs were coexpressed in the testis (Supplementary Table S30); and nine circadianly expressed larval lncRNAs were coexpressed in both the pineal gland and testis (Figure 7A, Supplementary Table S31). In contrast, under the LL condition, 64 circadianly expressed larval lncRNAs were coexpressed in the pineal gland (Supplementary Table S32), while 58 circadianly expressed larval lncRNAs were coexpressed in the testis (Supplementary Table S33); and 12 circadianly expressed larval lncRNAs under the LL condition were coexpressed in both the pineal gland and testis (Figure 7B, Supplementary Table S34).

We further analyzed the expression profiles of the circadianly expressed zebrafish larval lncRNAs from both DD and LL conditions that were coexpressed in the pineal gland and testis. Under the DD condition, we found that out of the nine lncRNAs coexpressed in the larvae, pineal gland, and testis, three lncRNAs, two lncRNAs and four lncRNAs peaked in the morning, evening and night, respectively (Figure 7A,C–E, Supplementary Table S31), while under the LL condition, out of the 12 coexpressed lncRNAs, we identified five lncRNAs, four lncRNAs and three lncRNAs peaking in the morning, evening and night, respectively (Figure 7B,F–H, Supplementary Table S34). Overall, we found that zebrafish larvae coexpress nine and twelve lncRNAs in both pineal gland and testis under DD and LL conditions, respectively, and we analyzed their expression profiles in the morning, evening and night.

We also investigated that how many coexpressing zebrafish larval lncRNAs encode peptides. In particular, we computationally predicted the 3D structures and functions of these lncRNAs coexpressing under both DD and LL conditions. Under the DD condition, 267 peptides were computed to be encoded by these nine lncRNAs coexpressed in zebrafish larvae, pineal gland, and testis (). In contrast, under LL condition, 226 peptides were predicted to be encoded by these 12 lncRNAs coexpressed in zebrafish larvae, pineal gland, and testis, respectively (). Supplementary Table S35 Supplementary Table S36

We further computationally predicted the 3D structures and functions of these coexpressing lncRNA-encoded peptides (Figure 8, Supplementary Tables S37 and S38). Out of the 267 peptides encoded by the nine coexpressing larval lncRNAs under the DD condition, we were able to predict 3D structures for three lncRNA-encoded peptides including 9_DANIO_RERIO_NEWGENE_1127_PEPTIDES_19 (CU682604.7↗), 6_DANIO_RERIO_NEWGENE_8560_PEPTIDES_3 (ENSDARG00000105693), and 8_DANIO_RERIO_NEWGENE_1627_PEPTIDES_2 (ENSDARG00000069440) (Figure 8A–C and Supplementary Table S37). In contrast, out of the 226 peptides encoded by the 12 coexpressing larval lncRNAs under the LL condition, we were able to predict the 3D structure for one lncRNA-encoded peptide 1_DANIO_RERIO_NEWGENE_3710_PEPTIDES_10 (XM_002667273.5↗) (Figure 8D). The 3D structures revealed the presence of α-helix, β-strand, and random coils (Figure 8A–D) among the coexpressing lncRNA-encoded peptides. Subsequently, we analyzed the 3D structures by investigating the known domains from Protein Data Bank (http://www.rcsb.org/↗, accessed on 1 October 2021) to determine the conserved protein domains for each of the 3D model (Supplementary Tables S37 and S38).

Interestingly, all the predicted 3D structures mapped to the known domains in the Protein Data Bank. For example, under DD condition, the peptides 9_DANIO_RERIO_NEWGENE_1127_PEPTIDES_19(CU682604.7↗), and 6_DANIO_RERIO_NEWGENE_8560_PEPTIDES_3 (ENSDARG00000105693) were mapped to the known domains 1vg5A and 1gxrA, respectively (Figure 8). As another example, under the condition, the peptide 1_DANIO_RERIO_NEWGENE_3710_PEPTIDES_10 (XM_002667273.5↗) was mapped to the known domains 1kpfA. Together, we predicted 3D structures for the four lncRNA-encoded peptides and determined their domains from the Protein Data Bank.

We investigated the conservation of circadianly expressed zebrafish larval lncRNAs with humans, mice, and fruit flies. In particular, we determined if some of these circadianly expressed zebrafish larval lncRNAs under DD and LL conditions, which are also coexpressed in the pineal gland and testis lncRNAs, are conserved with the aforementioned three species. In order to find the conservation of these lncRNAs, we determined human, mouse and fruit fly orthologs of these larval lncRNAs from both DD and LL conditions (Figure 9A,B, Supplementary Tables S39–S41). With the E-value-based similarity of sequences, we found that under the DD condition, 35 circadianly expressed zebrafish larval lncRNAs are conserved with humans (Supplementary Table S39), while only one circadianly expressed zebrafish larval lncRNA is conserved in mice (Supplementary Table S40), and only one circadianly expressed zebrafish larval lncRNA is conserved simultaneously in both humans and mice (Supplementary Table S41). Specifically, the larval lncRNA Danio_rerio_newGene_17148 (XM_039656533↗) has a human ortholog (AV748161.1↗) and a mouse ortholog (CF745161.1↗). We further divided the 35 larval lncRNAs conserved with humans into the morning, evening and night groups, and visualized the corresponding representative lncRNAs (Supplementary Table S39, Figure 9C–E). There were 12 lncRNAs, 9 lncRNAs, and 14 lncRNAs peaking in the morning, evening and night, respectively. Interestingly, out of the nine lncRNAs coexpressed in zebrafish larvae, pineal gland and testis, two lncRNAs, DANIO_RERIO_NEWGENE_18433 (ENSDARG00000012446) and DANIO_RERIO_NEWGENE_19201 (ENSDARG00000105693), are part of the 35 lncRNAs conserved with humans. However, the one larval lncRNA under DD condition conserved with mice does not belong to the nine coexpressing lncRNAs.

In contrast, under the LL condition, 42 circadianly expressed zebrafish larval lncRNAs are conserved with humans (Supplementary Table S42), while four circadianly expressed zebrafish larval lncRNAs are conserved with mice (Supplementary Table S43); and four circadian zebrafish larval lncRNAs are conserved simultaneously with both humans and mice (Supplementary Table S44). We also categorized 42 larval lncRNAs (LL) conserved with humans into the morning, evening and night groups and visualized the corresponding representative lncRNAs (Supplementary Table S42, Figure 9F–H). There were 8 lncRNAs, 26 lncRNAs, and 8 lncRNAs peaking in the morning, evening and night, respectively. Out of the 42 lncRNAs conserved with humans, three lncRNAs, DANIO_RERIO_NEWGENE_11679 (ENSDARG00000060849), DANIO_RERIO_NEWGENE_15482 (ENSDARG00000101175), and DANIO_RERIO_NEWGENE_2130 (ENSDARG00000098304), belong to the 12 lncRNAs coexpressed in the zebrafish larvae, pineal gland, and testis. However, none of the four lncRNAs conserved with mice belongs to the 12 coexpressing lncRNAs. Further, none of the circadianly expressed larval lncRNAs from either the DD or LL conditions is conserved with any in fruit flies. Taken together, our conservative analysis revealed that a higher number of zebrafish larval lncRNAs under both DD and LL conditions are more conserved with humans than with mice.

Subsequently, we analyzed the conserved lncRNA-encoded peptides from both DD and LL conditions and investigated their functions from the Protein Data Bank. Specifically, under DD condition, we identified 19, 8, and 4 peptides encoded by these conserved zebrafish larval lncRNAs, their human orthologs, and mouse orthologs, respectively (). The conservative analysis of the lncRNA-encoded peptides revealed one peptide conserved among zebrafish larvae, mice, and humans (). In contrast, under LL condition, we identified 132, 12, and 25 peptides encoded by the conserved zebrafish larval lncRNAs, their human orthologs, and mouse orthologs, respectively (). However, the conservative analysis of the lncRNA-encoded peptides revealed that no peptides are conserved among zebrafish larvae, mice, and humans. Together, we uncovered one lncRNA-encoded peptide conserved among these three species. Supplementary Tables S44–S50 Supplementary Table S47 Supplementary Tables S48–S50

We then performed multiple sequence alignments of the peptides encoded by this lncRNA and its orthologs to investigate its 3D structures and functions. However, there were no matching domain and 3D structures of the conserved peptides encoded by the conserved lncRNA. Overall, we uncovered one conserved peptide encoded by the lncRNA conserved among zebrafish larvae, mice, and humans.

The data supporting the reported results can be found in the. supplementary materials