Poly(ADP-ribose) polymerase 1 in genome-wide expression control in Drosophila

Poly(ADP-ribose) polymerase 1 (PARP-1) is a nuclear enzyme involved in DNA repair and transcription regulation, among other processes. Malignant transformations, tumor progression, the onset of some neuropathies and other disorders have been linked to misregulation of PARP-1 activity. Despite intensive studies during the last few decades, the role of PARP-1 in transcription regulation is still not well understood. In this study, a transcriptomic analysis in Drosophila melanogaster third instar larvae was carried out. A total of 602 genes were identified, showing large-scale changes in their expression levels in the absence of PARP-1 in vivo. Among these genes, several functional gene groups were present, including transcription factors and cytochrome family members. The transcription levels of genes from the same functional group were affected by the absence of PARP-1 in a similar manner. In the absence of PARP-1, all misregulated genes coding for transcription factors were downregulated, whereas all genes coding for members of the cytochrome P450 family were upregulated. The cytochrome P450 proteins contain heme as a cofactor and are involved in oxidoreduction. Significant changes were also observed in the expression of several mobile elements in the absence of PARP-1, suggesting that PARP-1 may be involved in regulating the expression of mobile elements.


Distribution of DEG types between the wild type and parp-1 −/− groups.
To determine the role of PARP-1 in gene expression regulation of DEGs, we first normalized the expression level of every DEG to the average expression level of all DEGs. The wild-type group was normalized by the average expression level of all DEGs in the wild-type group, and the parp-1 −/− group was normalized by the average expression level of all DEGs in parp-1 −/− . We then divided all DEGs into six clusters based on the difference of their expression level between the wild type and parp-1 −/− groups (Fig. 3A, Table 1). This clustering sorts DEGs according to their expression profile in both groups. The first cluster (C1) consists of 68 genes that had a normalized expression lower than average in both groups and that were downregulated in the parp-1 −/− group. The second cluster (C2) Scientific Reports | (2020) 10:21151 | https://doi.org/10.1038/s41598-020-78116-5 www.nature.com/scientificreports/ includes 115 genes with a lower than average normalized expression in both groups and that were upregulated in the parp-1 −/− group. The third cluster (C3) comprised 94 genes that had a lower normalized expression than the average in the parp-1 −/− group, but higher than the average in the wild type group. The fourth cluster (C4) includes 97 genes that had a higher expression than the average in the parp-1 −/− group and lower than the average in the wild-type group. The fifth cluster (C5) includes 77 genes that had a higher expression than the average in both groups and were downregulated in the parp-1 −/− group. The sixth cluster (C6) comprises the remaining 151 genes that had a higher expression level than the average in both groups and were upregulated in the parp-1 −/− group. Only 183 out of 602 genes (30.4%) had a lower than average normalized expression in both the wild type and parp-1 −/− groups (C1 and C2), whereas 419 out of 602 (69.6%) had a higher than average expression in at least one of our groups (C3-6). Because it is mathematically easier to attain a large fold difference when the initial level The last column represents the fold difference between the wild type and parp-1 −/− groups. The fold difference was measured as the negative expression level in the wild type larvae divided by the expression level in the parp-1 −/− larvae when the expression was higher in the wild type group (shown in blue) and the expression level in parp-1 −/− larvae divided by the expression level in wild type larvae when the expression was higher in the parp-1 −/− group (shown in red). (C) Distribution of the 602 DEGs, depending on the fold difference of their expression level between the wild type and parp-1 −/− groups. www.nature.com/scientificreports/ of expression is low, the fluctuation levels in expression of the DEGs belonging to C1 and C2 may not necessarily be biologically relevant. Changes in the expression levels of DEGs belonging to C3 or C4 in the presence or absence of PARP-1 are the most biologically relevant because they present higher expression than average in one of the groups and lower expression than average in the other group. Finally, genes in C5 and C6 are also relevant because they present higher expression than average in both groups, making it mathematically harder to attain a significant fold difference between them. Clustering DEGs allows us to distinguish between the more and less biologically relevant DEGs based on changes in expression levels between the two groups.
Functional classification of DEGs. To classify DEGs according to their function and to confirm that the functions of the downregulated DEGs differ from those of the upregulated DEGs, we performed gene set enrichment analysis using GSEA software (see "Materials and methods" and Supplementary Fig. 1). We selected the most significant Gene Ontology (GO)-term that included at least 4 DEGs (Fig. 3B). Interestingly, we found little overlap between the gene ontology (GO)-terms overrepresented among the up-and downregulated DEGs. This means that only a small proportion (4.2%) of downregulated DEGs was involved in a GO-term overrepresented among upregulated DEGs and only a small fraction (8.5%) of the upregulated DEGs was involved in a GO-term overrepresented among downregulated DEGs (data not shown). We identified nine significant GOterms that mainly include downregulated DEGs (Fig. 3B We also found five significant GO-terms that mainly included upregulated DEGs. DEGs involved in 'GO:0055114 Oxidation-reduction process' were overrepresented among upregulated DEGs with 31 out of 44 DEGs (70.4%) belonging to C4 and C6, respectively. DEGs involved in 'GO:0020037 Heme binding' and 'GO:0005975 Carbohydrate metabolic process' were overrepresented among upregulated DEGs with 14 out of 22 DEGs (63.6%) and 5 out of 9 DEGs (55.5%) belonging to C4 and C6, respectively. All DEGs involved in 'GO:0006030 Chitin metabolic process' were upregulated in parp-1 −/− larvae with 8 out of 8 DEGs (100%) belonging to C4 and C6. Finally, DEGs involved in 'GO:0055085 Transmembrane transport' were overrepresented among upregulated DEGs with 9 out of 21 DEGs (42.8%) belonging to C4 and C6.
Taken together, these results show that several processes are differentially regulated in the absence of PARP-1. For example, DEGs involved in the GO-term 'DNA-binding transcription factor activity' were all downregulated in the absence of PARP-1, whereas DEGs involved in the GO-term 'Heme binding' were all upregulated in absence of PARP-1.  Table 2). Nine of these transcription factors are involved in neuronal development or axonal growth. Among the twelve misregulated transcription factors, we found three C2H2 zinc finger transcription factors; two of them were involved in neuronal development or axonal growth. Specifically, the gene coding for Longitudinals lacking (Lola) (cluster 1) and the gene coding for Castor (Cas) (cluster 3) were both downregulated in the absence of PARP-1. Lola is a transcription repressor involved in axonal growth and guidance in the Drosophila www.nature.com/scientificreports/ central nervous system (CNS) 18 . Lola is also important for programmed cell death in ovaries 19 and in cell fate in eyes by antagonizing Notch induction 20 . Castor is involved in the temporal patterning of Drosophila CNS 21 .
We also found four basic helix-loop-helix (bHLH) transcription factors involved in neuronal development. The gene coding for Achaete (Ac) (cluster 1), the genes coding for two members of the Enhancer of Split family, E(spl)m7-HLH (cluster 3) and E(spl)mα-BFM (cluster 5), and the gene coding for Cycle (Cyc) (cluster 3) were all downregulated in the absence of PARP-1. Ac plays a role in the formation of neural precursors 22 by serving as an antagonist to Notch signaling 23 . E(spl)m7-HLH and E(spl)m α-BFM are also involved in Notch signaling. More specifically, E(spl)m α-BFM plays a role in Notch lateral inhibition 24 , and E(spl)m7-HLH acts within the Notch pathway to repress neural fate 25 , but it has also been reported to interact with Ac 26 . Cyc is a circadian clock protein, which mediates several processes, including the olfaction rhythms of the antennal neurons 27 and the interconnection of feeding and sleeping behavior 28 .
Three other transcription factors are involved in neuronal development. All were downregulated in the absence of PARP-1, including the gene coding for the high mobility group box transcription factor Soxneuro (SoxN) (cluster 3), the gene coding for the homeodomain transcription factor Reversed polarity (Repo) (cluster 3), and the gene coding for the Nuclear factor of activated T-cells transcription factor (NFAT) (Cluster 5). SoxN plays a role in the specification of neural progenitors in CNS 29 . This transcription factor is also important in the regulation of Wg/Wnt signaling activity 30 . Repo has been reported to participate in the maintenance of glial fate in the Drosophila nervous system 31 . NFAT plays an important role in several processes, including neuronal development and plasticity 32 .
The three remaining transcription factors downregulated in the absence of PARP-1 have not been reported to play a direct role in neuronal development. These included the gene coding for the C4zinc finger liganddependent transcription factor Ecdysone-induced protein 78C (Eip78C) (cluster 3), the gene coding for the Table 1. Recapitulation of the composition of the different DEG clusters used in this study. Clusters 1 and 2 include the DEGs with lower expression compared to the other genes in both the wild type and parp-1 −/− groups. Clusters 3 and 4 include DEGs with lower expression compared to the other genes in one group and a higher expression in the other group. Clusters 5 and 6 include DEGs with higher expression compared to other genes in both wild type and parp-1 −/− groups. www.nature.com/scientificreports/ basic leucine zipper transcription factor Slow border cells (Slbo) (cluster 5), and cg31365 coding for a C2H2 zinc finger transcription factor (cluster 1). Eip78C plays a role in regulating chromosome puffs 33 and ovarian germline stem cells 34 . Eip78C has also been reported to physically interact with SoxN and Cyc 35 . Slbo has been reported to participate in cell migration during oogenesis 36 . The transcription factor for the protein coded by cg31365 relies on predicted functionality based on its sequence similarity with BCL6B, a human transcriptional repressor, but has never been confirmed in vitro 17 . Human BCL6B is a tumor suppressor that inhibits hepatocellular carcinoma metastasis 37 . We identified four additional genes involved in the transcriptional process. Although not classified as transcription factors, they were still misregulated in the parp-1 −/− group; all were downregulated. These included the gene coding for the Negative elongation factor E (Nelf-E) (cluster 3), the gene coding for the histone acetyltransferase Chameau (Chm) (cluster 5), the gene coding for the Mediator complex subunit 7 (MED7) (cluster 5), and the gene coding for Modifier of mdg4 (Mod(mdg4)) (cluster 5). Nelf-E is a member of the NELF complex involved in pausing RNA polymerase II at several promoters, including hsp70 promoter 38 . Nelf-E is involved in the activation of several key developmental genes during embryogenesis 39 . Chm is known to act along with the polycomb group to maintain hox gene silencing 40 . Chm also serves as a cofactor in the JNK/AP-1-dependent transcription during metamorphosis 41 . MED7 is a member of the mediator complex, which serves as a hub for transcriptional signaling events 42 . Modifier of mdg4 has 31 alternative splice products reported to participate in a range of processes, including chromosome segregation and synapse formation 43,44 .
Taken together, these results suggest that PARP-1 promotes expression of multiple transcription factors, as well as several genes mediating transcription. The absence of PARP-1 suppresses expression of several transcription factors important in neuronal development and axonal growth.
DEGs coding for cytochrome P450 are upregulated in the absence of PARP-1. Among DEGs, we identified 22 genes involved in the 'GO:0020037 Heme binding' GO-term (Fig. 3B). All but one were upregulated in the parp-1 −/− group compared to wild type. The expression of these upregulated genes in the absence of PARP-1 varied from 2.56 for Cyp9f3Ψ to 1575.97 for Cyp6a17. We determined that 18 out of these 22 DEGs (81.8%) belonged to the cytochrome P450 family, three were related to the cytochrome b5 family, and the last one coding for peroxidase Globin 1 (Glob1) (cluster 6) was unrelated to the cytochrome family 45 (Table 3).
Cytochrome P450 proteins contain heme as a cofactor and are involved in oxidoreduction. These proteins are important in the clearance of several compounds and in hormone synthesis and breakdown 46 . Among DEGs related to the cytochrome P450 family, 14 out of 18 (77.8%) have at least a predicted function in the breakdown www.nature.com/scientificreports/ of toxic chemicals. Based on their sequence, seven have been predicted to play a role in the metabolism of insect hormones and the breakdown of toxic chemicals, although this function has still not been confirmed in vivo 17,47 . These included Cyp4s3 (cluster 2), Cyp9b1 (cluster 2), Cyp6a17 (cluster 4), Cyp4d8 (cluster 4), Cyp6d4 (cluster 6), Cyp28a5 (cluster 6), and Cyp6a13 (cluster 6). Cyp6a17 has been reported to play an important role in temperature preference behavior 48 . Seven have been shown to be involved in breakdown of toxic chemicals. Six are involved in Dichlorodiphenyltrichloroethane (DDT) resistance, including Cyp6a8 (cluster 2) 49 , Cyp6a2 (cluster 2) 49 , Cyp6w1 (cluster 4) 50 , genes coding for Cyp4p1 (cluster 6) and Cyp4p2 (cluster 5) that are organized in a cluster on DNA 51 , and the gene coding for Cyp9c1 (cluster 6) 52 . Finally, the gene coding for Cyp4e3 (cluster 4) plays a role in permethrin insecticide tolerance 53 . The last four DEGs coding for a member of the cytochrome P450 family (Cyp12e1, Cyp12a5, Cyp12c1, and Cyp9f3Ψ) do not have any predicted functions aside from heme-binding. Three DEGs involved in the GO-term 'Heme binding' had a cytochrome b5 heme-binding site protein signature; however, these DEGs are not directly linked to the cytochrome P450 family, and all of them were upregulated in the absence of PARP-1. Cytochrome b5 proteins are hemoproteins involved in electron transport 54 . These three genes include the gene coding for the Fatty acid 2-hydrolase (Fa2h) (cluster 2), the gene coding for the Cytochrome b5-related (Cyt-b5-r) (cluster 6), and cg5157 (cluster 2). Fa2h plays an important role in larvae survivability when a correct amount of sterol is lacking 55 . Cyt-b5-r is predominantly expressed in muscles 56 . The gene cg5157 is predicted to possess a cytochrome b5 heme-binding site, but its function is unknown 17 .
Taken together, these results suggest that PARP-1 inhibits the expression of several cytochrome-related proteins. Most of them contribute to resistance against toxic chemicals, such as DTT or permethrin, suggesting that PARP-1 may contribute to the organismal response to toxins.
Mobile elements tend to be upregulated in the absence of PARP-1. We also found eight mobile elements ( Table 4, first panel) among DEGs. Three were downregulated in the parp-1 −/− group, and five were upregulated, ranging from − 33.96 for the mobile element ivk to 2505.47 for the mobile element springer. Five of these mobile elements are non-LTR retrotransposons belonging to the LINE-like-element superfamily 57 . ivk and I-element belong to Clade I, whereas TART-element, He T-A, and juan belong to the Clade jockey 57 . The three remaining MEs are LTR retrotransposons belonging to the gypsy superfamily 57 . diver1 belongs to the Clade roo, transpac belongs to the Clade 17.6, and springer belongs to the Clade gypsy 57 . We also found ten MEs that were upregulated in the parp-1 −/− group, ranging in fold difference from 2.32 to 35.32 (Table 4, second panel). For four, the p values comparing the significance of change of expression between the wild type and parp-1 −/− groups relative to the standard deviation of all measurements were below 0.05, but the FDRs were above 15%. For the other six MEs, p values were greater than 0.05. Among those ten MEs, nine were LTR retrotransposons belonging to the gypsy superfamily, and the last one was a non-LTR retrotransposon belonging to the LINElike-element superfamily 57 . These results suggest that several mobile elements are misregulated in the absence of PARP-1 and that the majority of them were upregulated in the parp-1 −/− group.
In addition, we found that the arc2 gene was upregulated in the parp-1 −/− group compared to the wild-type group (fold difference of 2.98) ( Table 4, third panel). Arc proteins are composed of Group-specific antigen (Gag)like amino acid sequences typical of retroviruses and retrotransposons 58 . Interestingly, we also found that arc1 was upregulated in the parp-1 −/− group. However, the difference was not statistically significant by having a p value greater than 0.05. PARP-1 might be important for repressing retrotransposon and retrotransposons-like gene expression.
Finally, we identified a number of genes involved in retrotransposon regulation among DEGs (Table 4, fourth panel). Among them is a gene coding for the transcription factor-like Modifier of mdg4 (Mod(mdg4)). Mod(mdg4) is a member of the gypsy insulator complex reported to repress the mobility of P-element transposon 59 . Consistent with the upregulation of most mobile elements, Mod(mdg4) was downregulated in the parp-1 −/− group. However, the other members of this complex (Su(hw), CP190, and Topors) were not misregulated in the absence of PARP-1 60 . Finally, the gene coding for the Lola transcription factor was downregulated in the absence of PARP-1. Lola has been reported to repress the expression of retrotransposons in CNS 61 . The downregulation of lola expression may be involved in the upregulation of retrotransposons in the absence of PARP-1. Taken together these findings suggest that mobile elements are mainly upregulated in the absence of PARP-1.

Discussion
Our results show that the absence of PARP-1 triggers large-scale changes in the expression of 602 genes in vivo, suggesting that PARP-1 mediates the expression profile of those genes. Among DEGs with known function, several gene-ontology terms were overrepresented, including 'Signal transduction, ' 'DNA-binding transcription factor activity, ' 'Ubiquitin-protein transferase activity, ' 'Nervous system development, ' 'Oxidation-reduction process, ' 'Heme binding, ' 'Chitin metabolic process, ' and 'Transmembrane transport' . We found little overlap between the gene-ontology terms overrepresented among the up-and downregulated DEGs. Therefore, it appears that PARP-1 stimulates the transcription of genes responsible for certain functions and inhibits the transcription of other functional gene groups. Thus, genes involved in 'DNA-binding transcription factor activity, ' 'Ubiquitin-protein transferase activity, ' 'Signal transduction, ' 'Nervous system development' and 'Intracellular protein transport' were overrepresented among the downregulated DEGs, whereas genes involved in 'Oxidation-reduction process, ' 'Heme binding, ' 'Chitin metabolic process, ' 'Carbohydrate metabolic process' and 'Transmembrane transport' were overrepresented among the upregulated DEGs (Fig. 4).
More than a half of DEGs with known functions were not involved in GO-terms overrepresented among DEGs, suggesting that many standalone genes are controlled by PARP-1, even though several major functional groups of genes with transcription profiles are regulated by PARP-1. Most genes regulated by PARP-1 do not www.nature.com/scientificreports/ appear to share their functions with one another. In those instances when a functional group of genes is regulated by PARP-1, all members of the group are regulated in a similar way (i.e., up-or downregulated) (Fig. 4). We also checked in the literature the phenotype of the knockdown/knockout of the ten most downregulated DEGs in parp-1 −/− group to see if they are compatible with the phenotype observed in parp-1 mutant. Seven of them were knocked down using RNAi. For six of them, the flies are viable and do not exhibit any phenotype: CG7900 (− 387.5), alpha-Est2 (− 283.2), Obp57b (− 80.6) 62 , CG11893 (− 282.5), CG11034 (− 132.5) and CG44014 (− 90) 63 . The last one knocked down with RNAi is lethal before pupation (IntS12 (− 83.7)) 63 , similar to parp-1 knockout phenotype. Two of them were knocked out using gene trap and are lethal before pupation, Glut1 (− 88.5) 64 and PGRP-LE (− 107.7) 65 . Finally, the flies knocked out for Lama (− 110) using enhancer trap are viable 66 . All these phenotype are compatible with parp-1 knockout phenotype.
The transcription level of all twelve transcription factors found among DEGs was strongly downregulated in the absence of PARP-1. Eight of these transcription factors are directly involved in neuronal development and axonal growth. This finding is consistent with a study showing that post-mitotic neurons expressing RNAi against parp-1 present defects in axonal outgrowth and branch patterning in vitro 67 . Taken together, these results suggest that PARP-1 stimulates the expression of transcription factors, which mediate neuronal development and axonal growth. However, we did not detect any misregulation in the expression of known downstream targets for these transcription factors in the absence of PARP-1. It is possible that our approach, which focuses on large-scale Table 4. The first panel of the table lists mobile elements (MEs) misregulated in the absence of PARP-1. The second panel represents MEs were misregulated in the absence of PARP-1, but presented an FDR > 15% or a p value > 0.05. The third panel lists genes that act like endogenous retroviruses. The fourth panel lists genes reported to play a role in the regulation of MEs. Green shading marks DEGs with a fold difference > 2.0, p values < 0.05, and FDR < 15%. Orange shading marks DEGs with fold difference > 2.0, p values < 0.05, and FDR > 15%. Red shading marks DEGs with fold difference > 2.0, p values > 0.05, and FDR > 15%. DEGs downregulated in the parp-1 −/− group are marked in yellow. The number of insertions reported in Drosophila melanogaster is based on Kaminker et al. 57  www.nature.com/scientificreports/ differences in the transcription profile, did not detect tissue-specific misregulation. Therefore, genes expressed in several tissues, but only misregulated in neurons or neuronal progenitors in the absence of functional PARP-1, did not show fold difference at transcription levels sufficient to be recognized as DEGs in our study. Another possibility is that the absence of misregulation of downstream targets in parp-1 −/− group may be due to the fact that many developmentally regulated transcription factors serve as a pioneer factors and downstream activation or repression requires the involvement of other factors that act later in the development 68 . Aside from one DEG, all other DEGs shown to be related to the cytochrome P450 and Cytochrome-b5 families were upregulated in the absence of PARP-1. Among the 21 cytochrome-related DEGs, 14 have been reported as contributing to the clearance of toxic chemicals, including DTT and permethrin, suggesting that PARP-1 activity may reduce the resistance to toxic chemicals by limiting the expression of genes involved in the clearance of these toxic chemicals. Roles of most members of the cytochrome family listed in Table 3 remain poorly understood since seven only have predicted functions based on their DNA sequences. Apart from their role in toxic chemicals clearance, it is possible that all these genes are also involved in other processes, such as hormone synthesis/clearance.
Several mobile elements were also misregulated in the absence of PARP-1. Most were among the upregulated DEGs, consistent with one of our earlier studies that showed a de-repression of retrotransposons in the absence of PARP-1 69 . We also found that two known retrotransposon repressors were downregulated in the absence of PARP-1: Lola 61 and Mod(mdg4) 59 . We then cannot exclude the possibility that the retrotransposon de-repression observed in parp-1 −/− group may be a consequence of the downregulation of the expression of retrotransposon repressors such as Lola and Mod(mdg4) rather than a direct upregulation of retrotransposon expression due to the absence of PARP-1.
Possible mechanisms of action for PARP-1 in regulating expression. The main function of PARP-1 is the poly(ADP)ribosylation (PARylation) of target proteins. It is a post-translational modification involving the polymerization of ADP-ribose units. This modification is highly negatively charged and can lead to repulsion between the target proteins and DNA 7 . Two alternative mechanisms of PARP-1 involvement in the regulation of gene expression level are possible: (1) direct involvement through the regulation of chromosome condensation to facilitate or repress access to promoter of target genes or (2) indirect involvement through the PARylation of transcription factors/cofactors/insulators of the target genes. These two mechanisms are not mutually exclusive.
PARP-1 is involved in chromatin loosening through its activation by JIL-1 9 . Such PARP-1-mediated chromatin loosening could lead to activation of DEGs that were downregulated in the absence of PARP-1. Alternatively, because PARP-1 can compact chromatin in vitro 70 and because we have demonstrated that PARP-1 is enriched during the interphase in regions where gene expression is silenced in vivo 71 , PARP-1 could potentially repress transcription by maintaining a compact chromatin state.
PARP-1 could also play an indirect role in regulating DEG expression through the PARylation of transcription factors/cofactors important for their correct expression 72 . It has also been reported that PARP-1 is involved in the PARylation of insulators that lead to their repulsion from DNA 73 . The insulator is then unable to block the enhancer/promoter interaction with its target gene, leading to the activation of the repressed gene.

Materials and methods
Drosophila melanogaster strains used. Flies were cultured on standard cornmeal-molasses-agar media at 22 °C, unless otherwise indicated. All the fly stocks listed below are isogenic. The fly stocks were generated by the standard genetic methods or obtained from the Bloomington Drosophila Stock Center and the Exelixis Collection at the Harvard Medical School. Genetic markers are described in Flybase 17 . The parp-1 C03256 strain was generated in a single pBac-element mutagenesis screen 74 . Precise excision of parp-1 C03256 was carried out using pBac transposase on CyO chromosome. Balancer chromosomes carrying Kr::GFP, i.e., TM3, Sb, P{w + , www.nature.com/scientificreports/ Kr-GFP4} and FM7i, P{w1, Kr-GFP} 75 , were used to identify heterozygous and homozygous parp-1 C03256 . The genetic background is similar in both parp-1 C03256 and yellow white strain. Larvae from both genotypes were synchronized by placing adult flies in a fresh vial and letting them lay eggs for three hours before transferring them to another vial. 50 third instar larvae were collected for each biological replicate for both conditions without sex selection.
RNA isolation and microarray. Total RNA was purified using the RNeasy Mini kit (Qiagen, Valencia, CA) after its isolation using TRIzol reagent (Life Technologies, Inc., Grand Island, NY). Microarray services, including quality control tests of the total RNA samples by Agilent Bioanalyzer and Nanodrop spectrophotometry, were provided by the UPENN Molecular Profiling Facility. All protocols followed the NuGEN Ovation User Guide and the Affymetrix GeneChip Expression Analysis Technical Manual. Briefly, 100 ng of total RNA were converted to first-strand cDNA using the reverse transcriptase primed by poly(T) and random oligomers that incorporated an RNA priming region. Second-strand cDNA synthesis was followed by ribo-SPIA linear amplification of each transcript using an isothermal reaction with RNase, RNA primer, and DNA polymerase (Ovation RNA Amplification System V2, NuGEN Inc., San Carlos CA, USA), and the resulting cDNA was assessed by Bioanalyzer, fragmented and biotinylated (Encore Biotin Module, NuGEN Inc., San Carlos CA); 3.75 μg of labeled cDNA were added to Affymetrix hybridization cocktails. Target hybridization was performed on Gene-Chip Drosophila Genome 2.0 arrays (Affymetrix Inc., Santa Clara CA, USA), according to the manufacturer's procedures found in GeneChip Hybridization Oven 645, followed by washing and staining in GeneChip Fluidics Station 450. Data were acquired with the GeneChip Scanner 3000 7G (Affymetrix, Inc., Santa Clara, CA, USA). Three biological replicates per group were made (wild type and parp-1 knockout).
Microarray data analysis. Data analysis was performed using Partek Genomics Suite v6.5 to apply the GCRMA normalization algorithm. Differential expression analysis was conducted using Significance Analysis of Microarrays (SAM, v3.09) and 2-way ANOVA tool. SAM calculates a score for each gene based on the change in expression relative to the standard deviation of all measurements by computing a t-test based on the three biological replicates. SAM then performs a set of permutations to determine the false discovery rate (FDR) with an adjustment for multiple testing. The reported list of ranked genes has a 'delta value' which defines the threshold of false positive in the validated dataset, which was adjusted to a stringent FDR < 15% 76 . Affymetrix .cel files were imported into the Partek Genomic Suite (v6.5, Partek Inc., St. Louis, MO), and GCRMA normalization was applied. Data were exported to SAM (v3.09, Significance Analysis of Microarrays, Stanford University) to test for differential expression, yielding a fold change, a p value based on a t-test realized on the three biological replicates, and an FDR. The fluorescence intensity of each probe was recorded as a log 2 value. To determine the expression level for every gene, the formula 2 X was used, where X is the average fluorescence intensity of the three biological replicates. The ratio between the expression level in the wild type and the parp-1 knockout group, referred to as fold difference, was used to evaluate the changes of expression levels between the wild type and parp-1 knockout groups for every gene. A negative value corresponds to a gene that is upregulated in wild-Scientific Reports | (2020) 10:21151 | https://doi.org/10.1038/s41598-020-78116-5 www.nature.com/scientificreports/ type larvae compared to parp-1 knockout larvae, while a positive value corresponds to a gene upregulated in parp-1 knockout larvae. 602 genes were found to be differentially expressed using cutoffs of twofold increase or decrease with a False Discovery Rate < 15% and a p value for the two-tailed t-test < 0.05. Using a combination of p value and FDR allows us to identify Differentially Expressed Genes (DEGs), while minimizing the number of false positives 76 . The fold difference establishes an arbitrary threshold for a difference of expression level between the wild type and parp-1 knockout groups that we accepted as sufficient difference to consider a gene as misregulated between the two groups.
DEGs cellular functions and GO-terms overrepresentation analysis. The main cellular functions of the genes differentially expressed between the wild-type and parp-1 knockout groups was determined using PANTHER classification system 77,78 .
Overrepresentation of Gene Ontology-terms (GO-terms) among DEGs was determined using Gene Set Enrichment Analysis (GSEA) software, which has been described in Subramanian et al. and Mootha et al. 79,80 (see Supplementary Fig. 1 for details). The metric used to rank genes has been set up on "Signal2Noise". GOterms with fewer than 4 DEGs were excluded during the analysis. Only the GO-terms with p value < 0.05 and FDR < 25% were considered as overrepresented in one of our groups. To avoid redundancy, when a GO-term was included in another GO-term, only the GO-term representing the smallest subset of genes was selected.