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.

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.

In this section, the DE circRNA/miRNA/mRNA networks are analysed. Through the ceRNA mechanism, miRNAs can negatively regulate mRNA expression. Overall, 12 DE circRNAs, 21 DE miRNAs, and 30 DE mRNAs had interactions (Figure 2A, Supplementary Data Sheet S4). The binding sites between circRNA vs. miRNA (Supplementary Data Sheet S5) and miRNA vs. mRNA (Supplementary Data Sheet S6) are shown. According to the trends in the expression quantity changes, five DE circRNAs (dme_circ_0006913, dme_circ_0008173, dme_circ_0009500, dme_circ_0009667, and dme_circ_0010536) targeted the DE miRNAs with opposite expression trends (Figure 2B). Based on qPCR results, the expression patterns of four circRNAs (dme_circ_0008173, dme_circ_0009500, dme_circ_0009667, and dme_circ_0010536), and four miRNAs (dme_miR-289-5p, dme_miR-985-3p, dme_miR-286-3p, dme_miR-14-5p), four mRNAs (frizzled, CG31064, Abl and SERCA) were consistent with the RNA-seq data (Figure 2C). Importantly, the expression trends 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 conformed to the ceRNA mechanism.

According to the functional patterns of lncRNAs competing with mRNAs for binding to miRNAs, the interaction networks of DE lncRNAs/miRNAs/mRNAs were identified (Supplementary Data Sheet S7). Then, the binding sites between lncRNAs and miRNAs were included in Supplementary Data Sheet S8. In addition, DE miRNAs targeted DE mRNAs with opposite expression trends. Based on the DE genes in our database, 15 lncRNAs, 15 miRNAs, and 32 mRNAs had interactions (Figure 3A). Based on qPCR results, the expression patterns of two lncRNAs (XLOC_027736 and XLOC_189909), three miRNAs (dme_miR-289-5p, dme_miR-985-3p, and dme_miR-14-5p), and three mRNAs (frizzled, CG31064, and Abl) were consistent with the RNA-seq data (Figure 3A). Several specific lncRNA/miRNA/mRNA networks, including XLOC_027736/dme_miR-985-3p/Abl, XLOC_073604/dme_miR-994-3p-3p/mbl, XLOC_189909/dme_miR-985-3p/Abl (Figure 3B), were found. In the XLOC_027736/dme_miR-985-3p/Abl and XLOC_189909/dme_miR-985-3p/Abl networks, the expression trend of these genes was determined by qPCR analysis to be consistent with the RNA-seq data from our study (Figure 3C).

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.

In the GO annotations of the circRNA-associated networks, there were several ageing-related biological processes in GO terms, such as homeostatic process (10 DE genes), cellular homeostasis (7 DE genes), cation homeostasis (5 DE genes), ion homeostasis (5 DE genes), cell-cell signalling (9 DE genes), cell fate determination (5 DE genes), developmental growth (7 DE genes) (Figure 5A). Furthermore, 13 out of the first 30 GO terms of the lncRNA-associated networks were related to homeostatic processs, including cellular homeostasis, retinal homeostasis, ion homeostasis, calcium ion homeostasis, and cellular cation homeostasis (Figure 5B). In the ageing-related ceRNA networks dme_circ_0009500/dme_miR-985-3p/Abl, XLOC_027736/dme_miR-985-3p/Abl and XLOC_189909/dme_miR-985-3p/Abl in Drosophila, Abl was involved in multiple GO terms, such as regulation of cell morphogenesis, developmental growth, regulation of cell differentiation, and regulation of neuron differentiation. Furthermore, frizzled was involved in the positive regulation of developmental growth, developmental growth involved in morphogenesis, cell-cell signalling, cell fate determination, and the Wnt signalling pathway in the dme_circ_0009500/dme_miR-289-5p/frizzled network.

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, 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.

The binding sites of dme_circ_0009500/dme_miR-289-5p/CG31064, dme_circ_0009500/dme_miR-289-5p/frizzled, dme_circ_0009500/dme_miR-985-3p/Abl, XLOC_027736/dme_miR-985-3p/Abl, and XLOC_189909/dme_miR-985-3p/Abl were analyzed (Figure 6). Specificly, there were five binding sites between dme_circ_0009500 and dme_miR-289-5p with the higher binding free energy from -15.23 to -20.30 (Figure 6A; and Supplementary Data Sheet S5). Furthermore, 3′UTR sequence of frizzled had the two binding site with miR-289-5p, and also had the higher binding free energy -23.31 (Figure 6E and Supplementary Data Sheet S6).

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.

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.

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 independent-samples t-tests between groups.

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.