Research Paper Volume 12, Issue 3 pp 2453—2470

The whole transcriptome regulation as a function of mitochondrial polymorphisms and aging in Caenorhabditis elegans

Yuanjian Song1, *, , Yuechen Wang1, *, , Ying Li2, *, , Liang Wang3, , WenDa Zhang4, , Jing Cheng2, , Yao Zhu1, , Haoyu Zhang5, , Qiang Zhang1, , Haichen Niu1, , Yingwei zheng6, , Mengyu Liang7, , Mengqiong Deng7, , Hao Shi7, , Hao Wang7, , Fang Zhang8, , Zuobin Zhu1, ,

  • 1 Department of Genetics, Xuzhou Medical University, Xuzhou, China
  • 2 Medical Technology School of Xuzhou Medical University, Xuzhou, China
  • 3 Department of Bioinformatics, School of Medical Informatics and Engineering, Xuzhou Medical University, Xuzhou, China
  • 4 Department of Urology, Xuzhou Central Hospital, Xuzhou, China
  • 5 School of Marine Sciences, Nanjing University of Information Science and Technology, Nanjing, China
  • 6 Department of Biochemistry, Xuzhou Medical University, Xuzhou, China
  • 7 Clinical College of Xuzhou Medical University, Xuzhou, China
  • 8 Research Facility Center for Morphology, Xuzhou Medical University, Xuzhou, China
* Equal contribution

Received: November 2, 2019       Accepted: January 7, 2020       Published: February 4, 2020      

https://doi.org/10.18632/aging.102754
How to Cite

Copyright: © 2020 Song et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

Recently, mitochondrial-nuclear interaction in aging has been widely studied. However, the nuclear genome controlled by natural mitochondrial variations that influence aging has not been comprehensively understood so far. We hypothesized that mitochondrial polymorphisms could play critical roles in the aging process, probably by regulation of the whole-transcriptome expression. Our results showed that mitochondria polymorphisms not only decreased the mitochondrial mass but also miRNA, lncRNA, mRNA, circRNA and metabolite profiles. Furthermore, most genes that are associated with mitochondria show age-related expression features (P = 3.58E-35). We also constructed a differentially expressed circRNA-lncRNA-miRNA-mRNA regulatory network and a ceRNA network affected by the mitochondrial variations. In addition, Kyoto Encyclopedia of Genes and Genomes pathway analyses showed that the genes affected by the mitochondrial variation were enriched in metabolic activity. We finally constructed a multi-level regulatory network with aging which affected by the mitochondrial variation in Caenorhabditis elegans. The interactions between these genes and metabolites have great values for further aging research. In sum, our findings provide new evidence for understanding the molecular mechanisms of how mitochondria influence aging.

Introduction

Mitochondria functions as an essential energy provider for eukaryotic cells and could occupy up to 20% cell volume, which makes it vital in cell growth and aging. Although mounting evidence suggests that the mitochondrial dysfunction is a central event of aging processes [1], some studies showed that mitochondrial DNA polymorphisms are associated with human longevity [2, 3]. However, the mechanism behind the association remains unclear. Recently, nuclear genes have been reported to be capable of altering mitochondrial functions with severe biological dysfunction [4], and the mtDNA genotype vary among different nuclear allelic backgrounds could influence healthy ageing [5]. Most studies of regulatory mechanism of the mitochondrial function have focused on nuclear genes. However, a recent study showed that the mitochondrial-encoded peptide could regulate nuclear gene expression in response to metabolic stress [6]. Currently, mtDNA variants can significant increase or decrease the expressivity of nuclear genes causing diseases which leads to the postulation that mitochondrial-to-nuclear signaling might accelerate the process of aging [7]. In contrast, comprehensive understanding of how mitochondria regulate the nuclear genome and then influence aging is still lacking.

The longevity of C. elegans is an important feature of aging and is convenient for genetic analysis, as different populations of a species can live in a standardized and uniform environment. It is now known that natural mitochondrial variant in COX1 subunit of mitochondrial complex IV (CIV) and mitochondrial physiological functions show apparent differences in between two wild isolates of the Bristol strain (N2) and the Hawaii strain (CB4856) [8]. In addition, strain CB4856 has a comparatively shorter lifespan than strain N2 [8, 9]. Dingley et al. speculated that mtDNA and mitochondrial-nuclear genome (mit-n) mismatch may lead to a reduction in lifespan [8]. Our previous studies also confirmed that the degree of mit-n mismatch was linked with the mitochondrial function, leading to the quantitative trait of C. elegans lifespan [9]. So far, it is not clear whether the shorter lifespan of CB4856 is independently caused by mitochondrial dysfunctions or the nuclear gene expressions since the nuclear genomes between CB4856 and N2 are significantly different [10]. Comprehension of how naturally occurred variations in mitochondrial DNA causing quantitative variations in phenotypic traits is a major challenge of contemporary biology. To gain further insight into mitochondrial polymorphisms is the starting point in understanding aging process. In this study, we used the wild N2 strain and the trans-nuclear cybrid C. elegans CN30 to examine how the natural variations of COX1 in mitochondria interacts with the nuclear genes and decreases the longevity of C. elegans. The nuclear transfer strain CN30 was constructed by first crossing N2 males with CB4856 hermaphrodite, then by repeatedly crossing the offspring with N2 males up to 20 generations, and finally by inbreeding up to 30 generations to generate the final CN30 strain. The strains were verified by 20 SNPs that were evenly distributed on the six pairs of chromosomes and the mitochondrial genome. Since the N2 and CN30 worms have almost the same nuclear background and they were experimented with the standardized and uniform environment, the trans-nuclear cybrid strain CN30 is a good model to study mitochondrial function and its interaction with nuclear genome.

Complex and highly ordered trait should be subject to influenced by the complex regulatory networks which could be changed by the variations [11]. Accordingly, the aging and longevity traits that normally display quantitative genetic variations should be regulated by a highly ordered biological system. Thus, we postulated that all random variations with deleterious impacts such as random noise and entropy increment, etc. should be under negative selection and lead to increased aging. Natural variations in mitochondria not only generate the noises in the ordered molecular regulatory network but also supply the necessary energy for the effective operation of the biological system. In this study, we proposed a hypothesis that longevity or aging as a complex trait is likely to be controlled by nuclear gene expression network, which could be further regulated by natural variations in mitochondrial DNA. Thereafter, construction of gene expression regulation networks regulated by mitochondrial variations can facilitate the understanding of the regulatory mechanism on aging.

mRNAs expression profiles in C. elegans have been well studied in terms of aging [12]. In addition, accumulating evidence indicates that ceRNA network may play important roles in many disease processes [13, 14]. However, the interactions between the mRNA, ncRNA and ceRNA affected by the mitochondria variation in COX1 and the potential roles of the miRNA, lncRNA, circRNA and ceRNA in the pathogenesis of aging are still unclear and have not be characterized yet. Here we used next generation sequencing (NGS) to identify mRNA and ncRNA and evaluate their associations with mitochondria and aging. Differentially expressed mRNA, lncRNA, miRNA, circRNA, and ceRNA were identified between wild strain and cybrid strain along the aging process. Regulation of the transcriptional regulatory networks in the aging process was also investigated in this study.

The molecular and cellular networks associated with aging rely on a large number of proteins that collaborate together to make adaptive and contingent metabolic responses to environmental challenges [15]. So far, several metabolic studies of aging have been reported and a series of age-related metabolites were discovered [16, 17]. Here, we aim to dissect the associations between metabolic regulation and longevity regulated by the natural mitochondrial polymorphism, which will provide translational evidence to promote human health-span at genetic level. We applied a mass spectrometry-based approach (MS) to screen global endogenous metabolites in N2 wild strain and trans-nuclear cybrid strain. Integrative analysis of transcriptome and metabolome revealed that the biological process could be influenced by natural variation in mitochondria and lead to accelerate aging.

Through high-throughput sequencing and metabolomics analysis (Supplementary Figure 1), this study provided an insight into the roles of COX1 variations in aging process through whole-scale integrative transcriptomics and global metabolomics in C. elegans, which identified important potential aging targets for future studies in aging intervention and age-related disease prevention.

Both mitochondrial and nuclear genomes correlate with the quantitative traits of development and lifespan of Caenorhabditis elegans. (A) Average lifespans in strains N2, CB4856, N2-mitotype RIAILs, CB4856-mitotype RIAILs and CN30. (B) Lifespan of 34 N2-CB4856 recombinant inbred advanced intercross lines (RIAILs). (C) Motility of 0, 6, 9 and 12-day-old worms (darker color in the heatmap indicates weaker motility). (D) Histogram of normalized body length in control (gray) and in paraquat (black) conditions for 210 RIAILs. Statistical analysis was performed by using two-tailed unequal variant Student’s t-test (* P-value

Figure 1. Both mitochondrial and nuclear genomes correlate with the quantitative traits of development and lifespan of Caenorhabditis elegans. (A) Average lifespans in strains N2, CB4856, N2-mitotype RIAILs, CB4856-mitotype RIAILs and CN30. (B) Lifespan of 34 N2-CB4856 recombinant inbred advanced intercross lines (RIAILs). (C) Motility of 0, 6, 9 and 12-day-old worms (darker color in the heatmap indicates weaker motility). (D) Histogram of normalized body length in control (gray) and in paraquat (black) conditions for 210 RIAILs. Statistical analysis was performed by using two-tailed unequal variant Student’s t-test (* P-value < 0.05, ** P-value < 0.01).

Results

Both mitochondrial and nuclear genomes correlate with the quantitative traits of development and lifespan of C. elegans

Comparing to wild-type strain N2, a single nonsynonymous SNV occurs in COX1 of the wild-type strain CB4856 via the replacement of alanine with serine (p.A12S) [8]. To determine whether COX1 variation in CB4856 could independently reduce lifespan than N2, the mitochondrial origin and lifespan of the 34 recombinant inbred advanced intercross lines (RIAILs) was quantified. There was no significant difference in the averaged lifespans between CB4856 and N2 mitotypes RIAILs strains. The trans-nuclear cybrid strain CN30 which carried CB4856 mitochondria with N2 nuclear background was then generated, which was verified by 18 SNPs evenly distributed on the six pairs of chromosomes. The lifespan of CN30 was significantly shorter than those strains carrying CB4856 mitochondria but not N2 nuclear background (Figure 1A). In addition, among the 34 RIAILs, 13 strains had a shorter lifespan than CB4856, 5 strains had a longer lifespan than N2, and 16 strains showed intermediate lifespans (Figure 1B). Excessive shorter lifespan of CN30 compared to wild-type strains suggests that there may be some genetic incompatibility between these RIAILs.

Motility is an important physiological and functional indicator of health-span [18, 19]. In this study, motility of 0, 6, 9 and 12-day-old C. elegans was quantified according to liquid media spontaneity in wave-like channels. Decline in motility from day 6 was significantly accelerated for CN30 strain when compared with the other two wild-type isolates. The RIAILs that have partial mit-n mismatches also showed faster rates of motility reduction than wild-type isolates (Figure 1C). Previous studies have shown that growth rate (as measured by the length of offspring or optical density) varies substantially in 210 RIAILs with defined nuclear and mitochondrial genotypes [20]. Using these data, we identified 160 RIAILs with relatively lower growth rate compared to that of CB4856, 35 RIAILs with relatively higher growth rate compared to N2, and 12 RIAILs with intermediate growth rate (Figure 1D). Moreover, RIAILs growing in paraquat conditions showed the same trend as described above. These results confirmed that the mit-n mismatches could interfere the development and aging of C. elegans.

Variation in COX1 decrease mitochondrial mass with age

In this study, we examined how the variation in COX1 decrease mitochondrial mass with age. Oxygen consumption rate (OCR) of mitochondria is dependent on many biological reactions and has emerged as a molecular markers for mitochondrial health [21, 22]. Our results showed that OCR decreases with age in N2, CB4856 and CN30 worms, and OCR of CN30 have a greater reduction than N2 (Figure 2A). The relative mtDNA genome content of the two wild-type isolates and the cybrid CN30 strain were assessed by qPCR in 0, 6, 9 and 12 day old worms [9]. There was no difference of mtDNA content among the three strains in 0, 6, and 9 day old worms. However, in the 12-day old worms, N2 strain had a significant reduction in mtDNA content relative to CN30 (Figure 2B). These results suggested that the variation in COX1 led to a faster decline in OCR, which was not related to mitochondria copy number. Although the function of free radical compounds in aging is still controversial, aging is generally accompanied with an increase in reactive oxygen species (ROS) production [23]. Thus, we measured ROS content per worm and found that ROS content was significantly increased with age, especially in CN30 (Figure 2C). In addition, old CN30 worms have more mitochondrial copy number, lower OCR value, and higher ROS content than N2 worms, which suggested that variation in COX1 could reduce mitochondrial functions and lead to aging.

Variation in COX1 decreases mitochondrial mass with age. All measured mitochondrial indicators were normalized to the value of day 0 N2 worms. (A) Level variation of OCR with age. (B) Level variation of mtDNA with age. (C) Level variation of ROS with age. (D) ATP level variation with age. Data shown as the means ± SEM. Statistical analysis was performed by using two-tailed unequal variant Student’s t-test (* P-value

Figure 2. Variation in COX1 decreases mitochondrial mass with age. All measured mitochondrial indicators were normalized to the value of day 0 N2 worms. (A) Level variation of OCR with age. (B) Level variation of mtDNA with age. (C) Level variation of ROS with age. (D) ATP level variation with age. Data shown as the means ± SEM. Statistical analysis was performed by using two-tailed unequal variant Student’s t-test (* P-value < 0.05, ** P-value < 0.01).

ATP production levels were also examined to determine the quality of mitochondrial functions. Although ATP level decreases with age, the results did not show significant difference among N2, CB4856 and CN30 worms (Figure 2D). In specificity, mtDNA content of CN30 showed smaller downward trend with age-related to N2 strain. ATP levels of CN30 downward trend is more obvious (60% VS 49%). When combined these results together, we can speculate that an organism needs a certain amount of energy to accomplish its basic physiological activities, as COX1 can reduce mitochondrial function, more mitochondria are needed for energy, then more ROS produced by more dysfunctional mitochondria.

Variation in COX1 linking with nuclear genome SNPs and up-regulate or down-regulate mRNA, lncRNA, miRNA and circRNA profiles

The CN30 and CB4856 worms have the same mitochondria genotype with similar molecular functions but different lifespans, indicating that impaired mitochondrial function caused by COX1 variation alone cannot accelerate the aging process. Thus, nuclear genes associated with aging might be a function of mitochondrial polymorphisms in C. elegans. Our previous studies indicated that abundant loci and genes in the nuclear genome might be associated with mitochondria variations and the quantitative traits of lifespan in C. elegans [9]. By further comparing the transcriptomes of CN30 and N2 worms with the same nuclear background but different mitochondrial genotypes, we can study how mitochondria up-regulate or down-regulate nuclear genes and then influence the aging process.

Here, we detected which loci and genes were linked with mitochondrial variations by using the high-throughput sequencing data and the published SNPs chip data [10]. Single nucleotide polymorphisms in nuclear genome that might be linked with mitochondrial polymorphisms were searched for by using case/control association option of the PLINK software [30], which could estimate whether nuclear loci have epistasis effects with mitochondrial polymorphisms. In the 93 strains with npr-1 CB4856 genotype, 92 out of 1454 SNPs were found to have linkages with mitochondrial polymorphisms (P < 0.01, q-value <= 10%, Supplementary Table 1). The genotype of most SNPs associated with mitochondria (CB4856 or N2 genotype) matched with the mitochondria genotype (CB4856 or N2) (matched vs unmatched: 62 vs 30, P < 0.05, Supplementary Table 1). The results suggested that there exists a correlation between mitochondrial genome and nuclear genome.

The expression profiles of mRNAs and lncRNAs were generated by Illumina X ten in wild strain N2 and the trans-nuclear cybrid strain CN30. The expression profiles of miRNAs were measured by BGISEQ-500. RNAs with a threshold value of fold-change >2 and q-value ≤ 0.001 were identified as significantly differential expression genes. The expression levels of 801 mRNAs, 245 lncRNAs, and 26 miRNAs were significantly altered in CN30 compared to N2 controls (Figure 3A3C), among which 584 mRNAs, 53 lncRNAs, and 16 miRNAs were up-regulated, while 217 mRNAs, 192 lncRNAs, and 10 miRNAs were down-regulated. CircRNA which is a large class of noncoding RNA and plays important roles in gene expression regulation were also predicted, 76 circRNAs were significantly differentially expressed in between CN30 and N2 worms (Figure 3D). The differentially expressed mRNAs, lncRNAs, circRNA and miRNAs in the CN30 worms are all listed in Supplementary Table 2.

Influences of variation in COX1 on mRNA, lncRNA, circRNA and miRNA expression profiles. Heatmaps of the differentially expressed (A) mRNAs, (B) lncRNAs, (C) circRNAs, and (D) miRNAs in between CN30 and N2 worms. The data are depicted as matrices in which each row represents one mRNA, lncRNA, miRNA, or circRNA and each column represents one of the worm samples. Relative mRNA, circRNA, miRNA, or lncRNA expression is depicted according to the color scale shown at the top. Red and green colors represent relatively high and low expression values. The magnitude of deviation from the median is represented by color saturation.

Figure 3. Influences of variation in COX1 on mRNA, lncRNA, circRNA and miRNA expression profiles. Heatmaps of the differentially expressed (A) mRNAs, (B) lncRNAs, (C) circRNAs, and (D) miRNAs in between CN30 and N2 worms. The data are depicted as matrices in which each row represents one mRNA, lncRNA, miRNA, or circRNA and each column represents one of the worm samples. Relative mRNA, circRNA, miRNA, or lncRNA expression is depicted according to the color scale shown at the top. Red and green colors represent relatively high and low expression values. The magnitude of deviation from the median is represented by color saturation.

Age-dependent changes in gene expression

Among the 21473 genes sequenced by Illumina X ten platform, we identified 9762 genes with significantly differential expressions between young and elder worms. We then divided these genes into 4 clusters according similar time-dependent changing trend (Figure 4A4D). All the genes that regulate aging must change with aging in terms of expression profiles. We then examined whether genes affected by the variation of mitochondria could have a similar correlation with aging. Most of the genes that were associated with mitochondria (Figure 3A) showed age related expression feature (542/801 vs. 9762/21473, P-value = 3.58E-35, two-tailed Fisher exact test) (Figure 4E). These results suggest a role of the natural variation of mitochondria in aging via regulation of multiple nuclear genes that might link with aging. We also assessed the influence of age at whole genome expression level between N2 and CN30. In developing worms, by analyzing the gene expression level of 21473 genes, more genes showed a high expression level than adult and elder worms (Figure 4F). For 801 genes that were significant up-regulated or down-regulated by mitochondria, the average expression level was significantly lower than the whole genome gene expression level. For CN30 worms, by analyzing the gene expression level of 801 genes, the number of higher expression genes were more than N2 worms (646 vs 274) (Figure 4G). These results suggested that mitochondrial polymorphism lowered basal expression levels of relevant genes and participated in the regulation of aging.

Age-dependent changes in gene expression levels. (A–D) Clusters of genes associated with age. Genes with the same expression pattern were grouped into a cluster. Yellow and green colors correspond to low membership values, while red and purple colors correspond to high membership values. (E) Venn diagram of the differentially expressed genes. Each circle represents a set of genes. The left circles represent the genes associated with mitochondria variation. The right circles represent the genes associated with aging. The region superimposed by circles represent the genes both associated with mitochondria variation and aging. (F) Gene expression features at the whole genome level. High expression, intermediate expression, and low expression are represented by FPKM>=10, FPKM=1~10 and FPKMG) Gene expression features of 920 genes affected by mitochondria polymorphism.

Figure 4. Age-dependent changes in gene expression levels. (AD) Clusters of genes associated with age. Genes with the same expression pattern were grouped into a cluster. Yellow and green colors correspond to low membership values, while red and purple colors correspond to high membership values. (E) Venn diagram of the differentially expressed genes. Each circle represents a set of genes. The left circles represent the genes associated with mitochondria variation. The right circles represent the genes associated with aging. The region superimposed by circles represent the genes both associated with mitochondria variation and aging. (F) Gene expression features at the whole genome level. High expression, intermediate expression, and low expression are represented by FPKM>=10, FPKM=1~10 and FPKM<=1, respectively. (G) Gene expression features of 920 genes affected by mitochondria polymorphism.

Interactive network of mRNA, lncRNA, miRNA and circRNA, and ceRNAs

To better understand the functions of non-coding RNAs and the corresponding genetic regulatory networks affected by COX1, the mRNA-lncRNA-miRNA-circRNA network was visualized by software CytoScape [24]. The possible miRNA targets (Supplementary Table 3) and the lncRNA targets (Supplementary Table 4) were predict. Relationships between circRNA and lncRNA were also determined by the position of circRNA and lncRNA on chromosomes. Finally, according to the interactions among differently expressed mRNAs, lncRNAs, miRNAs and circRNAs, complex regulatory networks were generated (Figure 5A).

Interactive network of mRNA, lncRNA, miRNA and circRNA, and ceRNAs. (A) The regulatory network was constructed through visualizing the relationships between differentially expressed circRNA, lncRNA, miRNA and mRNA and their target genes. (B) The interaction network among the top 10 ceRNAs network.

Figure 5. Interactive network of mRNA, lncRNA, miRNA and circRNA, and ceRNAs. (A) The regulatory network was constructed through visualizing the relationships between differentially expressed circRNA, lncRNA, miRNA and mRNA and their target genes. (B) The interaction network among the top 10 ceRNAs network.

Competing endogenous RNA (ceRNA) can compete with other RNA transcripts for the same microRNA. Hence achievement of mutual communication and regulation. Thus, we predicted the candidate ceRNAs and the top 10 candidate ceRNAs which have the most correlation with other RNAs were selected (Supplementary Table 5). Then, the interaction network between the top 10 ceRNAs, shared miRNAs, and target RNAs was constructed (Figure 5B). Since most of the genes associated with mitochondria showed age related expression feature, studying on the composition of the gene regulatory network is also of great value for aging research.

Variation in COX1 affects metabolic activity

KEGG pathway analysis showed that 14 pathways were significantly altered by genes linked with mitochondria variation (q-value < 0.05), among which 9 pathways were related with metabolic activity (Figure 6A). For age-related genes, 18 pathways were significantly altered by genes linked with age (q-value < 0.05), 5 pathways of which were related to metabolic activity (Figure 6B). These results suggested that the variation in mitochondria regulates aging process mainly by influencing metabolic activity.

Variation in COX1 affects metabolic activity. (A) The significantly influenced KEGG pathways by mitochondria variation. X-axis: enrichment radio. Larger value represents greater degree of enrichment. Rich ratio indicates the number of differentially expressed genes located under the path term divided by the number of all annotated genes under the path term. Y-axis: pathway. Color: P-value. Circle size: number of genes. (B) Significantly influenced KEGG pathway by age. (C) Heat map of differentially expressed metabolites. The data were depicted as matrices in which each row represents one metabolite and each column represents one of the worm samples. Relative metabolite expression is depicted according to the color scale shown at the right side. Red and green colors represent high and low relative expression, respectively. The magnitude of deviation from the median is represented by color saturation. (D) Interactions between genes and metabolites both affected by the mitochondria variation. Pearson correlation coefficients that represents the association between genes and metabolites were depicted according to the color scale shown at the right side. Red and blue colors represent positive and negative correlations. The magnitude of deviation from the median is represented by color saturation.

Figure 6. Variation in COX1 affects metabolic activity. (A) The significantly influenced KEGG pathways by mitochondria variation. X-axis: enrichment radio. Larger value represents greater degree of enrichment. Rich ratio indicates the number of differentially expressed genes located under the path term divided by the number of all annotated genes under the path term. Y-axis: pathway. Color: P-value. Circle size: number of genes. (B) Significantly influenced KEGG pathway by age. (C) Heat map of differentially expressed metabolites. The data were depicted as matrices in which each row represents one metabolite and each column represents one of the worm samples. Relative metabolite expression is depicted according to the color scale shown at the right side. Red and green colors represent high and low relative expression, respectively. The magnitude of deviation from the median is represented by color saturation. (D) Interactions between genes and metabolites both affected by the mitochondria variation. Pearson correlation coefficients that represents the association between genes and metabolites were depicted according to the color scale shown at the right side. Red and blue colors represent positive and negative correlations. The magnitude of deviation from the median is represented by color saturation.

We further examined the metabolic changes in the elder worms by using mass spectrometer Xevo, G2-XS QTOF (Waters, UK) for mass spectrometry data collection. Identification of metabolites was based on database KEGG. The variable important for the projection (VIP) values (VIP≥1) of the first two principal components of the multivariable PLS-DA model was used to screen the differential metabolites by combining fold-change (≥1.2 or ≤0.8333) and q-value (<0.05). Analysis of 12 samples (9-day old) yielded 7133 total ions by positive ion mode and relative standard deviation (RSD) <=30%, among which 474 metabolites were differentially regulated between N2 and CN30 (Figure 6C). The metabolic pathways of differential metabolites were further annotated based on the KEGG database (Supplementary Table 6).

In order to integrate the data of differentially expressed genes and differentially present metabolites, we thoroughly analyzed these genes and metabolites involved in the regulation process, which formed a high-confident gene expression regulation network. Cor function in R language was used to calculate the Pearson correlation coefficient between 474 metabolites and 135 genes that were linked with age and were differentially expressed in between CN30 and N2 at both young and elder stages. According to the correlation coefficient, Pheatmap and RColorBrewer packages were used in order to generate the correlation clustering heatmap (Figure 6D) (Supplementary Table 7). These correlation analyses can provide important guidelines for future studies on aging regulation.

Discussion

Although mitochondrion only harbors a small set of genes compared with the organism's genome, it provides the necessary energy for cells and plays important roles in signal transduction, cell growth, and apoptosis, etc. Harman first proposed the "mitochondrial aging theory" in 1956, arguing that mitochondria are the core component of aging that is accompanied by an increase in mitochondrial DNA mutations [23]. Our previous studies have shown that the degree of mismatch between mitochondrial and nuclear genomes may be one of the important factors influencing the differences in lifespans among worm populations [9]. In this study, we showed that most individuals in the RIAILs hybrid population have lower lifespan and development level than their parents. We also confirmed the hypothesis proposed by Dingley et al. that mtDNA-nuclear genome mismatch may lead to reduced lifespan of C. elegans by using the trans-nuclear cybrid strain which have almost the same nuclear background with the N2 strain [8]. We further studied the various mitochondrial mass in different life stages of worms and explored how mitochondria variation regulated the nuclear genome and influenced the aging process.

The nuclear transfer strain CN30 was constructed by first crossing N2 males with CB4856 hermaphrodite, then by repeatedly crossing the offspring with N2 males up to 20 generations, and finally by inbreeding up to 30 generations to generate the final CN30 strain. The strains were verified by 20 representative SNPs that were evenly distributed on the six pairs of chromosomes and the mitochondrial genome. Since the N2 and CN30 worms have almost the same nuclear background and they were experimented with the standardized and uniform environment, the trans-nuclear cybrid strain CN30 is a good model to study mitochondrial function and its interaction with nuclear genome. In this study, we confirmed the hypothesis proposed by Dingley et al. that mtDNA-nuclear genome mismatch may lead to reduced lifespan of C. elegans [8]. We further studied the various mitochondrial mass in different life stages of worms and explored how mitochondria variation regulated the nuclear genome and influenced the aging process.

Strain CN30 has a reduced OCR compared to strain N2 in both young and elder worms, which suggests that variation in mitochondria can alter energy metabolism. This is consistent with previous studies [8]. However, there was no difference in ATP content between CN30 and N2 accompany age. Since the process of ATP production is aerobic, there must be more mitochondria with lower OCR in order to produce a certain amount of ATP. Based on this corollary, we further examined the copy number of mitochondrial DNA (mtDNA) in different life stages of worms. We found that mtDNA content has no significant decrease in developing worms. Meanwhile, N2 had a significant reduction in mtDNA content relative to CN30 for the elder worms. The cybrid strain with mit-n mismatch also posted significantly increased CIV activity [8]. By combining these results, we proposed that impaired mitochondrial function need mobilizes the maximum potential of producing ATP and increases quantity to supply the necessary energy for basic life activities. Lifespan is a complex biological trait that is bound to be regulated by a complex network involving many genes. We postulated that mitochondria should be able to not only provide energy for life activities but also supply energy for the gene regulatory networks. Thus, it is valuable to study the mechanism of mitochondrial regulation of aging by exploring the age-related gene regulatory networks affected by mitochondria polymorphisms.

Previous studies suggested that six quantitative trait locus (QTLs) in the worms were associated with lifespan [25]. In all the six QTLs, N2 genotypes have increased lifespan while CB4856 genotypes have decreased lifespan. Strain CNC30 with N2 nuclear background and CB4856 mitochondria have lower lifespan than strain CB4856, which suggests that mitochondrial genotypes also up-regulate or down-regulate nuclear gene expression and in turn decreases worm lifespan. Here we assumed that aging-related genes should be changed along with aging in terms of expression levels. Therefore, we studied the whole-transcriptome changes in gene expression affected by the natural variation of mitochondria and age of C. elegans. We identified a group of genes affected by mitochondria polymorphism, most of which showed changed expression levels with age. Since mitochondria consisting of about 1000 proteins and most of them are encoded by nuclear genes, previous studies mainly focused on how nuclear genes affect the aging process [4, 9, 26]. Our results suggested that many genes can be actually affected by mitochondria directly and were involved in aging regulation.

Cumulative evidence indicates that ceRNA and circRNAs can function as endogenous sponges to influence mRNA, lncRNA or miRNA activity [27, 28]. We found many ncRNAs were significant up-regulated or down-regulated by the mitochondria polymorphism in C. elegans, which include novel miRNAs and lncRNAs and predicted circRNAs and ceRNAs. The miRNAs-lncRNAs-mRNA-circRNAs gene regulatory network and the ceRNA network affected by mitochondria was constructed. In addition, due to vastly changed environmental stresses, organisms might make corresponding adjustments at molecular level [15] (especially mitochondria variation). Thus, facing the genetic pressure (especially mitochondria variation), gene regulation network that regulates complex traits would also be influenced. The application of predicted circRNAs and ceRNAs and the gene regulatory network analysis may provide a novel view of the interactions between mRNAs and ncRNAs. Our results provided the first assessment of gene regulation networks influenced by mit-n and identified certain known and novel mRNA, lncRNA, miRNA, circRNAs and ceRNAs that were significant up-regulated or down-regulated by mitochondria and involved in aging.

To better understand the biological functions and potential mechanisms of the genes affected by mitochondria and age, we also performed KEGG pathway analysis. The significantly relevant KEGG pathways were metabolic pathways. Combining with the fact that metabolic level could change significantly with age [16], these results indicated that mitochondria variation could change not only the metabolism of the mitochondria but also the metabolism of the whole organism, generating impacts on the aging process. We also examined elder C. elegans metabolism and found 474 metabolites changing between N2 and CN30. Most of these metabolites belong to the metabolic pathway. Multi-omics integration analysis between metabolomics and transcriptomics was then performed. Our study established the relationship between the natural mitochondria variation and aging from the perspectives of nuclear genome regulatory network and metabolomics.

Conclusions

This study provides important genetic evidence for the hypothesis that mitochondria regulate aging through its functions and effects on nuclear gene regulatory networks. To further verify this hypothesis, the roles of genes and metabolites in the gene regulation network should be investigated and clarified in the aging process. Since promotion of ATP synthesis can help maintain healthy life activities [29], analysis of transcriptome expression after the improved mitochondrial function is also needed for further study.

Materials and Methods

Strains and media

C. elegans RIAILs were gifts from Professor Leonid Kruglyak (Princeton University, Princeton). C. elegans were fed with E. coli OP50 at 20 °C on normal growth medium (NGM).

Association mapping

SNPs associated with mitochondrial genotype were analyzed by using software package PLINK (v1.07) with the quantitative trait association option [30]. SNPs in stronger linkage disequilibrium with other SNPs were not removed as all these SNPs can be used to discern the genomic extent of intervals associated with traits [10].

Lifespan assays

Worm lifespan assays were performed by following the manufacturing instructions and are briefly described as follows. First, 100 larva stage 4 worms were transferred to NGM containing E. coli OP50, and three culture dishes were set for each line. Then, worm survival condition was scored every day by observing the pharyngeal movement and the touch-provoked movement using platinum wire. Worms were considered to be dead if there was no pharyngeal movement and no touch-provoked movement. Each experiment was duplicated three times.

Motility assays

Synchronized worms were separately transferred to food-free 3 ml NGM medium at 20 °C for three minutes. Then, CCD camera was used to record the motility of worms in M9 buffer continuously and the number of head sway from left to right in 30 seconds was also recorded. All the strains of worms used in the study were examined at the same stages and each experiment was repeated three times.

Assay for the origin of nuclear and mitochondria

Mitochondrial genomes from the two wild worms N2 and CB5856 were compared, and two SNPs with different loci were randomly selected [9]. Then, the mitochondria DNA was extracted and amplified with two pairs of primers (SNP1: AGAATGATTTACGTTACCA/TTATT TTTTTGATTTT, A = N2, T = HW; SNP2; AGAATGA TTTACGTTACCA/TTATTTTTTTGATTTT, A = N2, T = HW). DNA sequences of the nuclear genomes were compared, and 18 SNPs that evenly distributed on the six pairs of chromosomes were selected. Nuclear DNA was extracted and amplified with 18 pairs of primers.

Measurements of relative mtDNA content per cell.

Real-time PCR was used to determine the relative copy number of mtDNA content per cell [8]. Since nuclear genomic DNA is relatively stable at different stages of life, it was used as a housekeeping gene [31]. The method is briefly described as follows. A single cleaned worm was picked up into a PCR tube containing 10 μ L 100 μ g/ml proteinase K. Then, the PCR tube was tightly sealed, and it was heated at 56 °C for 15 min and 95 °C for 10 min. Then, 10 μ L SYBR Green Supermix (BIO-RAD iQ TM SYBR Green Supermix, cat. No. 170-8882AP) were added to the PCR tube and amplified with primers specific to either nuclear DNA (F, TGGAACTCTGGAGTCACACC, R, C ATCCTCC TTCATTGAACG G) or mtDNA (F, GTTTATGCTGC TGTAGCGTG, R, CTGTTAA AGCAAGTGGACG AG) [31]. The 10 synchronized 0-day, 6-day, 9-day and 12-day old adult worms were examined for each strain and the mean values were used.

Assays for ATP levels ATP

Levels assays were performed by the manufacturing instructions of ATP bioluminescent assay kit (Sigma, product number FL-AA) [9]. ATP concentration in the synchronized 0-day, 6-day, 9-day and 12-day old adult worms was examined. The ATP content were normalized by the number of worms and each experiment was repeated 3 times.

Assays for ROS levels

The worms stained with DCFDA solution at a final concentration of 50 μM DCFDA for 30 min at 20°C [32]. Twenty synchronized 0 days, 6 days, 9 days and 12 days old worms were picked up onto microscope slides which coated with 3% agarose and anaesthetized with 2% sodium azide. Worms were photographed by fluorescence microscope (100× magnifications and 200 s exposure time). The fluorescence intensity of all images was calculated by software Image-Pro Plus.

Assays for the rate of oxygen consumption (OCR) levels

OCR levels for whole C. elegans in the synchronized 0-day, 6-day, 9-day and 12-day old adult worm were examined using an analyzer Seahorse XF24 [33]. All the Worms were washed 4 times in M9 buffer to segregate clean bacteria-free adult worms from bacterial debris. Worms were loaded about 30 worms to each well in M9 buffer. Three OCR readings were taken after 20 min equilibration. Data was normalized by the number of worms per well and each experiment was repeated three times.

Total RNA extraction and non-coding RNA library construction

Total RNA was extracted from the thousands of synchronized worms using Trizol (Invitrogen, Carlsbad, USA) according to the manufacturing instructions. RNA was qualified and quantified using a Agilent 2100 and NanoDrop bioanalyzer (Thermo Fisher Scientific, USA). Qualified lncRNA libraries were sequenced via pair-end sequencing on the Hiseq xten platform [34]. Qualified small RNA libraries via pair-end sequencing on the BGISEQ-500 platform (BGI-Shenzhen, China) [35].

lncRNA-seq data analysis

Raw reads with rRNA, low quality, joint contamination, and high content of unknown base Ns were filtered out (Supplementary Table 8). Then, clean reads were matched to the reference genome using HISAT [36] and assembled using StringTie [37]. After obtaining the new transcript, the mRNA and lncRNA were distinguished by using three software CPC [38], txCdsPredict, CNCI [39] and the database Pfam [40]. Scores should exceed the threshold (CPC_threshold = 0, CNCI_threshold = 0, txCdsPredict_threshold = 500). If transcripts can be compared to Pfam database, they were conserved to be mRNA, otherwise they were conserved to be lncRNA. At least three of the four criteria were consistent, and then it was confirmed that the transcript was an lncRNA or mRNA. NT, NR, KOG, KEGG, and SwissProt were used for annotation by software Diamond [41]. GO annotation was performed by using NR annotation results and the Blast2GO software [42]. Clean reads were compared to the reference sequence by using the Bowtie 2 software [43]. Then software of RSEM that a model generated by reads as used to calculate the genes expressions [44]. Transcripts and clean reads were allocated to different transcription by using maximum likelihood (ML) method. The chain-specific model used to distinguish the reads from the positive and negative strands could facilitate a more accurate quantitative analysis of gene expression. The expression levels of genes were standardized by FPKM. Differential genes between groups were analyzed via software DEGseq [45]. Genes with change fold >= 2 and q-value <= 0.001 were identified as with significantly different expression levels between groups in this study.

Small RNA-seq data analysis

The impurities of raw data are 5' primer contaminants, no insert tags, oversized insertion, low quality tags, poly A tags and small tags, among which no insert tags and 5' primer contaminants are defined as adaptor contaminants, oversized insertion manifests as missing 3' primer. Generally, adaptor contaminant is related to the sample itself and the density proportion between the adaptor and the sample. The higher the proportion is, the worse the contaminant is. Low quality tags which have more than four bases whose quality is less than ten and which have more than six bases whose quality is less than thirteen will be filtered to ensure that the tags to analyze are reliable. We will get rid of the above contaminant tags and low-quality tags from the FASTQ file and get the final clean tags. After filtering, clean tags were mapped to sRNA database. The characteristic hairpin structure of miRNA precursor can be used to predict novel miRNA. We use miRDeep 2 [46] to predict novel miRNA by exploring the secondary structure, the Dicer cleavage site and the minimum free energy of the unannotated sRNA tags which could be mapped to the genome. The small RNA expression level is calculated by using TPM [47].

The time series analysis of gene expression

The time series analysis of gene expression was carried by the software Mfuzz [48]. At the same time, genes with silimar expression may take part in the same biological function and help to discover regulatory genes. The membership value of Mfuzz which is used to evaluate whether the gene conforms to the trend of the gene cluster is the value between 0 and1. The membership value of the gene closer to 1, the more the trend of the gene is consistent with the trend of the gene cluster.

KEGG pathway analysis

KEGG pathway were used to divide the functional modules of the genes. If the candidate gene is significantly enriched in a functional module, the proportion of the candidate genes in the functional module is significantly higher than the overall genes. Each functional module has a p-value, and the smaller the p-value, the richer the candidate genes are in this functional module. Then q-value correction was performed on p-value, and the function of q-value <= 0.05 was set as significant enrichment.

lncRNA target prediction

The function of lncRNA is mainly achieved by acting on the target gene by cis or trans. lncRNA was determined as cis within 10k upstream of mRNA or 20k downstream of mRNA. Beyond this range, RNAplex was used to analyze the binding energy of lncRNA and mRNA [49]. If the binding energy was less than -30, it was considered to be trans. The two correlation coefficients of lncRNA and mRNA, spearman and Pearson, were calculated. Spearman_cor >= 0.6 and Pearson_cor >= 0.6 mean that there is a correlation between lncRNA and mRNA.

miRNA target prediction

In order to find possible miRNA targets, multiple software was used. Taking the intersection targets with appropriate filter conditions such as MFE, score for further analysis. Generally, we use RNAhybrid [50], miRanda [51] and TargetScan to predict miRNA target.

ceRNAs prediction

The significant miRNAs shared between candidate ceRNA and mRNA were calculated, and the significant ceRNAs were predicted by combining the ratio of the shared miRNA by ceRNA and miRNA, the ratio of the shared microRNA response element (MRE) by miRNA and ceRNA, the density of shared MRE by ceRNA and mRNA, and the evenness of the shared MRE by ceRNA and mRNA.

CircRNA prediction and annotation

This study used CIRI [52] and find_circ [53] to predict circRNA, and integrated the results based on the start and end positions of the circRNA (both the start and end positions are within 10 bases before and after the circRNAs are merged into one class). CIRI is a software developed based on the BWA-MEM genome alignment algorithm, which detects circRNA by scanning the generated SAM (Sequence Alignment/Map) file. The two parts before and after the same read need to be cross-aligned to the genome in order to prove that they are junction reads at both ends of the circRNA. Therefore, during the first pass, CIRI will filter the junction reads with paired chiastic clipping (PCC) signals and then limit the range of paired-end mapping (PEM). The GT-AG shear signal is further filtered for the junction reads. In the second scan of the SAM file, CIRI will detect new junction reads and eliminate false positive circRNAs caused by mismatches to homologous genes or repeats, thereby improving the sensitivity and accuracy of the software [52]. Please refer to Figure 3 for the specific process.

Gene regulatory network construction

The circRNA-lncRNA-miRNA-mRNA network was constructed and visually displayed by using CytoScape software [24]. The top 10 candidate ceRNAs which have the most correlation with other RNAs were selected for visualization in the interaction network between ceRNAs, shared miRNAs, and target RNAs. The edges of the network reflect the interactions among different expressions of circRNA, lncRNA, miRNA, mRNA, ceRNAand their target genes.

Metabolomics analysis

Whole-scale metabolites were detected via mass spectrometer Xevo, G2-XS QTOF (Waters, UK) and identified based on the database KEGG. Software Progenesis QI (version 2.2) (Waters, UK) and R package metaX were used for statistical analysis of mass spectrometry data [54]. Variable Importance in the Projection (VIP) of the first two principal components in the multivariable PLS-DA model was used to screen differential metabolites by combining fold-change and q-value values. The VIP-score is a quantitative measure that indicates the strength and explanatory ability of each metabolite on the classification discrimination of samples in each group. Screening conditions are VIP≥1, fold-change ≥1.2 or ≤0.8333 and q-value<0.05. Metabolites satisfy the three conditions were identified as differential metabolites.

Association analysis between metabolites and genes

The cor function in R language was used to calculate the Pearson correlation coefficient between the metabolites and the genes that are differentially expressed between CN30 and N2. According to the above correlation coefficient, Pheatmap and RColorBrewer packets are used for cluster analysis, the result of which was used for the correlation matrix heatmap visualization.

The parameter description of the bioinformatics software used in this study

The parameters of the software were listed in Supplementary Table 9.

Ethical approval

This article does not contain either human nor animal experiments.

Author Contributions

Z.B.Z and Y.L designed and wrote the manuscript. Y.J.S and Y.C.W conceived the experiment, analyzed and interpreted the data. J.C, L.W, Q.Z, H.C.N, Y.Z and Y.W.Z conceived the experiment. Y.H.Z, M.Y.L, M.Q.D, H.S, F.Z and H.W analyzed and interpreted the data. All authors discussed and commented on the manuscript.

Acknowledgments

The experiments in this article were completed in Public Experimental Research Center of Xuzhou Medical University, and medical genetics research center of central south university, we thanks the teachers for their support and help during the experiments. We thank E. Andersen and L. Kruglyak for the research materials. We thank BGI Shenzhen for metabolomics experiment, and BGI-Wuhan for lncRNA and sRNA sequencing.

Conflicts of Interest

All other authors declare that they have no competing interests.

Funding

This work was supported by the National Natural Science Foundation of China (81701390, 81902040, 81601198), the Key Research and Development Plan of Jiangsu Province (BE2016641), Natural Science Foundation of Jiangsu Province (BK20170250), Xuzhou science and technology innovation project (KC19057, KC16SH024), Research Foundation for Advanced Talents of Xuzhou Medical University (D2017003) and Postdoctoral Foundation of Jiangsu Province (2018K245C).

References

  • 1. Kauppila TE, Kauppila JH, Larsson NG. Mammalian Mitochondria and Aging: an Update. Cell Metab. 2017; 25:57–71. https://doi.org/10.1016/j.cmet.2016.09.017 [PubMed]
  • 2. Niemi AK, Hervonen A, Hurme M, Karhunen PJ, Jylhä M, Majamaa K. Mitochondrial DNA polymorphisms associated with longevity in a Finnish population. Hum Genet. 2003; 112:29–33. https://doi.org/10.1007/s00439-002-0843-y [PubMed]
  • 3. Ross OA, McCormack R, Curran MD, Duguid RA, Barnett YA, Rea IM, Middleton D. Mitochondrial DNA polymorphism: its role in longevity of the Irish population. Exp Gerontol. 2001; 36:1161–78. https://doi.org/10.1016/S0531-5565(01)00094-8 [PubMed]
  • 4. Hogarth KA, Costford SR, Yoon G, Sondheimer N, Maynes JT. DNM1L Variant Alters Baseline Mitochondrial Function and Response to Stress in a Patient with Severe Neurological Dysfunction. Biochem Genet. 2018; 56:56–77. https://doi.org/10.1007/s10528-017-9829-2 [PubMed]
  • 5. Latorre-Pellicer A, Moreno-Loshuertos R, Lechuga-Vieco AV, Sánchez-Cabo F, Torroja C, Acín-Pérez R, Calvo E, Aix E, González-Guerra A, Logan A, Bernad-Miana ML, Romanos E, Cruz R, et al. Mitochondrial and nuclear DNA matching shapes metabolism and healthy ageing. Nature. 2016; 535:561–65. https://doi.org/10.1038/nature18618 [PubMed]
  • 6. Kim KH, Son JM, Benayoun BA, Lee C. The Mitochondrial-Encoded Peptide MOTS-c Translocates to the Nucleus to Regulate Nuclear Gene Expression in Response to Metabolic Stress. Cell Metab. 2018; 28:516–524.e7. https://doi.org/10.1016/j.cmet.2018.06.008 [PubMed]
  • 7. McManus MJ, Picard M, Chen HW, De Haas HJ, Potluri P, Leipzig J, Towheed A, Angelin A, Sengupta P, Morrow RM, Kauffman BA, Vermulst M, Narula J, Wallace DC. Mitochondrial DNA Variation Dictates Expressivity and Progression of Nuclear DNA Mutations Causing Cardiomyopathy. Cell Metab. 2019; 29:78–90.e5. https://doi.org/10.1016/j.cmet.2018.08.002 [PubMed]
  • 8. Dingley SD, Polyak E, Ostrovsky J, Srinivasan S, Lee I, Rosenfeld AB, Tsukikawa M, Xiao R, Selak MA, Coon JJ, Hebert AS, Grimsrud PA, Kwon YJ, et al. Mitochondrial DNA variant in COX1 subunit significantly alters energy metabolism of geographically divergent wild isolates in Caenorhabditis elegans. J Mol Biol. 2014; 426:2199–216. https://doi.org/10.1016/j.jmb.2014.02.009 [PubMed]
  • 9. Zhu Z, Lu Q, Zeng F, Wang J, Huang S. Compatibility between mitochondrial and nuclear genomes correlates with the quantitative trait of lifespan in Caenorhabditis elegans. Sci Rep. 2015; 5:17303. https://doi.org/10.1038/srep17303 [PubMed]
  • 10. Rockman MV, Kruglyak L. Recombinational landscape and population genomics of Caenorhabditis elegans. PLoS Genet. 2009; 5:e1000419. https://doi.org/10.1371/journal.pgen.1000419 [PubMed]
  • 11. Huang S. New thoughts on an old riddle: what determines genetic diversity within and between species? Genomics. 2016; 108:3–10. https://doi.org/10.1016/j.ygeno.2016.01.008 [PubMed]
  • 12. Viñuela A, Snoek LB, Riksen JA, Kammenga JE. Genome-wide gene expression regulation as a function of genotype and age in C. elegans. Genome Res. 2010; 20:929–37. https://doi.org/10.1101/gr.102160.109 [PubMed]
  • 13. Wang Z, Xu P, Chen B, Zhang Z, Zhang C, Zhan Q, Huang S, Xia ZA, Peng W. Identifying circRNA-associated-ceRNA networks in the hippocampus of Aβ1-42-induced Alzheimer’s disease-like rats using microarray analysis. Aging (Albany NY). 2018; 10:775–88. https://doi.org/10.18632/aging.101427 [PubMed]
  • 14. Wang JD, Zhou HS, Tu XX, He Y, Liu QF, Liu Q, Long ZJ. Prediction of competing endogenous RNA coexpression network as prognostic markers in AML. Aging (Albany NY). 2019; 11:3333–47. https://doi.org/10.18632/aging.101985 [PubMed]
  • 15. Houtkooper RH, Williams RW, Auwerx J. Metabolic networks of longevity. Cell. 2010; 142:9–14. https://doi.org/10.1016/j.cell.2010.06.029 [PubMed]
  • 16. Davies SK, Bundy JG, Leroi AM. Metabolic Youth in Middle Age: Predicting Aging in Caenorhabditis elegans Using Metabolomics. J Proteome Res. 2015; 14:4603–09. https://doi.org/10.1021/acs.jproteome.5b00442 [PubMed]
  • 17. Goyal MS, Blazey TM, Su Y, Couture LE, Durbin TJ, Bateman RJ, Benzinger TL, Morris JC, Raichle ME, Vlassenko AG. Persistent metabolic youth in the aging female brain. Proc Natl Acad Sci USA. 2019; 116:3251–55. https://doi.org/10.1073/pnas.1815917116 [PubMed]
  • 18. Bansal A, Zhu LJ, Yen K, Tissenbaum HA. Uncoupling lifespan and healthspan in Caenorhabditis elegans longevity mutants. Proc Natl Acad Sci USA. 2015; 112:E277–86. https://doi.org/10.1073/pnas.1412192112 [PubMed]
  • 19. Vayndorf EM, Lee SS, Liu RH. Whole apple extracts increase lifespan, healthspan and resistance to stress in Caenorhabditis elegans. J Funct Foods. 2013; 5:1236–43. https://doi.org/10.1016/j.jff.2013.04.006 [PubMed]
  • 20. Andersen EC, Shimko TC, Crissman JR, Ghosh R, Bloom JS, Seidel HS, Gerke JP, Kruglyak L. A Powerful New Quantitative Genetics Platform, Combining Caenorhabditis elegans High-Throughput Fitness Assays with a Large Collection of Recombinant Strains. G3 (Bethesda). 2015; 5:911–20. https://doi.org/10.1534/g3.115.017178 [PubMed]
  • 21. Koopman M, Michels H, Dancy BM, Kamble R, Mouchiroud L, Auwerx J, Nollen EA, Houtkooper RH. A screening-based platform for the assessment of cellular respiration in Caenorhabditis elegans. Nat Protoc. 2016; 11:1798–816. https://doi.org/10.1038/nprot.2016.106 [PubMed]
  • 22. Gerencser AA, Neilson A, Choi SW, Edman U, Yadava N, Oh RJ, Ferrick DA, Nicholls DG, Brand MD. Quantitative microplate-based respirometry with correction for oxygen diffusion. Anal Chem. 2009; 81:6868–78. https://doi.org/10.1021/ac900881z [PubMed]
  • 23. Balaban RS, Nemoto S, Finkel T. Mitochondria, oxidants, and aging. Cell. 2005; 120:483–95. https://doi.org/10.1016/j.cell.2005.02.001 [PubMed]
  • 24. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003; 13:2498–504. https://doi.org/10.1101/gr.1239303 [PubMed]
  • 25. Doroszuk A, Snoek LB, Fradin E, Riksen J, Kammenga J. A genome-wide library of CB4856/N2 introgression lines of Caenorhabditis elegans. Nucleic Acids Res. 2009; 37:e110. https://doi.org/10.1093/nar/gkp528 [PubMed]
  • 26. Fang EF, Scheibye-Knudsen M, Chua KF, Mattson MP, Croteau DL, Bohr VA. Nuclear DNA damage signalling to mitochondria in ageing. Nat Rev Mol Cell Biol. 2016; 17:308–21. https://doi.org/10.1038/nrm.2016.14 [PubMed]
  • 27. Salmena L, Poliseno L, Tay Y, Kats L, Pandolfi PP. A ceRNA hypothesis: the Rosetta Stone of a hidden RNA language? Cell. 2011; 146:353–58. https://doi.org/10.1016/j.cell.2011.07.014 [PubMed]
  • 28. Das A, Gorospe M, Panda AC. The coding potential of circRNAs. Aging (Albany NY). 2018; 10:2228–29. https://doi.org/10.18632/aging.101554 [PubMed]
  • 29. Zheng DW, Xu L, Li CX, Dong X, Pan P, Zhang QL, Li B, Zeng X, Zhang XZ. Photo-Powered Artificial Organelles for ATP Generation and Life-Sustainment. Adv Mater. 2018; 30:e1805038. https://doi.org/10.1002/adma.201805038 [PubMed]
  • 30. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MA, Bender D, Maller J, Sklar P, de Bakker PI, Daly MJ, Sham PC. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. 2007; 81:559–75. https://doi.org/10.1086/519795 [PubMed]
  • 31. Artal-Sanz M, Tavernarakis N. Prohibitin couples diapause signalling to mitochondrial metabolism during ageing in C. elegans. Nature. 2009; 461:793–97. https://doi.org/10.1038/nature08466 [PubMed]
  • 32. Zhu Z, Han X, Wang Y, Liu W, Lu Y, Xu C, Wang X, Hao L, Song Y, Huang S, Rizak JD, Li Y, Han C. Identification of Specific Nuclear Genetic Loci and Genes That Interact With the Mitochondrial Genome and Contribute to Fecundity in Caenorhabditis elegans. Front Genet. 2019; 10:28. https://doi.org/10.3389/fgene.2019.00028 [PubMed]
  • 33. Luz AL, Smith LL, Rooney JP, Meyer JN. Seahorse Xfe 24 Extracellular Flux Analyzer-Based Analysis of Cellular Respiration in Caenorhabditis elegans. Curr Protoc Toxicol. 2015; 66:1–15. https://doi.org/10.1002/0471140856.tx2507s66 [PubMed]
  • 34. Pan Z, Fan Z, Ma J, Liu H, Shen L, He B, Zhang M. Profiling and functional characterization of circulation LncRNAs that are associated with coronary atherosclerotic plaque stability. Am J Transl Res. 2019; 11:3801–15. [PubMed]
  • 35. Hou X, Gong X, Zhang L, Li T, Yuan H, Xie Y, Peng Y, Qiu R, Xia K, Tang B, Jiang H. Identification of a potential exosomal biomarker in spinocerebellar ataxia Type 3/Machado-Joseph disease. Epigenomics. 2019; 11:1037–56. https://doi.org/10.2217/epi-2019-0081 [PubMed]
  • 36. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015; 12:357–60. https://doi.org/10.1038/nmeth.3317 [PubMed]
  • 37. Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015; 33:290–95. https://doi.org/10.1038/nbt.3122 [PubMed]
  • 38. Kong L, Zhang Y, Ye ZQ, Liu XQ, Zhao SQ, Wei L, Gao G. CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. 2007; 35:W345–9. https://doi.org/10.1093/nar/gkm391 [PubMed]
  • 39. Sun L, Luo H, Bu D, Zhao G, Yu K, Zhang C, Liu Y, Chen R, Zhao Y. Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. 2013; 41:e166. https://doi.org/10.1093/nar/gkt646 [PubMed]
  • 40. Finn RD, Coggill P, Eberhardt RY, Eddy SR, Mistry J, Mitchell AL, Potter SC, Punta M, Qureshi M, Sangrador-Vegas A, Salazar GA, Tate J, Bateman A. The Pfam protein families database: towards a more sustainable future. Nucleic Acids Res. 2016; 44:D279–85. https://doi.org/10.1093/nar/gkv1344 [PubMed]
  • 41. Buchfink B, Xie C, Huson DH. Fast and sensitive protein alignment using DIAMOND. Nat Methods. 2015; 12:59–60. https://doi.org/10.1038/nmeth.3176 [PubMed]
  • 42. Conesa A, Götz S, García-Gómez JM, Terol J, Talón M, Robles M. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005; 21:3674–76. https://doi.org/10.1093/bioinformatics/bti610 [PubMed]
  • 43. Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012; 9:357–59. https://doi.org/10.1038/nmeth.1923 [PubMed]
  • 44. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011; 12:323. https://doi.org/10.1186/1471-2105-12-323 [PubMed]
  • 45. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010; 26:139–40. https://doi.org/10.1093/bioinformatics/btp616 [PubMed]
  • 46. Friedländer MR, Chen W, Adamidi C, Maaskola J, Einspanier R, Knespel S, Rajewsky N. Discovering microRNAs from deep sequencing data using miRDeep. Nat Biotechnol. 2008; 26:407–15. https://doi.org/10.1038/nbt1394 [PubMed]
  • 47. ’t Hoen PA, Ariyurek Y, Thygesen HH, Vreugdenhil E, Vossen RH, de Menezes RX, Boer JM, van Ommen GJ, den Dunnen JT. Deep sequencing-based expression analysis shows major advances in robustness, resolution and inter-lab portability over five microarray platforms. Nucleic Acids Res. 2008; 36:e141. https://doi.org/10.1093/nar/gkn705 [PubMed]
  • 48. Kumar L, E Futschik M. Mfuzz: a software package for soft clustering of microarray data. Bioinformation. 2007; 2:5–7. https://doi.org/10.6026/97320630002005 [PubMed]
  • 49. Tafer H, Hofacker IL. RNAplex: a fast tool for RNA-RNA interaction search. Bioinformatics. 2008; 24:2657–63. https://doi.org/10.1093/bioinformatics/btn193 [PubMed]
  • 50. Krüger J, Rehmsmeier M. RNAhybrid: microRNA target prediction easy, fast and flexible. Nucleic Acids Res. 2006; 34:W451–4. https://doi.org/10.1093/nar/gkl243 [PubMed]
  • 51. John B, Enright AJ, Aravin A, Tuschl T, Sander C, Marks DS. Human MicroRNA targets. PLoS Biol. 2004; 2:e363. https://doi.org/10.1371/journal.pbio.0020363 [PubMed]
  • 52. Gao Y, Wang J, Zhao F. CIRI: an efficient and unbiased algorithm for de novo circular RNA identification. Genome Biol. 2015; 16:4. https://doi.org/10.1186/s13059-014-0571-3 [PubMed]
  • 53. Memczak S, Jens M, Elefsinioti A, Torti F, Krueger J, Rybak A, Maier L, Mackowiak SD, Gregersen LH, Munschauer M, Loewer A, Ziebold U, Landthaler M, et al. Circular RNAs are a large class of animal RNAs with regulatory potency. Nature. 2013; 495:333–38. https://doi.org/10.1038/nature11928 [PubMed]
  • 54. Wen B, Mei Z, Zeng C, Liu S. metaX: a flexible and comprehensive software for processing metabolomics data. BMC Bioinformatics. 2017; 18:183. https://doi.org/10.1186/s12859-017-1579-y [PubMed]