Altered non-coding RNA expression profile in F1 progeny 1 year after parental irradiation is linked to adverse effects in zebrafish

Gamma radiation produces DNA instability and impaired phenotype. Previously, we observed negative effects on phenotype, DNA methylation, and gene expression profiles, in offspring of zebrafish exposed to gamma radiation during gametogenesis. We hypothesize that previously observed effects are accompanied with changes in the expression profile of non-coding RNAs, inherited by next generations. Non-coding RNA expression profile was analysed in F1 offspring (5.5 h post-fertilization) by high-throughput sequencing 1 year after parental irradiation (8.7 mGy/h, 5.2 Gy total dose). Using our previous F1-γ genome-wide gene expression data (GSE98539), hundreds of mRNAs were predicted as targets of differentially expressed (DE) miRNAs, involved in pathways such as insulin receptor, NFkB and PTEN signalling, linking to apoptosis and cancer. snRNAs belonging to the five major spliceosomal snRNAs were down-regulated in the F1-γ group, Indicating transcriptional and post-transcriptional alterations. In addition, DEpiRNA clusters were associated to 9 transposable elements (TEs) (LTR, LINE, and TIR) (p = 0.0024), probable as a response to the activation of these TEs. Moreover, the expression of the lincRNAs malat-1, and several others was altered in the offspring F1, in concordance with previously observed phenotypical alterations. In conclusion, our results demonstrate diverse gamma radiation-induced alterations in the ncRNA profiles of F1 offspring observable 1 year after parental irradiation.

www.nature.com/scientificreports/ expression profiles, as well as their contribution to the establishment of gene expression patterns, disorders, and phenotypes under the influence of radiation, remain unknown. Gamma radiation-induced alterations at transcriptional, and DNA level can be inherited by the offspring of vertebrates. We recently demonstrated the inheritance of altered mRNA expression profiles, DNA methylation, and histone modifications patterns in zebrafish embryos after parental gamma irradiation [8][9][10] . The DNA methylation analysis along with the mRNA expression profile of F 1 embryos revealed pathways associated with gamma radiation response; such as molecular mechanisms of cancer, DNA damage response and cell death, along with pathways not previously seem involved in the response to gamma radiation, like signalling and retinoic acid receptor activation, and gonadotropin-releasing hormone (Gnrh) signalling. On the other hand, F 1 embryos showed enriched methylation at histone marks such as H3K4me3, H3K4me9, and H3K27me3, indicating alterations in chromatin structure and organization, as well as in the expression of developmental genes like hepatocyte nuclear factor 4 alpha (hnf4a) [8][9][10] .
Nevertheless, the inheritance of dysregulated profiles for most of the sncRNA classes in zebrafish offspring, as well as their involvement in the mechanisms underlying the response to gamma radiation, as a result of parental exposure, remains unclear.
Moreover, our previously published work on mRNA expression profile after parental gamma irradiation, using the same sample materials as for the current study, showed genes such as dicer, ago1, ago2, ago3b, ago4, and piwil2 to be down-regulated in F 1 offspring of zebrafish 9 . The products of these genes participate in the miRNA, siRNA, and piRNA biogenesis pathways, suggesting a potential impact of gamma radiation on their biogenesis and expression in the progeny of gamma-exposed parents.
In the present study, we aimed to investigate the effect of gamma radiation on the sncRNA expression profile in F 1 embryos 1 year after parental exposure through small RNA sequencing and subsequent gene expression analysis.

Results and discussion
Sequencing analysis. The sncRNA transcriptome in zebrafish embryos has previously been characterized and consist of several classes of small RNAs such as tRNA-halves, tRNA fragments, piRNAs, and miRNAs among others [11][12][13][14] . In this study, we have focused on sncRNA expression through small RNA-seq to analyse the sncRNA profile in 5.5 hpf F 1 offspring embryos 1 year after parental exposure to gamma radiation (8.7 mGy/h) during gametogenesis.
In average, 11,970,959 and 11,389,554 reads resulted from sequencing of F 1 -γ (8.7 mGy/h) and control samples (F 1 -C), respectively. All sequencing libraries had a Phred score higher than 20; however, sequencing reads were filtered by quality using a Phred score of 30. Approximately 95% of reads were recovered after quality filtering (Phred score > 30) over all libraries (data not shown). The read length distribution showed three main peaks present in all libraries, which corresponds to the theoretical sizes of miRNAs (22 nt), piRNAs (26-31 nt), and tRNA-halves (32-34 nt) (Fig. 1). Mapped reads, showed significant differences between offspring of gammaexposed parents (F 1 -γ) compared to controls (F 1 -C) (80.4% and 88.5%, respectively, p < 0.0001) (Fig. 2). An average of 47.6% of the mapped reads over all libraries aligned to multiple loci (> 5) on the zebrafish reference genome (Fig. 2). It is well known that reads aligning to features such as miRNAs, piRNAs, and tRNA-derived fragments align to multiple genomic regions, including repetitive elements 11,15,16 .
Mapped reads in both groups (F 1 -γ and F 1 -C) were extensively annotated to different genomic features (Fig. 3A,B). Reads mapping to piRNA were significantly enriched in the offspring of gamma-exposed parents (F 1 -γ) (1.4-fold, p < 0.001) as compared to the controls (F 1 -C). Whereas those mapping to lincRNA and snRNA were significantly depleted in the offspring from exposed parents (F 1 -γ) (2.1-and 4.1-folds, respectively p < 0.001). We found no differences with regard to the number of reads mapped to miRNAs or any other genomic feature. Differential miRNA expression in F 1 embryos from exposed parents. Normalization of miRNA expression datasets has been discussed in recent studies, and several methods have given controversial results in miRNA normalization [17][18][19] . In our case, TMM proved to be suitable to normalize the miRNA expression data, Figure 1. Reads length distribution of raw sequence reads after filtering by read length. Three biological replicate libraries per group of F 1 generation (5.5 hpf embryos) from parents exposed to 8.7 mGy/h of γ-radiation (F 1 -γ) and control parents (F 1 -C). www.nature.com/scientificreports/ and the exploration after normalization showed a clear effect of gamma radiation on miRNA expression in the offspring F 1 of exposed parents (F 1 -γ), compared to controls (F 1 -C) (Supplementary Figure S1). The differential expression analysis between offspring F 1 of gamma-exposed parents and controls, using a cut-off of > 30 read counts in all replicates, showed 22 DEmiRNAs (log2 FC > 0.6, p < 0.05, FDR < 0.05), from those, 55% were up-regulated ( Fig. 4A; Table 1).
Several of the DEmiRNAs listed in Table 1, such as miR-21, let-7g, and miR-150 have previously been reported as affected by ionizing radiation in directly exposed organisms including zebrafish embryos 4,20 . miR-21 is known to act as an oncomiR involved in cancer-related processes 21 , and overexpression of let-7g increase radiosensitivity in lung cancer 4 . In other experimental studies, miR-193b-3p, miR-23b, and members of let-7 family were affected following irradiation 4,21,22 . Besides, DEmiRNAs dre-miR-200a-5p and dre-miR-141-3p, members of family miR-8, are implicated in the control of pluripotency, cancer proliferation, and metastasis 23 . Interestingly, these two DEmiRNAs are oppositely modulated in the offspring of the gamma-exposed parents (Table 1), indicating different regulatory functions. In developmental studies in mouse and zebrafish, this family has been found coexpressed in epithelial and olfactory cells, respectively 24 . In summary, a large fraction of the DEmiRNAs found in offspring F 1 from gamma-expose parents is related to DEmiRNAs previously found in directly exposed models, indicating that alterations in miRNA expression might be inherited intergenerationally.
Ingenuity pathway analysis and miRNA target filter. To further explore the possible pathways under miRNA control we used our previously published mRNA dataset (GSE98539) 9 , generated from the same batch of embryos as used in this work, and used the IPA miRNA target filter to link differentially expressed genes (DEGs) to DEmiRNAs (Table 1). This resulted in 12 DEmiRNAs that were imported in IPA of which 11 were in the IPA knowledge base. We found 672 DEGs linked to the DEmiRNAs of which 380 showed an inverse relationship in expression rate with their counteracting DEmiRNA. The pathway analysis on these targets showed overrepresented pathways such as insulin receptor, NFkB, and PTEN signaling ( Table 2, Supplementary Data S1). let-7 and miR-21 showed the largest involvement in these specific pathways. Interestingly, IPA's molecule activity predictor revealed that the involved miRNAs followed the inverse expression relationship as specified in the miRNA target filter (Supplementary Data S1) resulting in predicted effects on apoptosis, transcription and inflammation (Supplementary Figures S2-S4). In directly exposed organisms, including zebrafish embryos, miR-125b acts as a negative regulator of P53, affecting apoptosis 25,26 . In this study, we did not find any modulation of miR-125b in the offspring from gamma-exposed parents (F 1 -γ). However, apoptosis was predicted as an outcome of the interaction of let-7g with a specific network of genes, which included dicer, ago1, ago2, and ago3 involved in the miRNA biosynthesis pathway (Fig. 5). Furthermore, most overrepresented disease functions were linked to cancer, gastrointestinal and hepatic disease and developmental disorders (Table 2) supporting the observed genomic instability found in larvae from exposed parents by our group 27 and indicating a vast involvement of miRNAs in DNA damage related response.

Figure 2.
General statistics of reads mapping to zebrafish genome (GRCz10). Three biological replicate libraries per group of F 1 generation (5.5 hpf embryos) from parents exposed to 8.7 mGy/h of γ-radiation (F 1 -γ) and control parents (F 1 -C). Total clean reads represent the number of reads after QC analysis and filtering. Filtered reads indicate the number of reads after size filtering. Total mapped reads indicate the number of reads mapped to zebrafish reference genome after size filtering. Multimapping reads < 5 include reads aligning to < 5 genomic locations. Multimapping reads > 5, indicate reads mapped to > 5 genomic locations. Unmapped reads show the amount of reads which could not be aligned to the genome. Asterisks represent significant differences p < 0.0001, Chi-square with Yates correction of continuity. www.nature.com/scientificreports/ Differential expression of piRNA clusters in F 1 embryos from exposed parents. In order to determine the expression of piRNAs, candidate piRNAs from both experimental groups (F 1 -γ and F 1 -C) were assigned into clusters. In total, 171 piRNA clusters were predicted from both groups overlapping in genomic position (matched clusters). In concordance with previous studies 11, 28 , matched clusters were found expressed from the sense strand (25.7%), antisense strand (48.5%), and bidirectional orientation (25.7%), ranging from 4.3 to 64.0 kb (Supplementary Data S2). After dataset normalization (Supplementary Figure S5), 11 piRNAs clusters were differentially expressed between offspring of gamma-exposed parents and controls, 6 up-regulated, and 5 down-regulated (log2 FC > 0.6, p < 0.05, Adjusted p value < 0.05; Fig. 4B, Table 3).
piRNA population within predicted clusters was also searched for piRNA signatures such as uridine (U) at 5′ end, which is recognized as a hallmark of the primary biogenesis pathway, and 10 nt 5′ overlap bias, accepted as a signature for the secondary piRNAs biogenesis pathway, known as ping-pong mechanism 29,30 . Besides the read length (Fig. 6A), found in concordance with the described theoretical size of piRNAs, more than 80% of Figure 3. Distribution of mapped reads (zebrafish genome reference GRCz10) onto genomic features. F 1 generation of embryos (5.5 hpf) from parents 1 year after exposure to 8.7 mGy/h γ-radiation (F 1 -γ) and control parents (F 1 -C). Values derived from three replicates in each group (n = 3). Asterisks denote statistical differences with p < 0.001, multiple t tests (FDR 99%). (A) miRNA total (total of reads mapping to microRNAs (miRNA)), miRNA danio_rerio (zebrafish miRNAs), miRNA other (miRNAs from other species), Mt_tRNA (mitochondrial transference RNA), piRNA (piwi-interacting RNA), no-annotation (mapped reads, which could not be assigned to any genomic feature), protein_coding (protein coding transcripts), rRNA (ribosomal RNA), tRNA (genomic transference RNA), others (sum of all reads mapping to genomic features presented in (B). www.nature.com/scientificreports/ piRNA reads in both experimental groups contained a U at their 5′end (Fig. 6B). In addition, 9.87% and 10.25% of piRNA candidate reads in F 1 offspring from exposed and controls parents respectively, showed the typical 10 nt 5′overlap (Fig. 6C). The significantly larger number of reads in the offspring from the gamma-exposed group (F 1 -γ), showing a 10 nt 5′ overlap compared to controls (p < 0.05), suggested a higher abundance of piRNAs in the group F 1 -γ derived from processed transposons through ping-pong mechanism. In zebrafish, Piwil1 process the activated transposons producing secondary piRNAs, which in turn are loaded into Piwil2, generating more piRNAs from primary piRNA transcripts. Thus driving the amplification loop (ping-pong mechanism) resulting in a population a piRNAs of opposite polarity with a 10 nt 5′overlap 31 . The analysis of ping-pong signatures showed the presence of a similar heterotypic 28 × 29nt peak in both F 1 -γ and control group (F 1 -C). Whereas a strong peak at 26 × 28nt was only found in the offspring form exposed parents (F 1 -γ) (average Z-score F 1 -C = 44.2, and F 1 -γ = 42.9; p < 0.01), indicating that piRNA population in the descendants from the gamma-irradiated parents (F 1 -γ) differs in size from that found in the offspring from the non-irradiated parents (F 1 -C) (Fig. 6D). Since there are only two piwi paralogs reported in zebrafish, piwil1 (piwi-like RNA-mediated gene silencing 1, ENSDARG00000041699) and piwil2 (piwi-like RNA-mediated gene silencing 2, ENSDARG00000062601), the same peak from 10 nt 5′ overlapping reads might be observed in both experimental groups. Besides, considering that, piwil2 was found down-regulated when analysing mRNA expression in a parallel study using siblings embryos as for those used in this study 9 , we then hypothesized that parental gamma radiation induced the activation of additional piwil2 paralog or isoform, which might be the cause for the presence of a different peak from 10 nt 5′ overlapping reads in the group F 1 -γ. There are four different piwil2 isoforms annotated in Ensembl (ENSDART00000090695.7, ENSDART00000162071.2, ENSDART00000134274.3, ENSDART00000136004.2). Nevertheless further research is necessary to prove this hypothesis.
Differentially expressed piRNA clusters are associated to transposable elements. Although piRNAs appear to be involved in transcriptional regulation and deadenylation of mRNAs 32 , its major role in the gonads is in controlling the expression of TEs, thus protecting the genome from their harmful effects, and in F 1 generation of embryos (5.5 hpf) from parents exposed to 8.7 mGy/h γ-radiation (F 1 -γ) and non-exposed control parents (F 1 -C). Differentially expressed genes are represented as red dots. Expression values are shown as log2 of fold changes (X-axis). Y-axis represents the negative log10 of the p values (n = 3). www.nature.com/scientificreports/ ensuring that genetic information is correctly passed down to the next generation 33,34 . Primary piRNA clusters are derived from genomic loci known for harboring TEs. Transcriptionally active TEs are processed by Piwil1 and Piwil2 producing secondary piRNAs (ping-pong mechanism), which are antisense to expressed TEs 31,33 . We intersected our set of piRNA clusters expressed in the offspring from the gamma-exposed group (F 1 -γ) and the controls (F 1 -C) (matched piRNA clusters) with TEs. As expected, we found a large number of matched piRNA clusters (85 out of 171) overlapping with 172 TEs (p < 0.0001), which were classified into six different orders; LTR, LINE, TIR, Crypton, Helintron and DNA/unknown (Supplementary Data S3). Afterward, we determined a significant association between DEpiRNA clusters and TEs (p = 0.0024), where 45% of DEpiRNA clusters overlapped with 9 TEs belonging to orders LTR, LINE, and TIR. In every case, DEpiRNA clusters appeared to  Table 3. Differentially expressed piRNA clusters. Pairwise comparison between generation F 1 from exposed (8.7 mGy/h γ-radiation) (F 1 -γ) and control parents (F 1 -C). www.nature.com/scientificreports/ be expressed from the complementary strand of the associated TEs ( Table 4), suggesting that these DEpiRNAs have been expressed in response to the activation of TEs.
In bovine embryos, piRNAs and TEs of the families LINE and ERV1 among others, that seemed to be associated considering their mapping locations, showed an inverse relationship in expression 35 . Additionally, TEs activation due to irradiation has been observed in yeast, plants, and Drosophila [36][37][38] . In concordance, TEs have been associated with DNA integrity disruption and double-stranded breaks in Drosophila 39 . Therefore, parallel experiments on the siblings of these embryos showed a significantly higher DNA damage rate in the offspring of the gamma exposed parents as compared to embryos from control (non-irradiated) parents 1 year after parental irradiation, which lacked of association with an effect from reactive oxygen species 27 .
Though piRNAs-mediated TEs silencing through recruitment of DNA methyltransferase has been proposed as a mechanism for controlling TEs activation in mammals 40 , we found no association between piRNA clusters Figure 6. Analysis of piRNA signatures in clusters predicted in F 1 generation of embryos (5.5 hpf) from parents exposed to 8.7 mGy/h γ-radiation (F 1 -γ) and non-exposed control parents (F 1 -C) (n = 3). A) Read length distribution. B) Per base distribution analysis. C) 5′10nt overlap distribution. Asterisk indicates significant differences (p < 0.05; paired t test). D) Read length analysis of reads with 5′ 10nt overlapping (ping-pong signatures; average Z score F 1 -C = 44.2, and F 1 -γ = 42.9; p < 0.01). Table 4. Genomic association between DEpiRNA clusters and transposable elements (TE) in F 1 offspring from parents exposed to gamma radiation. Fisher's exact test, p-value = 0.0024.

Class
Order Superfamily www.nature.com/scientificreports/ and DMRs obtained from our parallel study 8 . To our best knowledge, this is the first study reporting an association between altered expression of piRNA clusters and TEs in the generation F 1 of parents exposed to gamma radiation.
Differential expression of snRNA in F 1 embryos from exposed parents. We observed a significantly lower amount of reads mapping to snRNAs in F 1 offspring from gamma-irradiated parents (F 1 -γ) as compared to the F 1 offspring of non-irradiated parents (F 1 -C) (4.1-folds, p < 0.001) (Fig. 3B). Using a cutoff of more than 100 read counts, the expression of 25 snRNA genes was detected (Supplementary Data S4). Pairwise comparison on the normalized dataset (Supplementary Figure S6) showed 19 (76%) snRNA genes differentially expressed (DEsnRNA) (log2FC > 0.6, p < 0.05, FDR < 0.05). From those, 16 were down-regulated in the group F 1 -γ, whereas only three snRNAs were up-regulated ( Fig. 4D; Table 5). Small nuclear RNAs comprise a group of nuclear-localized sncRNAs, which are critical components of the spliceosome. From the five major spliceosomal snRNA (U1, U2, U4, U5, and U6), the snRNA U1 and U2 guide the association of the spliceosome to the 5′ and 3′splice sites respectively, and are essential for the early splicing process 41 . Interestingly, the vast majority of down-regulated snRNA genes (12 out of 16) in the group F 1 -γ were members of U1 and U2. In addition to U1 and U2, four members of U4 and U5 were found down-regulated. In total, 4 out of 5 major spliceosomal snRNAs were represented in our data as down-regulated genes ( Table 5).
We hypothesize that the down-regulation of the major spliceosomal snRNAs negatively affects the gene expression rate. Both snRNAs U1 and U2 participate in the regulation of transcription by increasing the formation of the first phosphodiester bond during transcription initiation, and also interact with components of the initiation complex such as TFIIH (transcription factor II human) [41][42][43] . On the other hand, gtf2h1 (general transcription factor 2 h polypeptide 1), a gene encoding for a subunit of TfIIh, was found down-regulated when analysing mRNA expression of siblings embryos as for those used in this study 9 . The later suggests that combined down-regulation of the major spliceosomal snRNAs and gtf2h1 can be the cause of the lower global gene expression rate in F 1 offspring of gamma irradiated parents from our parallel gene expression study 9 .
Differential expression of lincRNA in F 1 embryos from exposed parents. We found a significantly lower amount of reads mapping to lincRNA genes in the embryos from parents exposed to gamma radiation (F 1 -γ) than in controls (F 1 -C) (2.1-folds, p < 0.001) (Fig. 3B).
Aiming to compare the expression of lincRNAs between F 1 -γ and control (F 1 -C) embryos, the reads mapping to lincRNAs in both groups were quantitated. We detected the expression of 108 lincRNA genes with a cut-off of more than 100 read counts (Supplementary Data S5). After normalization of the quantitated dataset (Supplementary Figure S7), 44 lincRNAs were differentially expressed (DElincRNA) in the group F 1 -γ (Fig. 4C). We shortlisted the number of DElincRNAs to 21 genes using a more stringent fold change cut-off (log2FC > 0.6, p < 0.05, FDR < 0.05) that could be categorized into 14 (66.7%) down-regulated and 7 (33.3%) up-regulated lincRNAs in the offspring of the gamma irradiated parents (F 1 -γ) ( Table 6). Table 5. Differentially expressed snRNAs. Pairwise comparison between generations F 1 from exposed (8.7 mGy/h γ-radiation) (F 1 -γ) and control parents (F 1 -C). Transcript IDs were obtained from Ensembl database (http://www.ensem bl.org). www.nature.com/scientificreports/ To gain insight into the biological function of the DElincRNAs, we searched for DElincRNAs within our dataset with conservation and functional annotation in the ZFLNCRNA database 44 . Only five out of 21 DElin-cRNAs were found conserved in human or mouse representing the 23.8% (Table 7). Similarly, previous studies found only 5.3 to 8.9% of zebrafish lncRNAs conserved 45,46 .
Conserved genes are thought to exert similar biological functions across species. One of the most downregulated DElincRNAs found in the descendants from gamma-irradiated parents represents an ortholog of the Table 6. Differentially expressed lincRNAs. Pairwise comparison between generations F 1 from exposed (8.7 mGy/h γ-radiation) (F 1 -γ) and control parents (F 1 -C). Transcript IDs and transcript names obtained from Ensembl database (http://www.ensem bl.org) and ZFIN (The Zebrafish Information Network; http://zfin.org/) respectively.  www.nature.com/scientificreports/ well-characterised human lincRNA MALAT1 (metastasis associated lung adenocarcinoma transcript 1) (Tables 6,  7). MALAT1 has been shown to interact with several serine/arginine (SR) proteins, driving the distribution of splicing factors in the nucleus. The down-regulation of MALAT1 in HeLa cells resulted in decreased association of splicing factors to the nuclear speckle, producing alterations in the alternative splicing of endogenous pre-mRNA 47 . Besides, cellular responses to radiation exposure, DNA damage, and apoptosis can affect the alternative splicing process 48 . Therefore, the strong down-regulation of malat-1, and the down-regulation of major spliceosomal snRNAs suggest alterations in the alternative splicing in the generation F 1 from exposed parents. In addition, some lincRNAs can target miRNAs and counteract their function by reducing specific miRNA availability 49 . We used the experimental module from DIANA-LncBase (v.2) 50 to search for experimentally confirmed miRNA-lncRNA interactions. Human orthologs of miRNAs dre-miR-21-1-3p, dre-let-7g-5p, dre-miR-1-3p, and dre-miR-135c-5p, among others (Table 1), were found as targets of MALAT1 (Supplementary Figure S8). This set represent up-regulated miRNAs in our dataset in contrast to the down-regulation of malat-1. Our pathway analysis indicated the involvement of those DEmiRNAs in the regulation of pathways such as insulin receptor-signaling, NFKB signaling, and pTEN (Supplementary Figures S2, S3, S4).
On the other hand, in zebrafish embryos, the knockdown of malat-1 caused a high mortality rate and various types of phenotypical alterations 51 . This is in concordance with the expression profile of malat-1 in our data, suggesting that its down-regulation contributed to establishing altered phenotypes documented in the F 1 offspring from gamma-irradiated parents 9,27 . To the best of our knowledge the rest of conserved DElincRNAs represent uncharacterised members of this class, without available biological information.
We retrieved the co-expression networks of those DElincRNAs in our data with functional annotation in the ZFLNCRNA database (Table 7). This resulted in a subset of 59 genes representing unique entries. When intersecting this subset of genes with our previously published parallel mRNA-Seq data 9 , 32 genes (54.2%) had quantitated expression levels, whereas 27 genes (45.8%) were not found (Supplementary Data S6). The 32 genes detected were classified based on their expression levels as modulated (FDR < 0.05) or unmodulated (FDR > 0.05). The vast majority of modulated genes showed expression levels opposite to the DElincRNAs, while unmodulated genes had significantly lower expression levels as compared to modulated genes (p < 0.05) (Fig. 7A). This indicates that down-regulation of DElincRNAs result in increased expression of genes within the co-expression networks. The GO analysis of modulated genes under the category of biological processes showed significant enrichment in five GO terms (FDR < 0.05), such as ribosome biogenesis (GO:0042254), ribonucleoprotein complex biogenesis (GO:0022613), rRNA metabolic process (GO:0016072), rRNA processing (GO:0006364), and ribosomal large subunit biogenesis (GO:0042273) (Fig. 7B). Besides, the pathway ribosome biogenesis in eukaryotes (ID 03009) was found overrepresented (p value 2.406E−05, FDR 4.019E−03). The GO analysis of unmodulated genes did not show any GO terms significantly overrepresented.
The synthesis of ribosomes must be under strict control to guarantee correct cell growth and proliferation. Disturbances in this pathway result in altered cell cycle and proliferation, similar to that observed in cancer 52 . In addition, the ribosome biogenesis pathway is sensitive to radiation in directly exposed mice 53 . www.nature.com/scientificreports/ Based on publicly available data we have analysed the function of a small subset of DElincRNAs. We lack information about to the biological function of many other DElincRNAs in our data, limits our capacity to unravel all the implications of the dysregulation of this ncRNA class, in the descendants of parents exposed to gamma radiation. Despite not being fully understood, our results show the potential involvement of DElincR-NAs in the regulation of miRNAs responsive to gamma radiation, post-transcriptional processes, and protein synthesis.

Conclusions
Our results show the effect of gamma radiation on the non-coding transcriptome in the first-generation offspring of exposed parents. The exposure of adult zebrafish during gametogenesis to a dose rate similar to that observed in Chernobyl 60 days post-accident, produced alterations on the expression profile of different classes of ncR-NAs such as miRNAs, piRNAs, snRNAs, and lincRNAs, which were observable 1 year later in the generation F 1 . Moreover, pathway analysis, and ncRNA expression could be linked to previously observed gamma radiation effects on phenotype and gene expression of these embryos. Altogether, our results provide new knowledge on the involvement of ncRNAs in the response to gamma radiation, and contribute to better understand to what extent the adverse effects of gamma radiation are inherited by offspring. Furthermore, this work sets the bases for transgenerational studies, with a focus on mechanisms of inheritance, expression pattern maintenance, as well as the possible impact on population dynamics.

Material and methods
Zebrafish husbandry and exposures. Zebrafish strain AB wild type was obtained from the Norwegian University of Life Sciences (NMBU) zebrafish facility and kept according to standard operational procedures as previously described by Hurem et al. 7 . The exposures of fish to gamma radiation, mating, as well as the generation of embryos were carried out as previously described by Hurem et al. 27 . Briefly, adult zebrafish (6 months old) were exposed for 27 days to a 60 Co source at 8.7 mGy/h (5.2 Gy total dose), which represents a similar dose rate to that observed in Chernobyl 60 days post-accident 54 . Control fish were kept separately under similar experimental conditions. In both exposure and control groups, three replicates of 30 males and 30 females were used.
This study was conducted under the approval of the Institutional Animal Ethics Committee (IACUC) and the Norwegian Food Inspection Authority (NFIA), under permit number 5793. NMBU zebrafish facility, is licensed by the NFIA and accredited by the association for assessment and accreditation of laboratory animal care (AAALAC, license number: 2014/225976). The NMBU zebrafish facility and SOPs has AAALAC accreditation (No. 1036) and is approved by the National Animal Research Authority. All experiments were performed according to Norwegian Animal Welfare Act (2009) and the EU Directive 2010/63, following appropriate guidelines.
Embryo sampling. One year after exposure, fish were mated, and embryos were obtained as previously described by Hurem et al. 27 . First-generation embryos from both, exposed and control parents were harvested in 100 mm Petri dishes (Costar, Corning incorporated, USA) containing autoclaved system water. The embryos were incubated at a controlled temperature of 28 ± 2 °C until staging at 5.5 hpf (50% epiboly). The autoclaved system water was replaced every two hours. The staging of embryos was performed following Kimmel et al. 55 . Staged embryos from each experimental group were directly transferred into 12-well plate (Costar, Corning incorporated, USA), containing 3 mL of autoclaved system water at a controlled temperature of 28 ± 2 °C. Pools of 100 embryos from each replicate in both experimental groups were allocated in separate wells, and washed 3 times with 10 mL autoclaved system water (28 ± 2 °C) to remove debris from the previous incubation period. Afterward, each pool of embryos was transferred to 1.5 mL microfuge tubes (Thermo Fisher Scientific, Waltham, MA) and snap-frozen in liquid nitrogen. All samples were stored at − 80 °C for further analysis.

RNA isolation.
Batches of 100 embryos in triplicates were homogenized using Magnalyser Beads (Roche Diagnostics, Germany), and total RNA was isolated using Trizol following the instructions of the manufacturer (Thermo Fisher Scientific, USA). Exogenous synthetic kanamycin mRNA (Promega, USA) was added as spikein RNA to the Trizol (0.25 ng/mL). RNA integrity and concentration were assessed with Bioanalyzer using the RNA Nano LabChip Kit (Agilent Technologies, Santa Clara, CA) and Nanodrop 1000 (Thermo Fisher Scientific, USA) respectively. RNA samples with RIN values above eight were stored at − 80 °C until used for sequencing.
Small RNA sequencing. The total RNA samples were sent for custom sequencing (Novogen, Hong Kong) under Illumina platform (HiSeq 4000). Three single-end libraries (biological replicates) from the exposure group (F 1 -γ), and the control group (F 1 -C) were made following the manufacture's recommendations. In short, sequencing libraries were prepared from each sample using 1 μg of total RNA as input material (NEBnext Small RNA Library Prep Set for Illumina, New England Biolabs inc, USA). Adaptors at 5′and 3′-ends were ligated to generate non-directional cDNA libraries. After reverse transcription, cDNA was amplified using 12 PCR cycles. Amplified cDNA libraries were purified using AMPpure XP beads kit and size-fractionated using 6% polyacrylamide gel to obtain the fraction corresponding to sncRNAs (up to 150 bp). Libraries were checked for quality using Bioanalyzer 2100 system (Agilent Technologies, Palo Alto, USA) with DNA 1000 chip kit (Agilent Technologies). The libraries, with read lengths of 50 nt, were sequenced with targeted depth 10 million reads per library. Raw reads were deposited at Gene Expression Omnibus Database under accession number GSE152189.
Bioinformatics analysis. Raw reads were trimmed from adapter sequences and quality assessed using Trim Galore! v0.3.7 56 www.nature.com/scientificreports/ and QC filtering, raw reads were filtered by length using a combination of on terminal awk, unix, and perl commands. Only reads within the range 18-36 nt were retained. Length filtered reads were then mapped to zebrafish reference genome GRCz10 (http://www.ensem bl.org), classified into genomic features, and counted using Unitas (v1.5.2) 30 under default settings. We used miRDeep2 v0.0.5 58 to quantify miRNAs using the built-in script quantifier.pl, the option -W was used to weight read counts by their number of mapped loci. Options -p, -m, -r and -t were used to indicate miRNA precursor reference, mature miRNA reference, read files in FASTA format, and species, respectively. The rest of the options were kept as default. Precursor and mature miRNA references were obtained from miRBase v21 (http://www.mirba se.org) 59 , and updated with other reported miRNAs 60 .
Reads mapping to piRNA producing loci (piRNA cluster database released 01.08.2018, http://www.small RNAgr oup-mainz .de) and mapped reads that were not assigned to any known genomic feature were subsequently used for piRNA analysis. We first determined ping-pong signatures using Unitas (v1.5.2) 30 under the option -pp, to determine the 5′ 10 bp overlaps among all potentially piRNA mapped reads.
Once ping-pong signatures were determined, we proceeded to group piRNA reads by clusters to predict piRNA producing loci using proTRAC (v4.2.4) 34 . As pre-processes, we used the perl script TBr2-duster.pl (v2.1) 61 under default values to collapse mapped reads and remove low complexity sequences. Afterward, reads were remapped using sRNAmapper (v1.0.5) 61 to generate the input files for downstream processes. All options were left to default settings except for -alignment, which was set to best. Due to the presence of multimapping reads, the script reallocate.pl (v1.1, http://www.small RNAgr oup-mainz .de) was used to reallocate multimapping reads based on the calculation of estimated expression rates of uniquely mapping reads. The output files obtained from the reallocation step were used as input for proTRAC (v4.2.4) 34 under default values. A zebrafish repeat masker annotation for GRCz10 was obtained from http://www.repea tmask er.org and passed onto proTRAC.
Bedtools package (v2.29.1) 62 with the subcommand intersect was used to intersect the genomic coordinates of predicted piRNA clusters from the three replicates of each experimental group. The option -f was set to 0.50 to report any overlap of piRNA clusters between replicates within each group by more than 50%. Any piRNA cluster not found in all replicates was removed and not considered for further analysis.
Aiming to associate piRNA clusters to transposable elements (TEs), Bedtools package (v2.29.1) 62 under the same settings mentioned above, was used to intersecting the genomic coordinate of differentially expressed piRNA clusters (DEpiRNAs) in both experimental groups, and the previously downloaded repeat annotation reference. We also used previously published data from our group in an attempt to link piRNA clusters to previously published DNA methylation changes observed in the same batch of embryos used in this study 8 . Therefore, we used the genomic coordinates of the DEpiRNAs and overlapped those with observed differentially methylated regions (DMRs) (GSE100470) using Seqmonk 63 .
Differential expression analysis and statistics. Statistical differences in global read counts for reads mapped to genomic features between experimental groups were obtained by multiple t-tests followed by False Discovery Rate (FDR) calculation (Benjamini, Krieger, and Yekutieli method, Q = 1%.).
Genomic features such as miRNA, piRNA clusters, snRNA, and lincRNA (long intergenic non-coding RNA), were analysed for differential expression. All differential expression analyses were performed following the same analytical sequence as follows; data exploration, normalization, data exploration post-normalization, and pairwise comparison between experimental groups (F 1 -γ and F 1 -C; n = 3). Pre-normalization data log2 transformed was explored for descriptive statistics such as minimum, first quartile, median, third quartile, and maximum also the similarity among samples was determined by correlation and hclust (Ward method) analysis to determine the distance between samples. Multidimensional Scaling Plot (MDS-plot) was used to analyze the variances, except for piRNA clusters where Principal Component Analysis (PCA) was utilized. miRNA, lincRNA, and snRNA expression datasets were TMM normalized (trimmed mean of M-values, edgeR v3.24.3, Bioconductor) 64 , whereas piRNA clusters dataset was RPKM normalized to correct for differences in clusters length. Post-normalization data exploration was conducted similarly to pre-normalization, in both cases, the statistical package R v3.0.2 65 was used.
Normalized datasets were used for differential expression analysis. In the case of miRNAs, lincRNAs, and snR-NAs, the statistical analysis was based on the pairwise comparison between treatment (F 1 -γ) and control (F 1 -C) (n = 3). Experimental groups were compared using the exact test under Bioconductor package edgeR v3.24.3 64 . As for piRNA clusters, F 1 -γ and F 1 -C groups were compared through a pairwise comparison using the function voomWithQualityWeights followed by a linear model and eBayes from Bioconductor package limma (v3.36.3) 66 .
In all analyses, the FDR was set up to 95%. Only ncRNAs with significant FDR (p > 0.05) and log2 FC > 0.6 were considered as differentially expressed.
Fischer's exact test was used to search for overrepresentation by comparing the number of piRNA clusters overlapping all methylated regions and piRNA clusters merely overlapping differentially methylated regions (DMRs), as well as TEs overlapping with all piRNA clusters and those overlapping to DEpiRNA clusters. Kruskal-Wallis test followed by Dunn's multiple comparison (alpha 0.05, confidence level 95%) was used to compare the expression of lincRNAs, to genes classified as modulated, and unmodulated within the lincRNA co-expression networks.
Functional analyses. In order to use the miRNA target filter in Ingenuity Pathway Analysis (IPA, Qiagen, USA), the DEmiRNAs were manually converted to their human orthologue via miRBase 59 . Only miRNAs with a 100% match in seed sequence and a maximum of 2 mismatches in the mature sequence were used in the downstream analysis (Table 1). Within IPA, the miRNA list was linked to our previously published differentially expressed mRNAs (DEGs) (GSE98539) 9 . The miRNA target filter is based on computational as well as experi-  44 under default parameters (e-value 0.001, word size 11, sensitivity normal) to search for conservation and functional annotations. DElincRNAs not annotated as conserved in the database were then BLAST against the nucleotide database in NCBI (https ://www.ncbi.nlm.nhi.gov) to search for similar sequences. We used human and mouse as the target organisms. The megablast algorithm was selected and all parameters were left to default values (Supplementary Data S8).
We retrieved the co-expression network for all the DElincRNAs with functional annotation in ZFLNCRNA (gene ontology and KEGG pathway) and classified the genes as modulated or unmodulated by comparing the obtained gene list to our previously reported differentially expressed mRNAs (GSE98539) 9 . Expressed genes were considered as modulated with FDR < 0.05, whereas genes with FDR > 0.05 were registered as unmodulated (Supplementary Data S7). Later, gene ontology and pathway analyses on the generated gene sets were performed using WebGestalt (http://www.webge stalt .org) 67 based on over-representation analysis. Only GO terms within the biological process category and pathways from KEGG were considered. DIANA-lncBase (v2) 50 was used to search for experimentally confirmed miRNA-lncRNA interactions using the experimental module.
RT-qPCR validation. RT-qPCR was performed to validate the expression level of 9 (random selected) out of 22 annotated DEmiRNAs. For cDNA synthesis, 1 μg of total RNA from each sample was used as input material (miScript II RT kit, Qiagen, USA). Following manufacturer's instruction, 4 μL of 5 × miScript HiSpec Buffer, 2 μL of 10 × nucleotide triphosphate mix, 2 μL miScript reverse transcriptase mix, and 10 μL of RNAse free water, were mixed to reach a total volume of 20 μL. Reactions were incubated first for 60 min at 37 °C followed by 5 min at 95 °C to inactivate the reverse transcriptase. qPCR (miScript PCR system, Qiagen, USA) reactions were set up in a total volume of 25 μL containing 12.5 μL 2 × QuantiTec SYBR Green PCR Master Mix, 2.5 μL 10 × miScript Universal Primer, 2.5 μL miRNA forward primer (5 mmol/L), 2.5 μL template cDNA (diluted ten times in RNAse free water) and 5 μL of RNAse free water. The polymerase was activated during 15 min at 95 °C followed by a 3-step amplification program as follows: denaturation 15 s 94 °C, annealing 30 s 57 °C, extension 30 s 70 °C. Finally, 95 °C 10 s, 65 °C 60 s, and 97 °C 1 s were used to obtain the dissociation curves. The ramp rate was adjusted to 1 °C/s. Measurements were performed in a LightCycler96 (Roche Diagnostics, Switzerland) from three biological replicates per experimental group (F 1 -γ and F 1 -C). Linreg (v2017.0) 68 was used to determine primers efficiency and to calculate the expression level of miRNAs as measured by the threshold cycle values (Ct). Quantification values derived from five technical replicate per biological replicate. Relative quantification was normalized to a synthetic spike-in kanamycin RNA and expression levels were determined as equation 2 −ΔΔCT . Forward primers were designed with CLC-Workbench v6.0 (Qiagen, USA) ( Table 8).
The obtained mean relative miRNA expression values (F 1 -γ vs F 1 -C) were compared to mean relative miRNA expression values for the same miRNAs from RNA-seq, and a Spearman's correlation coefficient was calculated (p < 0.05) (Graphpad Prism 7 v7.0, La Jolla, USA).