Age-related ceRNA networks and the mRNA/protein correlation in Drosophila ageing

As Drosophila is a classic model organism, understanding of the regulatory networks has great signicance in revealing the genetic mechanisms of ageing and human diseases. Competing endogenous RNA (ceRNA)-mediated regulation is an important mechanism of circular RNA (circRNA) and long non-coding RNA (lncRNA) regulation in ageing and age-related disease. However, extensive analyses of the multi-omics (circRNA/miRNA/mRNA, lncRNA/miRNA/mRNA, and mRNA/protein) characteristics of adult Drosophila ageing have not been reported. Our results determine the ceRNA networks and the mRNA/protein correlations of key differentially expressed genes in adult Drosophila aging, and provide a solid foundation for understanding mechanisms of human ageing and ageing-related diseases.


Conclusions
The discovery of these ceRNA networks and the mRNA/protein correlations in adult Drosophila ageing will provide new information for research on human ageing and age-related disease. Background Drosophila melanogaster is a classic and modern genetic model that has been used for more than 100 years to study various aspects of the life sciences. In particular, fruit ies have been widely utilized to study ageing [1] and multiple human diseases, such as cancer [2], neurodegenerative disease [3], obesity and diabetes [4], sterile in ammation [5], and regeneration [6]. D. melanogaster has also been utilized to study complex behavioural and developmental biology topics, including exercise [7], courtship [8], and foraging [9]. Recently, mounting evidence has suggested that Drosophila is an outstanding model for studying ageing and age-related disease [10][11][12]. Ageing is a physiologic/pathologic process featuring declines in normal physiological function and progressive impairment of cellular functions [13,14]. The ageing phenomenon has been conserved during biological evolution; even yeast and other single-celled eukaryotes experience ageing [13]. Thus, in-depth study of the regulatory mechanism of ageing in Drosophila will be helpful for the study of human ageing and disease.
As Drosophila is a workhorse model organism, thoroughly studying the characteristics of adult Drosophila from a multiomics perspective is necessary. However, extensive analyses of the multiomics (circRNA/miRNA/mRNA, lncRNA/miRNA/mRNA, and mRNA/protein) characteristics of adult Drosophila ageing have not been reported. In the present study, we investigated the circRNA/miRNA/mRNA axis, lncRNA/miRNA/mRNA axis, and mRNA/protein 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 will improve knowledge regarding the gene regulation network of adult Drosophila ageing and provide a solid foundation for understanding mechanisms of human ageing and ageing-related diseases.

Overview of multiomics data
The DE genes and proteins in Drosophila between day 7 and day 42 were screened, and their networks were analysed (Table 1). A total of 537 DE mRNAs, 348 DE proteins, and 43 DE lncRNAs were obtained from a previous study in our laboratory (Dataset S1). A total of 6003 circRNAs and 226 miRNAs were identi ed at day 7 and day 42 (Dataset S2). Ultimately, 29 DE circRNAs and 30 DE miRNAs were found (Dataset S1). The merged sequences of novel circRNAs (Dataset S3, Dataset S4, and Dataset S5) and lncRNAs (Dataset S6) are shown in supplementary les. The DE circRNAs and miRNAs between 7 and 42-day-old Drosophila were analysed. Between day 7 and day 42, 29 DE circRNAs in Drosophila were identi ed, including 21 upregulated and 8 downregulated circRNAs at day 42 ( Table 2). The circRNAs were derived from different source genes. Furthermore, 24 DE miRNAs (11 upregulated and 13 downregulated) were screened between day 7 and day 42 ( Figure 1A). Dme-miR-9a-3p and dme-miR-985-3p were identi ed as canonical speci c fruit y miRNAs. Then, dme-miR-956-3p, dme-miR-284-3p, and dme-miR-289-5p were identi ed as non-canonical miRNAs. The rest of 19 DE miRNAs were canonically conserved miRNAs.
Furthermore, functional annotation was carried out for the source genes of DE circRNAs and DE miRNAs.
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, transcription factor binding, and transcription regulator activity. In addition, different circRNAs were observed to originate from the same mRNA transcript. For example, dme_circ_0008175 and dme_circ_0008173 were originated from Nlg1 gene, and dme_circ_0009514 and d me_circ_0009500 were derived from pan gene. The Gene Ontology (GO) annotations of ageing-related DE miRNAs were investigated, and 17 DE miRNAs were clustered with 20 GO terms ( Figure 1B), including cell fate commitment, programmed cell death, developmental growth, and regulation of circadian rhythm terms.  Figure 2A). According to the trends in expression quantity changes, three DE circRNAs (dme_circ_0006913, dme_circ_0008173, and dme_circ_0009500) targeted the DE miRNAs with opposite expression trends ( Figure 2B). Based on qPCR results, the expression patterns of two circRNAs (dme_circ_0008173 and dme_circ_0009500), four miRNAs (dme-miR-289-5p, dme-miR-985-3p, dme-miR-286-3p, and dme-miR-14-5p), and four mRNAs (frizzled, CG31064, Abl, and SERCA) were consistent with RNA-seq data ( Figure 2C). Importantly, the expression trend of dme_circ_0009500/dme-miR-14-5p/SERCA, dme_circ_0009500/dme-miR-289-5p/frizzled, and dme_circ_0009500/dme-miR-289-5p/CG31064 conformed to the ceRNA mechanism.

mRNA/protein correlation in Drosophila ageing
To reveal the post-transcriptional regulatory network of gene expression in Drosophila at day 7 and day 42, the expression trends between proteins and mRNA were analysed. A relatively weak correlation (Pearson value 0.049) between mRNA and proteins was observed. A total of 2753 mRNAs and mapped proteins were found in Drosophila at day 7 and day 42 (Dataset S8). 16 pairs of interacting DE mRNAs/proteins were identi ed, including 10 pairs with the same expression trends and 6 pairs with the opposite expression trends ( Figure 5A, Dataset S8). GO annotation was performed for the 16 pairs of interacting mRNAs/proteins ( Figure 5A). FBgn0035542 (Diabetes and obesity regulated, DOR) was enriched for the RNA processing term, FBgn0012042 (Attacin-A, AttA) was enriched for the extracellular region term, and FBgn0036948 (CG7298) was enriched for the carbohydrate derivative metabolic process term; the rest of the proteins/mRNAs were enriched for the response to biotic stimulus and organonitrogen compound metabolic process terms. Evidently, the mRNA levels of FBgn0001624 (discs large 1, dlg1), FBgn0052699 (lysophosphatidylcholine acyltransferase, LPCAT), FBgn0002707 (meiotic 9, mei-9), FBgn0086907 (cytochrome c distal, Cyt-c-d), and FBgn0265101 (suppressor-of-G2-allele-of-skp1, Sgt1) were not changed between day 7 and day 42; however, the proteins encoded by these mRNAs were DE between day 7 and day 42. Furthermore, dlg1, LPCAT, mei-9, Cyt-c-d, and Sgt1 were clustered into different KEGG pathways ( Figure  5B).
10 mRNA genes were randomly selected for expression level detection via qPCR at day 7 and day 42. FBgn0034335 (Glutathione S transferase E1, GstE1), FBgn0034480 (CG16898) and FBgn0035542 (DOR) exhibited the same expression patterns as observed in the RNA-seq data ( Figure 5C). When compared to 7 day-old fruit ies, GstE1 and CG16898 were upregulated in 42 day-old fruit ies, while DOR was downregulated at day 42.

Discussion
Increasing evidence has suggested that ceRNA networks play key roles in a variety of biological processes, such as cancer [31,32], Alzheimer's disease [33], skeletal muscle myogenesis [34], and ageing [35,36]. Drosophila models are used to study multiple biological phenomena and various human diseases [7,[37][38][39][40]. Thus, identi cation and analysis of gene regulatory networks in Drosophila are very important and are necessary to improve the applicability of research on ageing and age-related diseases. In this study, DE circRNA/miRNA/mRNA networks, DE lncRNA/miRNA/mRNA networks, and the mRNA/protein correlation were analysed in 7 and 42-day-old Drosophila.
Thus, the miRNA regulatory network is far-reaching [21]. Previous studies have veri ed 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 [21,60] and a number of cancer types [61][62][63]. 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 [21,[64][65][66]. In our study, 60% (18 out of 30 total miRNAs) of DE miRNAs between day 7 and day 42 were found to affect fruit y ageing, including dme-miR-14 (shortlived, mir-14 Δ1 ) [67] and dme-miR-276a/b-5p (short-lived, miR-276a KO ) [68,69], etc (Table S1). The results indicate that these miRNAs may act as sponges of circRNAs and lncRNAs involved in regulation of Drosophila ageing.

Several of the DE miRNAs in those networks have been predicted to interact with both DE circRNAs and DE
lncRNAs, such as dme-miR-14-5p, dme-miR-276a-5p, and dme-miR-276b-5p. Based on our study ndings, dme-miR-14-5p may target XLOC_100429 and thereby regulate SERCA and Vha100-4 mRNA expression. Recent studies have demonstrated that mir-14 is an abundant miRNA in S2 cells [75], and it is one of the most engaged miRNAs in AGO1-RNP complexes [70]. Furthermore, the expression trend of dme_circ_0009500/dme-miR-14-5p/SERCA network determined by qPCR was consistent with the RNA-seq data in our study. miR-14-mutant ies present defects related to apoptosis, the stress response, survival, and metabolism [76] and are more sensitive to nutrient deprivation than wild-type ies [67]. Drosophila miR-14 affects ageing by regulating insulin production and metabolism through its target, sugarbabe [67]. In male ies, miR-14 mutation can reduce survival during nutrient deprivation [67]. Therefore, miR-14 may be a sponge of dme_circ_0009500 and XLOC_100429 to regulate Drosophila ageing. The speci c regulatory mechanism will be the subject of future research.
GO annotations of the DE circRNA/miRNA/mRNA and DE lncRNA/miRNA/mRNA networks demonstrated that those DE circRNAs and lncRNAs may play an important part in Drosophila aging via various biological processes. Apparently, multiple GO terms of DE circRNA/miRNA/mRNA networks are related to aging, such as cellular homeostasis [77] and developmental growth [78]. Speci cally, multiple GO terms are involved in regulation of sequestering of calcium ion, including inorganic ion homeostasis, cation homeostasis, and cellular calcium ion homeostasis, etc. Ca 2+ dyshomeostasis is associated with several aging-related neurodegenerative diseases, such as Alzheimer's disease (AD), Huntington's disease (HD), Parkinson's disease (PD), and amyotrophic lateral sclerosis (ALS) for altered Ca 2+ buffering capacity, altered regulation of Ca 2+ channels and pumps, and altered neuronal excitability [79]. Previous study demonstrated that aging was closely linked to dysregulation of Ca 2+ homeostasis resulting in a chronically elevated level of cytosolic Ca 2+ in experimental models of neuronal aging [80,81]. As Drosophila aging-related DE mRNA targets of lncRNA XLOC_100429, Vha100-4 and SERCA were involved in inorganic ion homeostasis, cation homeostasis, ion homeostasis, cellular cation homeostasis, and cellular ion homeostasis. The activity of Vha100-4 and SERCA protein were related to the restoration of resting cytoplasmic Ca 2+ concentrations and termination of such intracellular Ca 2+ transients [82]. The results suggest that those DE circRNA/miRNA/mRNA and DE lncRNA/miRNA/mRNA networks may play a signi cant role in in Drosophila aging.
In our present study, the mRNA/protein correlation in Drosophila ageing was analysed. The relatively weak correlation (Pearson value 0.049) may have been due to different post-translational protein regulatory mechanisms, such as protein modi cation (e.g., phosphorylation and ubiquitination) [83]. Explaining the mechanism of the same expression patterns between mRNAs and proteins is relatively simple. Ten pairs of DE mRNAs and proteins with the same expression were identi ed in this study. The mRNA and protein expression levels of AattA (short-lived with allele AttA GD17307 ) [84], FBgn0034511 (GNBP-like3, short-lived with alleles GNBP-like3 GD503 , GNBP-like3 KK106085 , and GNBP-like3 NIG.13422R ) [84], and FBgn0038074 (Gnmt, long-lived with allele Gnmt UAS.cOa ) [85] were high in Drosophila at day 7 and day 42, which affected the lifespan of the ies. In contrast, in the present study, DOR exhibited relatively low expression in ies on day 7 and day 42. DOR is antagonized by insulin signalling via the product of foxo in the fat body and forms part of a feed-forward mechanism whereby ecdysone potentiates its own signalling [86], which also regulates autophagosome formation and protein degradation via ecdysone signalling [87,88]. Therefore, the mRNA/protein correlation analysis provides an important foundation to parse the genetic process of Drosophila ageing.

Sample collection and preparation
Female ies (Dahomey WT ) were bred under a 12 hr on/off light cycle at 25°C in 50% humidity. Sample collection of adult ies at day 7 and day 42 was performed as described by Yang et al. (2016) [20]. RNA-seq (of mRNA, lncRNA, and circRNA) and proteome sequencing was conducted on two biological replicates, while miRNA sequencing was conducted on three biological replicates. All samples were stored at -80°C until use.
LncRNA, mRNA, and protein data for Drosophila at day 7 and day 42 RNA-seq data of lncRNAs and mRNAs from Drosophila at day 7 and day 42 was carried out by Yang et al. (2016) [20]. The proteins from Drosophila at day 7 and day 42 were sequenced by Zhu (2017) [89]. These original data were from our laboratory.
CircRNAs in Drosophila at day 7 and day 42 The circRNA analysis was based on the RNA-seq data for Drosophila at day 7 and day 42 from Yang et al. (2016) [20]. Overall, 95.79% clean reads were obtained from 26.7 GB of raw sequence data (SRP073695) and then aligned to the D. melanogaster genome from FlyBase (Dmel_Release_6, http://FlyBase.org/). The nd_circ [90] and CIRI2 [91] software tools were utilized to identify circRNAs. Then, the overlapping Drosophila circRNAs identi ed by both software programs were selected. The input data for 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 [92] based on a negative binomial distribution. The P value was adjusted using Hochberg and Benjamini's methods [93] to control the error discovery rate. A P value < 0.05 was considered to indicate a DE circRNA.

MiRNAs in Drosophila at day 7 and day 42
TRIzol Reagent (Invitrogen, CA) was used to extract total RNA from fruit ies 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 for accurate measurement of 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 England Biolabs, USA) following the manufacturer's recommendations. Subsequently, an Agilent 2100 Bioanalyzer and an Agilent DNA 1000 chip kit (Agilent, part #5067-1504) were utilized for accurate assessment of the quality and concentration of the sequencing libraries. The libraries were sequenced using an Illumina NextSeq 500. MiRNA Fragment Sequencing service provided by Aksomics company.
Clean reads were generated from the raw sequence data from the Illumina NextSeq instrument through realtime base calling and quality ltering. The clean reads were recorded in FASTQ format and contained the 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 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 [94]. The expression pro les of miRNAs were determined based on the counts of the reads mapped. The DE miRNAs were screened with the R package EdgeR based on the count values [95]. A fold change cutoff of 1.5 and a P value cutoff of 0.05 were applied only when replicates were used for screening of DE miRNAs. GO enrichment analysis of the DE miRNAs was implemented with the GOseq R package [96].
CeRNA analysis of lncRNA/circRNA-miRNA-mRNA pathways LncRNAs, circRNAs, and miRNAs showed signi cantly different expression levels between day 7 and day 42 and thus were analysed. The potential ceRNAs were searched based on the sequences of the lncRNAs, circRNAs, and mRNAs. O ine software of MiRanda (http://cbio.mskcc.org/miRNA2003/miranda.html) was utilized to predict miRNA binding seed sequence sites, and overlap of the same miRNA binding sites on lncRNAs-miRNAs and circRNAs-mRNAs was taken to indicate a lncRNA/circRNA-miRNA-mRNA interaction.
Network analysis of mRNA/protein correlation DE mRNAs and proteins between day 7 and day 42 were analysed. Expression regulation analysis was performed by integrating the mRNA and protein data with the Perl language. The Python program Matplotlib was used to integrate the mRNA and protein expression levels, and the Kyoto Encyclopedia of Genes and Genomes (KEGG) Orthology-Based Annotation System (KOBAS) online database [97] was used to analyse the statistical enrichment of the DE mRNAs and proteins in KEGG pathways (http://www.genome.jp/kegg/). GO enrichment analysis of DE mRNAs and proteins was implemented with the GOseq R package [96]. GO terms with corrected P values less than 0.05 were considered signi cantly enriched for the DE genes.

qPCR detection
To evaluate the quality of ceRNA networks and mRNA/protein correlation, 8 miRNAs, 19 mRNAs, 6 lncRNAs, and 3 circRNAs were selected to qualify the expression using qPCR. The samples from fruit ies at 7 days (three biological duplicates) and 42 days (three biological duplicates) were used to isolated total RNA using TRIzol™ LS Reagent (ThermoFisher) according to the manufacturer's instructions. Total RNA (1 mg) from each sample was reverse transcribed with random primers using a RevertAid First-strand cDNA Synthesis Kit (ThermoFisher) according to the manufacturer's protocol, which was utilized to detect the expression level of mRNA, lncRNA, and circRNA. In addition, total RNA (1 mg) for miRNA expression detection in each sample was reverse transcribed using a TaqMan™ MicroRNA Reverse Transcription Kit (ThermoFisher) according to the kit instructions. All gene primers (Dataset S9) were designed using 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. Gene expression levels were calculated by the cycle threshold values that were identi ed as 2-DDCT. Ribosomal protein L32 (rp 49) was used as the reference gene to calculate the relative mRNA, miRNA, lncRNA, and circRNA levels.
In summary, we screened and identi ed DE circRNAs and miRNAs in Drosophila between days 7 and 42. Then, we used the DE mRNAs, circRNAs, miRNAs, lncRNAs, and proteins between day 7 and day 42 to analyse y age-related circRNA/miRNA/mRNA and lncRNA/miRNA/mRNA networks and the mRNA/protein correlation in the context of Drosophila ageing. These regulatory networks have the potential to reveal the genetic mechanisms of y ageing and to eventually enhance our understanding of human ageing and agerelated diseases.

Declarations Supplementary Information
The online version contains supplementary material available at https://doi. org/xxxx/xxxx. Ethics approval and consent to participate Not applicable.

Consent for publication
Not applicable.

Availability of data and materials
Raw read data of miRNA sequencing is available in the Sequence Read Archive (https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=search_obj, accession number SRP311738).

Competing interests
The authors declare that they have no competing interests.

Funding
The study was supported by the National Natural Science Foundation of China (Grant No. NSFC31771338). The funding bodies had no role in the design of the study; the collection, analysis, and interpretation of the data; or the writing of the manuscript.  Expression and GO annotations of DE miRNAs in Drosophila between day 7 and day 42. A, the expression level of 30 DE miRNAs in Drosophila between day 7 and day 42; "↑" and "↓" mean upregulation and downregulation on day 42 compared to day 7 in RNA-seq data.B, GO annotations of 17 DE miRNAs in Drosophila between day 7 and day 42. Grey, no target genes of miRNA were enriched in GO terms; Color, the target genes of these miRNAs were signi cantly enriched in GO terms (P≤0.01), and the darker presents the higher of the enrichment degree. DE circRNA/miRNA/mRNA networks of Drosophila at day 7 and day 42. A, circRNA/miRNA/mRNA interaction networks; yellow, circRNAs; red, mRNAs; green, miRNAs. B, speci c circRNA/miRNA/mRNA networks; "↑" and "↓"mean upregulation and downregulation on day 42 compared to day 7 in RNA-seq data, respectively. C, data on relative expression of circRNA, miRNA, and mRNA detected by qPCR analysis. DE lncRNA/miRNA/mRNA networks in Drosophila at day 7 and day 42. A, lncRNA/miRNA/mRNA interaction networks; yellow, lncRNAs; red, mRNAs; green, miRNAs. B, speci c lncRNA/miRNA/mRNA networks; "↑" and "↓" indicate upregulation and downregulation on day 42 compared to day 7, respectively. C, data on relative expression of lncRNA, miRNA, and mRNA by qPCR analysis. 7 d, 7 days; 42 d, 42 days.  Figure 5A and 5B. A positive number indicates upregulation and a negative number means downregulation in 42 day when compared to 7 day. The "+" in Figure 5C means the same expression patterns between qPCR results and RNA-seq data.