Age-related ceRNA networks in adult Drosophila ageing

As Drosophila is an extensively used genetic model system, understanding of its regulatory networks has great significance in revealing the genetic mechanisms of ageing and human diseases. Competing endogenous RNA (ceRNA)-mediated regulation is an important mechanism by which circular RNAs (circRNAs) and long non-coding RNAs (lncRNAs) regulate ageing and age-related diseases. However, extensive analyses of the multiomics (circRNA/miRNA/mRNA and lncRNA/miRNA/mRNA) characteristics of adult Drosophila during ageing have not been reported. Here, differentially expressed circRNAs and microRNAs (miRNAs) between 7 and 42-day-old flies were screened and identified. Then, the differentially expressed mRNAs, circRNAs, miRNAs, and lncRNAs between the 7- and 42-day old flies were analysed to identify age-related circRNA/miRNA/mRNA and lncRNA/miRNA/mRNA networks in ageing Drosophila. Several key ceRNA networks were identified, such as the dme_circ_0009500/dme_miR-289-5p/CG31064, dme_circ_0009500/dme_miR-289-5p/frizzled, dme_circ_0009500/dme_miR-985-3p/Abl, and XLOC_027736/dme_miR-985-3p/Abl XLOC_189909/dme_miR-985-3p/Abl networks. Furthermore, real-time quantitative PCR (qPCR) was used to verify the expression level of those genes. Those results suggest that the discovery of these ceRNA networks in ageing adult Drosophila provide new information for research on human ageing and age-related diseases.


Introduction
Drosophila melanogaster is an extensively used genetic model system that has been used for more than 100 years to study various aspects of the life sciences. In particular, fruit flies have been widely utilized to study ageing (Gubina et al., 2019) and human diseases, such as cancer (Enomoto et al., 2018), neurodegenerative disease (Cha et al., 2019), obesity and diabetes (Musselman et al., 2019), sterile inflammation (Nainu et al., 2019), and regeneration (Fox et al., 2020). D. melanogaster has also been utilized to study complex behavioural and developmental biology topics, including exercise (Watanabe and Riddle, 2019), courtship , and foraging (Khodaei and Long, 2019). Recently, mounting evidence has suggested that Drosophila is an outstanding model for studying ageing and age-related diseases (Surguchov et al., 2019;Brenman-Suttner et al., 2020). Ageing is a physiologic/ pathologic process featuring declines in normal physiological functions and progressive impairment of cellular functions (Stallone et al., 2019). The ageing phenomenon has been conserved during biological evolution; even yeast other single-celled eukaryotes experience ageing (He et al., 2019). Thus, in-depth study of the regulatory mechanism of ageing in Drosophila can inform the study of human ageing and disease.
Current studies suggest that non-coding RNAs (ncRNAs) are involved in organismal ageing Tower, 2019). With the development of sequencing technology, dynamic changes in the transcriptome [including changes in mRNA, long non-coding RNA (lncRNA), microRNA (miRNA), and circular RNA (circRNA) (Yang et al., 2016;Perry et al., 2017;Barter et al., 2019;Kinser and Pincus, 2020); the proteome ; and the metabolome (Song et al., 2017)] have been described in the context of Drosophila ageing. Most miRNAs, mRNAs, and the proteins in fruit flies are evolutionarily conserved up to humans and regulate similar signalling pathways across organisms, such as the NF-κB, AMPK, mTOR, P53, PGC1α, and FoxO pathways (Zha et al., 2019;Kinser and Pincus, 2020;. Large studies have demonstrated that the circRNA/miRNA/mRNA and axis the lncRNA/miRNA/mRNA axis play vital roles in ageing and age-related disease (Ruan et al., 2020;. As Drosophila is a workhorse model organism, thoroughly studying the characteristics of adult Drosophila from a multiomics perspective is necessary. Understanding how these conserved protein genes regulate ageing through ceRNA mechanisms in Drosophila is important. However, extensive analyses of the multiomics (circRNA/miRNA/mRNA lncRNA/ miRNA/mRNA) characteristics of adult Drosophila ageing have not been reported. In the present study, we investigated the circRNA/miRNA/mRNA the lncRNA/miRNA/mRNA axis in adult Drosophila at two age points (day 7 and day 42) and determined the regulatory network of key differentially expressed (DE) genes in Drosophila ageing. The results provide knowledge on the gene regulation network of adult Drosophila ageing and a solid foundation for understanding the mechanisms of human ageing and age-related diseases.

Overview of multiomics data
The DE genes and proteins in Drosophila between day 7 and day 42 were identified, and their networks were analysed (Table 1). A total of 537 DE mRNAs and 43 DE lncRNAs were obtained from a previous study in our laboratory (Supplementary Data Sheet S1). A total of 6,003 circRNAs and 226 miRNAs were identified at day 7 and day 42 (Supplementary Data Sheet S1). Ultimately, 29 DE circRNAs and 24 DE miRNAs were found (Supplementary Data Sheet S1). The merged sequences of novel circRNAs (Supplementary Data Presentation S1-S3) and lncRNAs (Supplementary Data Sheet S3) are shown in supplementary files.

DE circRNAs and miRNAs in Drosophila between day 7 day 42
The DE circRNAs and miRNAs between 7 and 42-day old flies were analysed. Between day 7 and day 42, 29 DE circRNAs in Drosophila were identified, including 21 upregulated and 8 downregulated circRNAs at day 42 ( Table 2). The circRNAs were derived from different source genes. These source genes were found to be involved in multiple molecular functions (Table 2). Evidently, the biological processes of the short lifespan-related source genes arm and pan were involved in the Wnt signalling pathway and had similar molecular functions, such as binding, protein binding, and transcription factor binding, transcription regulator activity. In addition, different circRNAs were observed to originate from the same mRNA transcript. For example, and dme_circ_0008175 and dme_circ_0008173 originated from the Nlg1 gene, and dme_circ_0009514 dme_circ_0009500 were derived from the pan gene.
Furthermore, 24 DE miRNAs (11 upregulated and 13 downregulated) were identified at day 42 compared to day 7 (Supplementary Data Sheet S1). Dme_miR-9a-3p and dme_miR-985-3p were identified as canonical specific fruit fly miRNAs. Then, dme_miR-956-3p, dme_miR-284-3p, and dme_miR-289-5p were identified as non-canonical miRNAs. The remaining 19 DE miRNAs were canonically conserved miRNAs. Then, functional annotation was carried out for the DE miRNAs. The Gene Ontology (GO) annotations of age-related DE miRNAs based on their targets were investigated. 11 upregulated and 13 downregulated miRNAs targeted to 3,292 mRNAs and 2,951 mRNAs, respectively (Supplementary Data Sheet S1), which mainly enriched in biological process (1,314 terms and 1,233 terms, respectively) (Supplementary Data Sheet S1). The bar plot shows the top ten enrichment score value of the significant enrichment terms ( Figures  1A, B), such as multicellular organism development, nervous system development, neurogenesis, development process, and cell differentiation. Specifically, 20 DE miRNAs were related to Drosophila aging based on previous reports, including 11 upregulated and 9 downregulated miRNAs in 42 days when compared to 7 days ( Figure 1C). 13 DE miRNAs involved in the age-signalling pathways by target genes on post-transcriptional level ( Figure 1C). Their tagets involved in regulation of ROS detoxification, autopage, circadian rhythm, apoptosis, and immunity biological process. Otherwise, three out of 24 total miRNAs were conserved in Drosophila, human and mouse, containing dme-miR-8-5p, dme-miR-133-3p, and dme-miR-10-5p.

Functional annotation of DE circRNA/ lncRNA-associated networks
GO functional annotation of DE circRNAs/lncRNAs/mRNAs was carried out based on 74 target mRNA genes (Supplementary Data Sheet S9). The first 30 GO terms based on the lowest p values are listed ( Figure 4). In circRNA/mRNA GO terms, there were just two major GO categories in the circRNA-associated networks, including biological processes (29 GO terms) and cellular components (1 GO term, perinuclear region of cytoplasm). Similarly, the GO annotations of the lncRNA-associated networks consisted of 25 GO terms in biological processes and 5 GO terms in the cellular component.

Tissue expression pattern analysis
The tissue-specific expression patterns of the DE lncRNAs, circRNAs, miRNAs, and mRNAs were analysed in the head, ovary, gut, and fat body, the qPCR results of which were consistent with the RNA-seq data in the DE circRNA/miRNA/ mRNA networks and lncRNA/miRNA/mRNA networks ( Figure 5). The results showed that the tissue-specific expression patterns of dme_circ_0009500/dme_miR-289-5p/CG31064, dme_circ_ 0009500/dme_miR-289-5p/frizzled, and dme_circ_0009500/dme_ miR-985-3p/Abl were mainly expressed in the head. Specifically, FIGURE 2 DE circRNA/miRNA/mRNA networks of Drosophila at day 7 day 42. (A), circRNA/miRNA/mRNA interaction networks; yellow, circRNAs; red, mRNAs; green, miRNAs. (B), Specific circRNA/miRNA/mRNA networks; "↑" "↓" indicate upregulation downregulation on day 42 compared to day 7 in RNA-seq data, respectively. (C), Data on the relative expression of circRNA, miRNA, mRNA detected by qPCR analysis. The relative gene expression levels were calculated by the cycle threshold values that were identified as 2 −ΔΔCT . The ribosomal protein L32 (rp 49) gene was used as the reference gene to calculate the relative mRNA, miRNA, circRNA levels. "*" above the bars indicates a significant difference at the 0.05 level, "**" indicates a significant difference at the 0.01 level. (+) represents the qPCR results were consistent with the RNA-seq data.
Frontiers in Genetics frontiersin.org dme_circ_0009667 was mainly located in the ovary, and dme_miR-14-5p was mainly expressed in the gut. Similar results were found in the lncRNA/miRNA/mRNA networks. The tissue patterns of the XLOC_027736/dme_miR-985-3p/Abl and XLOC_189909/dme_ miR-985-3p/Abl networks were also mainly expressed in the head.

Discussion
Drosophila is an ideal model for genetics, and the multiple agerelated researches were carried out by Dahomey strain (Mołoń et al., 2020;De Groef et al., 2021). Previous study has reported the trend of wild-type female lifespan (Dahomey, Canton S, Oregon R) changes were similar (Sanz et al., 2010). Furthermore, several studies just used one wildtype strain in multiple-omics research (Shi et al., 2020;Wang et al., 2022). Thus, the wild-type female of Dahomey strain was utilized in our manuscript. Increasing evidence has suggested FIGURE 3 DE lncRNA/miRNA/mRNA networks in Drosophila at day 7 day 42. (A), lncRNA/miRNA/mRNA interaction networks; yellow, lncRNAs; red, mRNAs; green, miRNAs. (B), Specific lncRNA/miRNA/mRNA networks; "↑" "↓" indicate upregulation downregulation on day 42 compared to day 7, respectively. (C), Data on the relative expression of lncRNAs, miRNAs, mRNAs by qPCR analysis. 7 days, 7 days; 42 days, 42 days. The relative gene expression levels were calculated using the cycle threshold values that were identified as 2 −ΔΔCT . The ribosomal protein L32 (rp 49) gene was used as the reference gene to calculate the relative mRNA, miRNA, lncRNA levels. "*" above the bars indicates a significant difference at the 0.05 level, "**" indicates a significant difference at the 0.01 level. (+) represents the qPCR results were consistent with the RNA-seq data.
Frontiers in Genetics frontiersin.org that ceRNA networks play key roles in a variety of biological processes, such as cancer (Wang et al., 2018;Abdollahzadeh et al., 2019), Alzheimer's disease (AD) , skeletal muscle myogenesis (Yue et al., 2019), and ageing Chen et al., 2020). miRNAs are the core molecules of the ceRNA regulatory system (Yue et al., 2019;Chen et al., 2020). Thus, screening DE miRNAs is essential for research on ceRNA networks related to Drosophila ageing. miRNAs generally induce mRNA degradation or repress translation of target transcripts through sequence-specific binding to the transcript 3′UTR (Bushati and Cohen, 2007;Chen et al., 2019). Each transcript can be targeted by multiple miRNAs, and each miRNA can target hundreds of different transcripts (mRNA, circRNA, and lncRNA transcripts) (Dori and Bicciato, 2019;Zhou et al., 2019;Kinser and Pincus, 2020). Thus, the miRNA regulatory network is far-reaching (Kinser and Pincus, 2020). Previous studies have verified that miRNAs are important small regulatory ncRNA molecules that control a fairly large number of biological processes; their important functions have generated interest in their use as biomarkers and their roles as regulators of ageing (Lai et al., 2019;Kinser and Pincus, 2020) and a number of cancer types (Liang et al., 2020;Mishan et al., 2020;. Recent studies have also suggested that miRNAs are involved in the regulation of age-associated processes and pathologies in multiple mammalian tissues, including the brain, heart, bones, and muscles John et al., 2020;Kinser and Pincus, 2020;Ullah et al., 2020). In our study, more than 80% (20 out of 24 total miRNAs) of DE miRNAs between day 7 and day 42 could affect fruit fly ageing. Among these miRNAs, miR-14 as a cell death suppressor, regulates fat metabolism, insulin production and metabolism through its targets (Xu et al., 2003;Nelson et al., 2014). In addition, knockout of five DE miRNAs (miR-133, miR-284, miR-286, miR-318, miR-956, and miR-988) decreased the Drosophila lifespan, and miR-286 KO increased female lifespan (Chen et al., 2014). The consreved miR-8 was the homolog of vertebrate miR-200 family. It is worth noting that miR-8 acted through U-shaped to activate PI3K and thereby promoted fat cell growth cell-autonomously and the Insulin-like Receptor signaling pathway (Soler Beatty et al., 2021;Chen et al., 2022). Then, two targets of miR-200 (dme-miR-8 homolog) in human were JAGGED1 (JAG1) (Xue et al., 2021) and epidermal growth factor receptor (EGFR) (Xue et al., 2019). In a lung cancer metastasis, overexpression of the miR-200 decreased JAG1 protein levels and impeded cell growth (Xue et al., 2021). Furthermore, Dmel\ EGFR is orthologous to human gene ERBB2, Frontiers in Genetics frontiersin.org which has also been implicated in multiple cancers (Xue et al., 2019).
The results indicate that these miRNAs may act as sponges through the circRNA/miRNA/mRNA and lncRNA/miRNA/mRNA networks involved in the regulation of Drosophila ageing. In a previous study, global identification of functional miRNA-mRNA interactions in Drosophila was performed (Wessels et al., 2019). However, circRNA/miRNA/mRNA and lncRNA/miRNA/mRNA interactions in fly ageing have remained unclear.
Frontiers in Genetics frontiersin.org Furthermore, the lncRNAs XLOC_027736 and XLOC_189909, the dme_circ_0009500 and the mRNA Abl were predicted to competitively bind miR-985-3p, which was mainly expressed in the Drosophila head. Abl phosphorylates cell adhesion and cytoskeletal proteins and acts as a scaffold in a signalling complex to regulate both epithelial and nervous system morphogenesis (Zhu and Bhat, 2011;Liu and Wu, 2014). The age-related phenotype associated with the Abl mutant led to a shorter lifespan for flies with three specific alleles (Abl l2 , Abl l3 , and Abl GD1344 ) (Belote et al., 1990;Huang et al., 2007;Neely et al., 2010). Furthermore, human ABL1 (Drosophila Abl homolog) protein kinases play many important roles in neuron development, maintenance, and signalling (Manley et al., 2022). In future studies, it will be necessary to determine the biological role of the lncRNAs XLOC_027736 XLOC_189909, the circRNA dme_ circ_0009500, and miR-985-3p in the fly ageing process.
In our study, DESeq2 and EdgeR were used to analyse the DE circRNAs and miRNAs, respectively. DESeq2 and EdgeR are efficient tools for differential analysis of RNA-seq data with the more than 80% overlapping, both of which use the negative binomial distribution (Liu et al., 2021). Then, DESeq2 can more accurately identify DE genes for small samples to reduce false positives (Love et al., 2014). Then, two biological replicates were utilized to analysis DE circRNAs. Thus, DESeq2 was chosen to identify the DE circRNAs. The above results suggest that our data are also available. To date, DE circRNAs and miRNAs cannot be reanalysed only by DESeq2 or EdgeR in the present study. This would bring some questions to our subsequent analysis. For example, some DE circRNAs and miRNAs may be lost, which may result that part ceRNA networks could not be recognized. More experiments will be needed to verify these networks in our future study.
In this study, several DE miRNAs were identified between day 7 and day 42, such as the dme_miR-289, dme_miRNA-14, and the conserved miR-8. Then, the potential ceRNA networks may play a role in Drosophila aging, for example, dme_circ_0009500/dme_miR-289-5p/frizzled and XLOC_027736/dme_miR-985-3p/Abl networks. Therefore, the result of DE miRNAs, circRNA/miRNA/mRNA, and lncRNA/miRNA/mRNA interactions provides an important foundation to parse the genetic process of Drosophila ageing.

Sample collection and preparation
Female flies (Dahomey WT ) that had mated with male flies for 48 h after hatching were bred under a 12 h on/off light cycle at 25°C in 50% humidity. Sample collection of adult flies at day 7 and day 42 was performed as described by Yang et al. (2016). RNA-seq (mRNA, lncRNA, and circRNA) was conducted on two biological replicates, while miRNA sequencing was conducted on five biological replicates. All samples were stored at −80°C until use.
4.2 LncRNA, mRNA, circRNA data from Drosophila at days 7 and 42 RNA-seq data for the lncRNAs and mRNAs of Drosophila at day 7 and day 42 were obtained by Yang et al. (2016). The circRNA analysis was based on the RNA-seq data for Drosophila at day 7 and day 42 from Yang et al. (2016). Overall, 95.79% clean reads were obtained from 26.7 GB of raw sequence data (SRP073695) then aligned to the D. melanogaster genome from FlyBase (Dmel_ Release_6, http://FlyBase.org/). The find_circ (Memczak et al., 2013) and CIRI2 (Gao et al., 2015) software tools were utilized to identify circRNAs. Then, the overlapping Drosophila circRNAs identified by both software programs were selected. The input data for the circRNA differential expression analysis were readCount data obtained from the circRNA expression level analysis. Then, paired differential expression analysis of circRNAs between day 7 and day 42 was conducted with DESeq2 (Varet et al., 2016) based on a negative binomial distribution. The p-value was adjusted using Hochberg and Benjamini's methods (Hinkley et al., 2018) to control the error discovery rate. A p-value <0.05 was considered to indicate a DE circRNA. These original data were from our laboratory.

MiRNAs in Drosophila at day 7 and day 42
TRIzol Reagent (Invitrogen, CA) was used to extract the total RNA from fruit flies at day 7 and day 42. Agarose gel electrophoresis was performed to verify the integrity of the total RNA samples. A NanoDrop ND-1000 instrument was used to accurately measure the concentrations and protein contamination of the total RNA samples. miRNA sequencing libraries were generated using rRNA-depleted RNA with a NEBNext ® Ultra ™ Multiplex Small RNA Library Prep Set Kit for Illumina ® (New EnglBiolabs, United States) following the manufacturer's recommendations. Subsequently, an Agilent 2,100 Bioanalyser and an Agilent DNA 1000 chip kit (Agilent, part #5067-1504) were utilized to accurately assess the quality and concentration of the sequencing libraries. The libraries were sequenced using an Illumina NextSeq 500. MiRNA fragment sequencing was performed by the Aksomics company.
Clean reads were generated from the raw sequence data from the Illumina NextSeq instrument through real-time base calling and quality filtering. The clean reads were recorded in FASTQ format and contained read information, sequences and quality encoding. Subsequently, the 5′-and 3′-adapter sequences were trimmed from the clean reads by Cutadapt, and reads with lengths shorter than 14 nt or longer than 40 nt were discarded. The trimmed reads were collapsed into FASTA format. The raw data has been uploaded to NCBI database (PRJNA716466). The trimmed reads that did not map to mature or precursor tRNA sequences were aligned with an allowance of only one mismatch to miRNA reference sequences with miRDeep2 (Friedlander et al., 2012). The expression profiles of miRNAs were determined based on the counts of the reads mapped. The DE miRNAs were identified with the R package EdgeR based on the count values (Robinson et al., 2010). A fold change cut-off of 1.5 and a p-value cut-off of 0.05 were applied only when replicates were used for screening DE miRNAs. The gene prediction of DE miRNA integrates two algorithms, miRanda (Enright et al., 2003) and TargetScan (Garcia et al., 2011). GO enrichment analysis of the targets of DE miRNAs was implemented with the GOseq R package (Young et al., 2010). The lncRNAs, circRNAs, and miRNAs showed significantly different expression levels between day 7 and day 42 and were thus analysed. The potential ceRNAs were searched based on the sequences of the lncRNAs, circRNAs, and mRNAs. The offline software MiRanda (Enright et al., 2003) was utilized to predict miRNA binding seed sequence sites, and overlap of the same miRNA binding sites on lncRNAs/circRNAs-miRNAs and miRNAs-mRNAs was taken to indicate a lncRNA/circRNA-miRNA-mRNA interaction. Then, the "clusterProfiler" R package was utilized to perform Gene Ontology (GO) enrichment of ceRNA networks based on the mRNAs.
The DE mRNAs associated with ageing and the corresponding ceRNA networks (including the circRNAs, lncRNAs, miRNAs, and mRNAs) were selected to detect the expression level by qPCR. In total, 9 miRNAs, 10 mRNAs, 5 lncRNAs, and 5 circRNAs were selected to qualify the expression level. Samples (three biological duplicates) from 7-and 42-day-old fruit flies were used to isolate the total RNA using TRIzol ™ LS reagent (Thermo Fisher) according to the manufacturer's instructions. The total RNA (1 μg) from each sample was reverse transcribed with random primers using a RevertAid First-strand cDNA Synthesis Kit (Thermo Fisher) according to the manufacturer's protocol, which was utilized to detect the expression levels of mRNA, lncRNA, and circRNA. In addition, the total RNA (1 μg) used for miRNA expression detection in each sample was reverse transcribed using a TaqMan ™ MicroRNA Reverse Transcription Kit (Thermo Fisher) according to the kit instructions. All gene primers (Supplementary Data Sheet S10) were designed using the Primer 5.0 software and purchased from Sangon Biotech. The SYBR green method was used for qRT-PCR with TransStart ® Green qPCR SuperMix (TransGen Biotech) following the manufacturer's instructions. All tested samples were three biological replicates. The relative gene expression levels were calculated using cycle threshold values and the 2-ΔΔCT method. Ribosomal protein L32 (rp 49) was used as the reference gene to calculate the relative mRNA, miRNA, lncRNA, and circRNA levels. Differential expression levels were compared by independentsamples t-tests between groups.

Tissue-specific expression pattern
Based on the same expression pattern between day 7-and day 42-day-old flies, as determined by the RNA-seq data and qPCR tests on 7-day and 42-day-old flies, the tissue-specific expression model of four circRNAs (dme_circ_0008173, dme_circ_000950, dme_circ_0009667, and dme_circ_0010536), four miRNAs (dme_miR-289P, dme_miR-985-3P, dme_miR-286-3P, and dme_miR-14-5P), four mRNA (SERCA, frizzled, Abl, and CG31064), and two lncRNAs (XLOC_027736 and XLOC_ 189909) were selected to analyse the tissue specificity of the head, ovary, gut and fat body of female fruit flies at 7 days. Each sample included three biological replicates. The expression levels of these genes were detected by qRT-PCR as described above. Differential expression levels were compared by independent-samples t-tests between groups.

Data availability statement
The original contributions presented in the study are publicly available. This data can be found here: https://www.ncbi.nlm.nih. gov/ (Accession number PRJNA716466).
Author contributions DY and FX: writing data analysis; JL: data analysis; XF: quality control data verification; QN, YL, MZ, and TY: data visualization and quality control; DY: revision and typesetting; ZH and MY: design of the study revision.

Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.