Transcriptome Analysis on the ALAN-Induced Zebrafish Ovary: Desynchronization of Life Processes and Initiation of Diseases
Article information
Abstract
Objective
The effect of artificial light at night (ALAN) on “transcriptome” is prominent owing to its capacity for “desynchronization” of organismal physiology. Light influences the circadian rhythm. This study aims to explore the ALAN-induced ovarian transcriptome of zebrafish for desynchronization of life processes.
Methods
Four experimental conditions were set up for female zebrafish: one normal 12-hour light and 12-hour dark (LD) cycle, and three continuous exposures to ALAN for one week (LLW), one month (LLM), and one year (LLY). The whole transcriptome data analysis of the ALAN-exposed samples was then compared with the normal sample using RNA-Seq, followed by exploratory analyses.
Results
The analysis revealed two different patterns of expression of genes where LLW and LLM differ with LLY samples in comparison to LD. Compared to LD, downregulation of the predicted hub genes was observed in all treatments; ribosome and oxidative phosphorylation pathways were enriched. LLY vs. LD contrast depicts the enrichment of three more pathways—RNA polymerase, adhesion junction, and signaling. The gene ontology (GO) enrichment portrays more prevalent biological processes in LLW vs. LD and LLM vs. LD than in LLY vs. LD. Contrast-wise disease annotation represents neoplasms as the most prevalent; disease enrichment denotes the major class of neoplasm, carcinoma, coupled with intellectual disability, global developmental delay, and seizures.
Conclusion
Our study displayed desynchronization of various genes and pathways leading to the initiation of diseases, for the first time in zebrafish. Even though, this data shows that ALAN is a serious threat, further research is needed to determine the intensity and the duration of ALAN which might cause potential repercussions.
INTRODUCTION
Long recognized as a man-made disruption, artificial light at night (ALAN) has become more prevalent, depriving Earth’s surface of its typical light cycles. The repercussions of change in the light-dark cycle on the environment and human health have been the subject of research [1,2]. Although light is giving ecosystem services to all living beings in the earth for their regulation of cell cycles, development, and growth; longer exposure to artificial CIMlighting has serious consequences [3]. It has been pointed out that longer exposure to ALAN impacts ecosystem functioning and causes strong hormonal responses, which in turn alters physiological processes, daily routines, and other life characteristics [2]. Further, it was also evident that ALAN disrupts invertebrates through changes in their community structure, life cycle, and species interaction, leading to an effect on food web dynamics across ecosystem boundaries, nutrient flow, pollination, and natural control of pest insects [4]. More recently, it has been indicated that ALAN may have contributed to the emergence of COVID-19 [5]. Overall, these investigations demonstrate that ALAN disturbs seasonal or circadian rhythms, endangering every living organism, and eroding biodiversity and wildlife, leading to decreased ecosystem services from the global ecosystem or biosphere.
In the last decade, transcriptomic profiling was performed in a variety of situations, including hypoxia, salt, drought, environmental pollution, and UV radiation [6]. Few transcriptomic studies were carried out to check the repercussions of ALAN in several organisms like birds [7,8], mice [9], coral reefs [10], and common toad tadpoles [11]. The tadpole study revealed that ALAN downregulates the expression of most genes both during the day and at night, and that light at night affects physiological pathways for immunity and lipid metabolism [11]. These studies have significantly enhanced our understanding of the molecular machinery needed for organismal adaptability. Therefore, it was evident that employing a transcriptomic approach to examine the physiological changes through gene expression is one of the best ways to study the consequences of ALAN.
The biology of reproduction in many species, including fish, is influenced by photoperiod and circadian cycles. ALAN also alters the expression of melatonin, a hormone that is largely produced at night, interferes with the synthesis of sex hormones, energy metabolism, antioxidant defences, and immune system components [12,13]. Improper exposure to ALAN has a negative impact on the circadian system, causing short-term effects on sleep and cognition as well as long-term endocrine disruptions that lead to obesity, cardiovascular disease, mood disorders, diabetes, hormonal imbalances during pregnancy, and cancer [14-18]. It has been found that ALAN primarily affects the internal circadian timekeeping system by changing the expression of clock genes [15], as well as other genes whose expression is controlled by photoperiod and involved in processes like reproduction, immunity, and metabolic detoxification [12]. Therefore, a comprehensive study is immediately warranted to analyze the effect of ALAN on a living being.
Zebrafish (Danio rerio) is a reliable vertebrate model for circadian research, originally developed for their utility in developmental biology studies, due to their high fecundity and relatively low cost of maintenance, experimental use of this model has evolved to include a broad range of phenotypes from embryonic to adult stages as well as a wide array of molecular assays [19]. Research on the effects of ALAN on zebrafish reveals multifaceted impacts spanning ovarian physiology, cognitive functions, and transgenerational behavior. In our previous study [15], we uncovered significant chronodisruption in zebrafish ovarian physiology under continuous light exposure, leading to desynchronized clock-associated genes and reduced melatonin levels, ultimately linked to ovarian tumor development. Lucon-Xiccato et al. [20] demonstrated ALAN’s detrimental effects on zebrafish learning abilities, slowing habituation learning, and altering behavioral covariances. Meanwhile, Li et al. [21] highlighted ALAN’s transgenerational impact, showing short-wavelength ALAN’s ability to induce anxiety-like behaviors and affect offspring behavior despite no direct exposure. These findings collectively underscore the profound implications of light pollution on zebrafish reproductive health, cognitive functions, and behavioral traits, urging further attention to mitigate its broader ecological consequences. Moreover, the zebrafish genome was fully sequenced in 2013 and contains homologues of around 70% of human genes, including 80% of those implicated in diseases, which makes zebrafish a reasonable model to study human diseases [22].
This enhances the ability to conduct human comparative transcriptome-based research using zebrafish with proper genomic annotations and identifiers. Previously our group employed a transcriptomic approach for the first time to show the development of thecoma and granulosa cell tumor, induced by a few candidate genes in ALAN-exposed zebrafish or any vertebrate ovary for the first time [15].
This study aims to demonstrate the level of differentially expressed genes (DEGs), the gene and pathway ontology, and diseases, through meticulous bioinformatics analysis of the complete ALAN-exposed ovarian whole transcriptome data.
METHODS
Animals and housing
The third generation of wild-type zebrafish around 6–7 months old, with a body length of 4.0±0.3 cm and weight of 0.4±0.15 g was utilized from the Institute of Bioresources and Sustainable Development (IBSD) zebrafish facility of, Imphal, Manipur, India. Fish were maintained in 50-L glass aquaria (3 fish per aquarium) under normal (12 h light:12 h dark) photic conditions (light intensity was fixed at 300 lux by a household fluorescent tube) [15]. Everyday light was on in the morning (at 6:00 am) and turned off in the evening (at 6:00 pm), maintained by the timer (Frontier digital timer, New Taipei City, Taiwan). The adequate water temperature for zebrafish (28°C±0.50°C) was maintained by using glass submersible aquarium immersion heaters (100 W, RS Electrical, Noida, India) placed in each aquarium. A biological filter (E-Jet, P.R.C) was used for the aeration and recirculation of water. The pH, hardness, and other parameters of water were maintained under standard conditions [23]. Food was given thrice a day. In the morning (9:00 am; ZT03) and midday (1:00 pm; ZT07) with commercial floating type small pellets (Perfect Companion Group Co. Ltd., Thailand). Live Artemia nauplii (cultured from Artemia cysts, Ocean Star International, Snowville, UT, USA) was given in the evening (5:00 pm; ZT11). Fish care and study schedule were done by following international standards [24]. Ethical clearance was obtained from the Institutional Animals Ethical Committee constituted as per the recommendations of the Committee for Control and Supervision of Experiments on Animals (CPCSEA), Government of India (IBSD/AE/2013/AB/001/01).
Experiment designs and sampling
As mentioned earlier, the adult female zebrafish (identical in age, body weight, and gonadal status) were randomly divided into four experimental groups in different photic conditions: one normal photoperiod (LD=12 h light:12 h dark) and three continuous light conditions 1) continuous light for one week, LLW (7 days); 2) for one month, LLM (30 days); and 3) for one year, LLY (365 days) [15]. All the samplings were performed at same time (12:00 PM). The ovary was removed by dissecting the body cavity respectively [25-27]. Sample for whole transcriptome sequencing was stored in RNAlater (Life Technology, Carlsbad, CA, USA), were stored at -80°C.
Isolation, qualitative and quantitative analysis of total RNA
Total RNA was isolated from the sample using TRIzol (Invitrogen, Carlsbad, CA, USA) method for RNA extraction. The quality of the total RNA was checked on 1% denatured agarose gel and quantified using a Qubit fluorometer.
RNA sequencing with Illumina 2×150 PE library
SENSE mRNA-Seq Library Prep Kit V2 (Lexogen, Vienna, Austria) was used as per the protocol to generate the libraries. This kit contains an extremely sensitive bead-based poly(A) selection phase which can eliminate even traces of rRNA, tRNA, and nonpolyadenylated RNA. After the selection of the desired RNA, the magnetic bead-bound poly(A) RNA was randomly hybridized with the starter/stopper heterodimers. These heterodimerize with Illumina-compatible linker sequences. One-step reverse transcription coupled with ligation reaction extended one starter to the other hybridized heterodimer. Now, a newly synthesized cDNA insert was ligated to the stopper. RNA fragmentation was not required because the insert size is dependent on the distance between starter/stopper binding sites. It eliminated spurious secondstrand synthesis which arose from the 5' ends of fragments and gave a platform for the excellent strand-specificity. Finally, a second strand was synthesized to release the library from the oligodT beads, and magnetic beads again purified the library, followed by polymerase chain reaction (PCR) amplification, and finally, the complete sequences required for cluster generation were introduced.
RNA-Seq data processing and identification of DEGs
Raw paired-end FASTQ files containing the RNA-Seq data were processed for a quality check using FastQC (http://www. bioinformatics.babraham.ac.uk/projects/fastqc/). Low-quality reads and adaptor sequences were removed using Trimmomatic (v.0.33) [28]. For genome-guided alignment, an index was created using Danio rerio (GRCz11) reference genome for HISAT2 (v2.2.0) [29], and alignment was performed. The aligned reads were then assembled by Cufflinks (v2.2.1) [30]. Based on our experiment and objective, we have created three contrasts i.e., LLW vs. LD, LLM vs. LD, and LLY vs. LD for the detection of differentially expressed genes/transcripts. Fragments per kilobase of exon per million mapped fragments (FPKM) values which reflect the normalization of the raw data for reading depth and gene length, were utilized to detect and quantify the DEGs by using the Cuffdiff tool of the Cufflinks package. A parametric distribution test to detect significant differences if they truly exist in the transcripts was checked using R software (v4.2.1; R Foundation for Statistical Computing, Vienna, Austria) with Tukey’s honestly significant difference (HSD) test based on log2(FPKM+1) means of all samples. Overall genomic features were assigned to the transcripts by using the featureCounts program [31] with Danio rerio gene structure from Ensembl GRCz11 v106 annotation. R software with ‘dplyr’ (v1.0.10) and ‘ggplot2’ (v3.3.5) libraries were used with log2FC (fold change) value ≥|2| (two-fold change in expression, either up or down-regulated) and significant p-value (p≤ 0.05) threshold for the identification and depiction (volcano and Venn plots) of significantly dysregulated DEGs across contrasts. The DEGs which were exclusively present in one sample but absent in others (genes with on/off signal) were also obtained for downstream analysis. Further, a principal component analysis (PCA) was also computed with blind normalization of the top 1,000 most variable transcripts to inspect the forming of clusters among the samples.
Network construction, functional module analysis, and identification of hub genes
Respective significant dysregulated DEGs of each contrast were imported into Cytoscape (v3.9.1) by using the stringApp (v1.7.1) plugin and with STRING (v11.5) (https://string-db.org/) database for creating the protein-protein interaction (PPI) network. Functional module/cluster analysis was performed using molecular complex detection (MCODE v2.0.2) [32] to get the densely connected regions and seed nodes/genes (the genes from where the cluster originates) within the PPI network. Some clusters had more than one seed and some had none, hence all seeds were taken from all the clusters of all three contrasts to analyze their expression changes in all three contrasts and depicted through a heatmap. For the calculation of hubs (genes with high connectivity) in the network, the cytoHubba (v0.1) plugin of Cytoscape was used. This tool helps in exploring biological networks through 13 different topological algorithms to identify important nodes including hub nodes and how they are interconnected with other genes. In our study, the top 50 hubs were predicted from each of the 13 algorithms out of which 4 are local-based methods: 1) Degree method (Deg), 2) Maximum Neighborhood Component (MNC), 3) Density of Maximum Neighborhood Component (DMNC), 4) Maximal Clique Centrality (MCC); and 7 are global-based methods: 1) Closeness (Clo), 2) EcCentricity (EC), 3) Radiality (Rad), 4) BottleNeck (BN), 5) Stress (Str), 6) Betweenness (BC), 7) Edge Percolated Component (EPC), details of which are mentioned by Bader and Hogue [32], and selected the hubs resulted in more than any of the 6 algorithms for downstream analysis. An interaction network of all the hubs predicted in each contrast was created using the STRING database and described through an organic layout for downstream analysis. To see the change in gene expression of all the hub and seed genes in all three contrasts a log2FC based heatmap was also constructed.
Pathway annotation and enrichment
Pathway annotation was carried out by taking 145,943 total DEGs to the Kyoto Encyclopedia of Genes and Genomes (KEGG) Automatic Annotation Server also known as KAAS, through BLAST comparisons against the KEGG GENES database. The BBH (bi-directional best hit) option was used to assign KO (KEGG orthology) terms (representative gene data set Danio rerio). The KO database (http://www.genome.jp/kegg/ko.html) was used for pathway mapping. Enrichment of the pathways in significantly dysregulated DEGs of different contrasts was performed using Database for Annotation, Visualization, and Integrated Discovery (DAVID v2021; https://david.ncifcrf.gov/) and REACTOME pathway database (https://reactome.org/), where all the genes of Danio rerio were taken as default background. The REACTOME database is additionally considered due to its core data model of reactions, where different entities participate to form a network of biological interactions and pathway groups. Each contrast’s significant DEGs were mapped separately with the pathways and followed the enrichment statistics of Fisher’s exact test and Benjamini-Hochberg false discovery rate (FDR). The enriched pathways were predicted based on the threshold of p<0.05, gene count >2, and FDR <0.05.
Gene ontology enrichment
Gene ontology (GO) term enrichment was done using R (v4.2) Bioconductor (v3.15) packages ‘clusterProfiler’ with Danio rerio OrgDb (“org.Dr.eg.db”) for all the respective dysregulated DEGs from all the contrasts. First, the DEGs were converted to Entrez gene ID (not found for less than 5% DEGs from all contrasts) and a default Benjamini & Hochberg adjusted p-value of 0.05 was used to identify the significant GO terms. Enrichment analysis was performed for a combination of all three GO categories, biological process (GO:BP), cellular component (GO:CC), and molecular function (GO:MF) as well as individual categories. Enrichment plots with the top 15 GO terms for combined categories and top 20 GO terms for individual categories were produced using built-in functions. Further, directed acyclic graphs (DAG) were also created for all three individual GO categories of each contrast to represent the hierarchical structure of the GO terms in them.
Disease association and enrichment
Following Howe et al. [33] zebrafish reference genome, 70% of protein-coding human genes are related to genes found in the zebrafish, and 84% of genes known to be associated with human disease have a zebrafish counterpart. So, we used the DisGeNET (v7.0) discovery platform, which is the largest publicly available collection of genes associated with human diseases [34], for enrichment analysis. All the dysregulated DEGs were subjected to manual curation to find their respective human orthologous using the ZFIN database (https://zfin.org/) March 2021 update and taken separately for their respective disease association and enrichment analysis using the R package ‘disgenet2r’ (v0.99.2) provided by DisGeNET. The same process was repeated for hub genes too. Few genes in our DEG dataset were automatically unused due to their absence in the DisGeNET collection. The complete disease dataset collection present at DisGeNET was used for genedisease class mapping, and heatmaps were generated for the respective gene-to-disease classes (collection of diseases of the same kind) using built-in functions. A set of genes associated with less than five diseases in all the contrasts were further taken to create a signature. Disease enrichment was performed against the whole collection of genes present at DisGeNET as background for the Fisher test. The p-values resulting from the multiple Fisher tests are corrected for false discovery rate using the Benjamini-Hoch-berg method. Significant diseases were marked with a p-value threshold of 0.05.
RESULTS
Whole transcriptome profiling and DEGs
High-quality paired-end reads for control (LD: 28,529,534 reads) and treatments (LLW: 26,459,892; LLM: 28,420,636; LLY: 28,541,418 reads) were obtained from the RNA sequencing experiment and found 88.62% (LD), 87.53% (LLW), 88.60% (LLM), 88.20% (LLY) mapped reads to Danio rerio (GRCz11) reference genome and a total of 145,943 genes with FPKM >0 were found in all the samples. The FPKM distribution showed the sample quality and their grouping by the way of mean differences (ANOVA; p<0.0001) through Tukey’s HSD test (at 95% confidence level) creating two groups, one for control sample LD, and the other for the treated samples LLW, LLM, and LLY (Supplementary Figure 1). Altogether, these results stipulated that the quality of the data generated could satisfy further analysis. The genomic features predicted in the expressed transcripts showed so many other overlapping biotypes with the annotated zebrafish genome (Supplementary Figure 2). The PCA of the expressed transcripts showed that LLM and LLW samples formed a cluster, whereas LD and LLY samples each were isolated separately, far from all other samples (Supplementary Figure 3). This suggests the effect of ALAN is different depending on the amount of exposure to each sample. Nonetheless, we identified 777, 758, and 802 significantly dysregulated transcripts in LLW, LLM, and LLY groups, respectively, in comparison to the LD sample based on the thresholds of p<0.05 and a log2FC≥|2| (Figure 1A). A scaled heatmap of all DEGs from each contrast in comparison to the other two contrasts showed major differences in their expression (Supplementary Figure 4). Most of the genes were found to be significantly upregulated in all three treated groups whereas the number of dysregulated genes were higher in LLY vs. LD contrast compared to the other two. The union of significantly dysregulated DEGs showed 311 commonly expressed genes in all three contrasts, whereas in paired comparison we found 141 in LLW vs. LD and LLM vs. LD, 97 in LLM vs. LD and LLY vs. LD, and 104 in LLW vs. LD and LLY vs. LD, respectively (Figure 1B). The sample-wise trend in significantly dysregulated DEGs showed a decrease in downregulated genes towards continuous light exposure over a year. For upregulated DEGs, we observed a U-shaped trend over time: higher at LLW (446) and LLY (504) with a dip at LLM (444), suggesting a dynamic response to continuous light exposure in comparison to LD sample (Figure 1C).
Differential expression analysis. A: Volcano plot of differentially expressed genes (DEGs) with a threshold of log2FC ≥|2|, p≤0.05 of (i) LLW vs. LD, (ii) LLM vs. LD, and (iii) LLY vs. LD samples. B: Union of significantly dysregulated DEGs. C: Sample wise trend in significantly dysregulated DEGs.
DEG network analysis
The PPI network for each set of dysregulated DEGs was constructed using the STRING database in Cytoscape and their functional cluster analysis was carried out. Almost all the seed genes were found to belong to ribosomal proteins. No common seed genes were found in all three contrasts, but we found zgc:56596 (orthologous to human UBL4A [ubiquitin-like protein 4A]) and zp2.6 (orthologous to human ZP4-zona pellucida glycoprotein 4) genes common in LLW vs. LD and LLM vs. LD contrasts. The number of clusters formed and all seed genes of all the clusters were mentioned in Supplementary Table 1, whereas their network pathway enrichment was presented in Supplementary Table 2. The pathways enriched were found same for all three contrasts as ribosome, protein export, oxidative phosphorylation, metabolic pathways, and cardiac muscle contraction, but in LLY vs. LD, two extra pathways were enriched i.e., RNA polymerase and cytosolic DNA-sensing pathway. Further, the prediction of hubs using cytoHubba resulted in 49 (6 up, 43 down), 44 (8 up, 36 down), and 52 (12 up, 40 down) hub genes in LLW vs. LD, LLM vs. LD, and LLY vs. LD contrasts, respectively, shows most of them are downregulated. In all three contrasts, most hub genes enriched in ribosome and oxidative phosphorylation pathways, but only in LLY vs. LD hub genes, we find three extra pathways enriched, namely, RNA polymerase, adherence junction, and FoxO signaling pathway (Supplementary Figure 5). Further, the gene expression in both hub and seed genes in terms of log2FC have shown drastic changes in all three contrasts (Figure 2). These results show a significant change in gene expression level in long exposure to ALAN in LLY samples compared to LD samples.
Heatmap of log2FC values for hub genes (A) seed genes (B).
Functional association and enrichment analysis
The impact of ALAN on the biological processes was analyzed by annotation and enrichment tools/servers like KEGG, DAVID, and R Bioconductor packages. Although the total DEGs (145,943) mapping to the KAAS server showed 41,177 pathway terms, only the terms related to five major classified functional hierarchies namely, metabolism, genetic information processing, environmental information processing, cellular processes, and organismal systems were taken for the analysis, rest unclassified terms are discarded. Amongst the classified pathway terms (Figure 3), more than 500 counts were found in signal transduction (1,901), endocrine system (952), immune system (915), and signaling molecules and interaction (905), suggesting that signal transduction and endocrine system systems were stimulated due to ALAN exposure.
Pathway annotation of all differentially expressed genes (DEGs) using KEGG Automatic Annotation Server (KAAS), different color represents different groups of pathways mapped.
The pathway enrichment analysis for significantly dysregulated DEGs in DAVID showed oxidative phosphorylation, ribosome, RNA polymerase, ubiquitin-mediated proteolysis, RNA degradation, spliceosome, and salmonella infection pathways common to the three contrasts (Figure 4A). Homologous recombination, and glycosylphosphatidylinositol (GPI)-anchor biosynthesis pathways, were enriched only in LLW vs. LD DEGs. 2-Oxocarboxylic acid metabolism, citrate cycle (TCA cycle), biosynthesis of amino acids, mitophagy-animal, glycerophospholipid metabolism, carbon metabolism, cell cycle, and protein processing in the endoplasmic reticulum were only enriched in LLM vs. LD. Finally, in LLY vs. LD, notch signaling, lysine degradation, and vascular smooth muscle contraction pathways were found enriched. Further, the REACTOME-based pathway analysis showed signaling-related pathways were more enriched in LLY vs. LD contrasts pathway (Supplementary Figure 6). The separate enrichment test for genes exclusively present in one sample and absent in another were shown many pathways were absent in the control sample LD but present in treated samples, especially metabolic pathways are the main term found enriched only in LD sample and not at all enriched in the treated samples (Figure 4B).
Pathway enrichment comparative analysis using DAVID server. A: Contrast-wise significantly dysregulated differentially expressed genes (DEGs). B: DEGs exclusively present in one sample only while absent in other, also called the signal off/on genes.
The GO categories were widely used to describe or classify the genes according to three aspects, biological process (GO:BP), cellular component (GO:CC), and molecular function (GO:MF). GO:BP corresponds to the specific purposes that the organism is genetically designed to achieve, like translation or signal transduction, etc. GO:CC represents the sites relative to cellular structures (like mitochondrion compartments, or stable macromolecular complexes [ribosome]) in which a gene product carries out a function. GO:MF terms describe molecular-level activities, such as transport or methyltransferase or binding activity performed by gene products. The enrichment of GO terms in significant dysregulated DEGs of LLW vs. LD and LLM vs. LD contrasts showed a large number of biological processes in comparison to LLY vs. LD, and the opposite was found for molecular function which is more in LLY vs. LD and less in LLW vs. LD and LLM vs. LD (Figure 5). Apoptotic signaling pathway and homeostasis of number of cells are enriched in in LLM vs. LD and LLY vs. LD contrasts but not in LLW vs. LD. Terms related to ATP metabolism or synthesis are absent in the year-long exposed sample (Figure 5). Further, the DAG plot of GO terms showed more biological processes in LLW vs. LD followed by LLM vs. LD and LLY vs. LD. Conversely, the enriched cellular components and molecular functions were higher in LLY vs. LD and LLM vs. LD than LLW vs. LD (Figure 6).
Gene ontology (GO) enrichment with separated GO categories (BP, CC, MF) of significantly dysregulated DEGs of each contrast with top 20 enriched terms. A: LLW vs. LD. B: LLM vs. LD. C: LLY vs. LD contrasts.
Directed acyclic graph of GO:BP, GO:CC, and GO:MF enrichment of significantly dysregulated differentially expressed genes (DEGs). A: LLW vs. LD. B: LLM vs. LD. C: LLY vs. LD contrasts.
Disease association and enrichment analysis
Gene-disease association (GDA) and disease enrichment in the significantly dysregulated DEGs were studied using the current DisGeNET collection of genes and diseases (v7.0), which contains 1,135,045 GDAs, between 21,671 genes and 30,170 diseases, disorders, traits, and clinical or abnormal human phenotypes. The whole collection of diseases was made the background for the enrichment analysis. Although in the zebrafish genome, 70% of protein-coding human genes are related to genes found in the zebrafish, manual curation was done to all the DEGs to find their counterpart orthologous in humans using the Zebrafish Information Network (ZFIN) database. Many disease classes were found associated with the DEGs (Supplementary Tables 3-5), out of them, disease group “neoplasms” which deals with cancer types, and phenotype pathological conditions, signs, and symptoms classes were found to be abundant in all the contrasts (Figure 7). The enrichment analysis showed “malignant neoplasm of prostate” as the top enriched disease in all contrasts with a significant change in ratio towards the year data (LLW vs. LD: 0.25, LLM vs. LD: 0.25, LLY vs. LD: 0.3). Prostate carcinoma, tumor progression, melanoma, and other cancer-related diseases were found common in all top enriched diseases in all three contrasts, but the disease “malignant neoplasm of lung” is only enriched in LLY vs. LD contrast (Figure 8).
Gene-MeSH disease class heatmap of significantly dysregulated differentially expressed genes (DEGs) of LLW vs. LD (A), LLM vs. LD (B), and LLY vs. LD contrasts (C).
Disease enrichment analysis of significantly dysregulated differentially expressed genes (DEGs) of LLW vs. LD (A), LLM vs. LD (B), and LLY vs. LD contrasts (C).
Furthermore, a merged set of all DEGs from all the contrasts was also tested to check overall enriched diseases, which showed, mammary neoplasms, neuroblastoma, malignant neoplasm of the pancreas, intellectual disability, carcinoma, and other cancerrelated diseases to be the top enriched diseases (Figure 9A). We have also marked a gene set of 28 genes, which were associating themselves with less than five diseases in all the contrasts. Their log2FC values were plotted as a heatmap to analyze their expression changes (Figure 9B). We believe these 28 genes can be tested in the future for a valuable prognosis marker for the diseases developed (Supplementary Table 6). Additionally, the enrichment of hub genes for probable diseases showed similar kinds of diseases in LLW vs. LD and LLM vs. LD contrasts but found “carcinogenesis” largely enriched only in the LLY vs. LD contrast (Figure 10).
Disease enrichment analysis of all differentially expressed genes (DEGs). A: Disease enrichment analysis of merged total significant dysregulated DEGs of all contrasts. B: Heatmap with log2FC values of significant DEGs encoding ≤5 diseases during the disease enrichment analysis.
Disease enrichment analysis of hub genes. A: LLW vs. LD. B: LLM vs. LD. C: LLY vs. LD contrasts.
DISCUSSION
Our study comprising four transcriptomes of zebrafish ovary exposed to different photic conditions showed that ALAN induces a strong disruption of gene expression starting from one day to a year. Earlier studies indicated that ALAN disrupts sleep and circadian rhythm, raising the risk of cancer, diabetes, cardiovascular disease, obesity, mood disorders, and age-related macular degeneration [18]. A recent study employing the RNA-Seq method to study tadpole larvae exposed to ALAN found that their capacity to deal with immunity is reduced [11]. The impact of ALAN on the transcriptome level has not yet been well investigated in zebrafish, particularly biological rhythms, (controlled transcriptionally) [35]. In this study, the full spectrum of diseases that may develop from exposure to ALAN, is elaborated for the first time (Supplementary Tables 3-5). This helps to predict sustainability of this organism in the ecosystem under rapidly changing environmental parameters. The U-shaped trend in upregulated DEGs suggests an intriguing aspect of light exposure. Initially, LLW shows an acute stress response with high upregulation. At LLM, organisms begin adapting, reducing the number of upregulated genes as they seek homeostasis. By LLY, chronic stress from prolonged exposure increases gene upregulation again, reflecting long-term effects and potential cumulative cellular damage. This pattern highlights the complex, nonlinear nature of biological responses to extended environmental changes, such as continuous light exposure disrupting circadian rhythms and inducing varied stress responses over time. However, the majority of DEGs were found upregulated in all three ALAN exposed conditions in comparison to the normal condition, leading to changes in many physiological pathways like translation, signaling, metabolism, ATP synthesis, methyltransferase, and oxidase activities, causing overall chronodisruption in zebrafish ovary (Figures 3-5 and Supplementary Figure 6). Pathways like RNA degradation, RNA polymerase, signaling, and ABC transporters began enriched in the one-week sample and accelerated over one year of exposure (Figure 4A). The Notch signaling system, which controls cellular differentiation, survival, apoptosis, and proliferation [36], and the Lysine degradation pathway which controls the inborn errors of metabolism [37], were found to be dysregulated in the year-long ALAN exposed samples (Figure 4A). This might have a significant impact on the ability of the organ to function properly. Interestingly these were not the top enriched pathways in the week and month-long ALAN-exposed samples (Figure 4A). The decrease in several genes contributing to oxidative phosphorylation following the year-long exposure to ALAN reveals that the energy needed to maintain metabolism, the physical structure of the cells, and the transportation of molecules and ions are affected and disrupted (Figure 4A). The previous study [38] demonstrated this disruption as the result of the intracellular environment. This was further confirmed by the enrichment of the ATP metabolic process in week and month samples but not in the LLY sample (Figures 4 and 5). Moreover, the activation of the apoptotic signaling pathway and homeostasis in the number of cells in LLM and LLY samples shows that cell death is evident and, hence disrupting the total organ (Figure 5). The GO enrichment provided additional support for the findings of pathway enrichment, indicating a decline in significant biological processes and an increase in significant molecular functions as the ALAN exposure was continued (Figure 5). Furthermore, genes expressed specifically in treated or normal samples demonstrated a disruption in either metabolic pathways or signaling pathways, respectively (Figure 4B). These two processes are bidirectional, and therefore, cells can keep track of the conditions of crucial metabolic pathways and exert feedback control over signal transduction networks [39], allowing the cells to coordinate their metabolic processes with their growth and division. The disruption of these pathways amply demonstrates that the organ is progressing toward a serious diseased condition. In addition, the treated samples showed enrichment in the TGF-beta and GnRH signaling pathways, indicating significant changes in growth and reproduction processes. TGF-beta signaling has complex roles; it can promote growth, inhibit cell multiplication, trigger cell specialization or death, and give rise to various diseases, depending on the situation [40]. Meanwhile, GnRH signaling is crucial for reproductive function, affecting the release of hormones that regulate reproduction (Figure 4B). Our previous findings [27] also support this data, moreover, the effect of ALAN on the brain-pituitary-gonadal axis implies the degree of lethality of this exposure. According to the hub genes pathway enrichment of adherence junction and FoxO signaling system only happened in the upregulated genes of the year-long exposed sample. ALAN’s negative effects were further reduced by preventing cell cycle arrest, stress tolerance, and apoptosis (Supplementary Figure 5). Disruption of the FoxO family of transcription factors results in the sick condition of the ovary because they are essential for the integration of growth factor signalling, oxidative stress, and inflammation [36,41].
The LLY group, being older due to 12 months of continuous light exposure, may likely exhibits age-related gene expression changes. These changes may interact with the effects of ALAN, contributing to the observed dysregulation of genes associated with neoplasms. As zebrafish age, there is a documented disruption in the expression of genes critical for stem cell maintenance and tissue homeostasis, which can lead to increased susceptibility to diseases, including various cancers [42,43]. Though zebrafish are considered as adult from 90 days to 2 years old, and the age of the fish in the present study falls within this range, future research with age-appropriate control might be useful in understanding the exact DEGs [44]. In our study, the significant disruption of various DEGs and pathways in the treated samples indicated the emergence of numerous clinical states, therefore a disease annotation and enrichment analysis using human orthologous data in the entire DisGeNET database was carried out. Previously it was studied that, exposure to ALAN suppresses melatonin production, disrupts circadian rhythms, and is linked to increased risks of hormone-related cancers such as breast, prostate, and colorectal cancers due to its effects on biological pathways like hormone signaling, cell proliferation, DNA repair, and inflammation [45]. Additionally, ALAN is significantly correlated with increased risks of all forms of cancer, including lung, breast, colorectal, and prostate cancers, across 158 countries. This study found positive associations between ALAN exposure and cancer rates, even after adjusting for confounders such as population, electricity consumption, air pollution, and forest coverage [46]. A global study further highlighted higher breast cancer incidence in countries with elevated light pollution levels, suggesting that light pollution might disrupt circadian rhythms and melatonin production, thereby contributing to cancer risk [47]. As anticipated, neoplasms are the most common type of sickness class annotated, followed by the phenotype of pathological states, signs, and symptoms (Figure 7). The most prevalent disorders enriched in all treated samples included malignant prostate cancer (which was previously confirmed by a study [48]), glioblastoma, melanoma, tumor progression, and other cancer-related diseases confirming our previous study [15]. We have taken all significantly dysregulated DEGs in all the contrasts to get a comprehensive picture of the disease enrichment in the treated samples (Figure 9A). We found major classes of neoplasms in the mammary and pancreas, neuroblastoma and glioblastoma in both children and adults, carcinoma, along with intellectual disability, global developmental delay, and seizures. This shows that ALAN’s effects are worse than previously thought. Additionally, genes that were all annotated with less than 5 diseases in the DisGeNET database were chosen to create a signature of genes whose change in expression level may result in these types of diseased states (Figure 9B). It was also noted that, two genes, zgc:56596 and zp2.6, were commonly dysregulated in both the LLW vs. LD and LLM vs. LD contrasts. The zgc:56596 gene is orthologous to the human UBL4A gene which is known to play a role in protein degradation, a critical process for maintaining cellular homeostasis and responding to stress [49,50]. Its dysregulation suggests that ALAN exposure might disrupt protein homeostasis and degradation pathways, contributing to cellular stress responses observed in our study. Similarly, the zp2.6 gene, orthologous to the human ZP4 gene, is involved in reproductive processes, specifically in the structure of the zona pellucida, which surrounds the oocyte and is essential for fertilization [51]. The dysregulation of zp2.6 highlights potential impacts of ALAN on reproductive health, aligning with other findings in our study that suggest reproductive processes are significantly affected by continuous light exposure.
In our transcriptomics analysis of four ALAN exposure scenarios, we identified 777 dysregulated transcripts in the LLW group, 758 in the LLM group, and 802 in the LLY group, compared to the LD sample. The functional enrichment analysis revealed that signaling pathways were more prevalent while ATP-related processes became scarce in the year-long ALAN exposure. Several disease groups that affect the basic functioning of the organism, such as the neural, immune, endocrine, reproductive, digestive, lymphatic, and cardiovascular were recognized by the transcriptomic changes. Major diseases including reproductive problems, global developmental delay, seizures, dwarfism, and cancer (prostate, mammary, lung, kidney, glioblastoma, melanoma, and lymphoma) were found to be triggered by ALAN. It is important to note that while these associations were observed in our study are based on zebrafish models and require validation in human studies to establish direct correlations. Although further studies are required, our study concluded that there is a possibility of ALAN globally affecting the physiology of the organism, initiating various disease-inducing pathways, and may cause numerous diseases by affecting the global gene expression pattern in zebrafish ovary.
Supplementary Materials
The online-only Data Supplement is available with this article at https://doi.org/10.33069/cim.2024.0009.
Sample comparable quality estimation using the FPKM distribution (log2(FPKM+1)), predicted through Tukey’s HSD test showing significant difference between control (LD) and treated samples (LLW, LLM, and LLY).
Transcript feature count revealing different biotypes overlapping with the zebrafish genome.
PCA plot of top 1,000 most variable transcripts.
Heatmap of all significant DEGs, from the point of view of 777 DEGs of LLW vs. LD with other contrasts (A); 758 DEGs of LLM vs. LD with other contrasts (B); and 802 DEGs of LLY vs. LD with other contrasts (C).
Interaction network of hub genes and the enriched pathways in them. LLW vs. LD (A), LLM vs. LD (B), and LLY vs. LD contrasts (C).
Pathway enrichment comparative analysis using REACTOME database.
Details of cluster analysis of dysregulated DEGs
Pathway enrichment of top two clusters of significant dysregulated DEGs of each contrast, other clusters didn’t enrich any pathways
LLW vs LD disease association table, with each sheet containing the encoded diseases
LLM vs LD disease association table, with each sheet containing the encoded diseases
LLY vs LD disease association table, with each sheet containing the encoded diseases
28 gene showing less than 5 disease association in all the contrasts presented sheetwise
Notes
The authors have no potential conflicts of interest to disclose.
Availability of Data and Material
The datasets generated during and/or analyzed during the current study are available from the corresponding author upon reasonable request. The RNA-Seq raw files are submitted at ArrayExpress with accession no. E-MTAB-13218.
Author Contributions
Conceptualization: Asamanja Chattoraj. Data curation: Rajendra Kumar Labala. Formal analysis: Rajendra Kumar Labala, Zeeshan Ahmad Khan, Gopinath Mondal. Investigation: all authors. Methodology: Rajendra Kumar Labala, Zeeshan Ahmad Khan. Project administration: Asamanja Chattoraj. Resources: Asamanja Chattoraj. Software: Rajendra Kumar Labala. Supervision: Asamanja Chattoraj. Validation: Asamanja Chattoraj, Rajendra Kumar Labala, Zeeshan Ahmad Khan, Gopinath Mondal. Visualization: Rajendra Kumar Labala. Writing—original draft: Asamanja Chattoraj, Rajendra Kumar Labala, Zeeshan Ahmad Khan, Gopinath Mondal. Writing—review & editing: all authors.
Funding Statement
None
Acknowledgements
The authors are grateful to the Director, IBSD, India, and DBT (BT/407/NE/U-Excel/2013) for funding. The authors are also thankful to Mr. L. Umakanta Singh, Mr. S. Surjit Singh, Mr. W. Rahul, and Mr. H. Bishorjit for their continuous support in collection and maintenance of fish in laboratory conditions. Authors are thankful to the anonymous reviewers for their valuable suggestions to improve the manuscript. Our sincere gratitude to every Indian taxpayer for their contribution towards the funding for research.