Long-Term Exercise Alters the Profiles of Circulating Micro-RNAs in the Plasma of Young Women

Objective: The objective of this paper was to study the effects of long-term exercise on circulating microRNAs (miRNAs) in human plasma. Methods: Whole blood was collected from 10 female elite athletes with at least 5 years of training experience in a Synchronized Swimming Group (S group) and 15 female college students without regular exercise training (C group). Plasma miRNAs were then isolated, sequenced, and semi-quantified by the second-generation sequencing technology, and the results were analyzed by bioinformatics methods. Results: We found 380 differentially expressed miRNAs in the S group compared with the C group, among which 238 miRNAs were upregulated and 142 were downregulated. The top five abundant miRNAs in the 380 miRNAs of the S group are hsa-miR-451a, hsa-miR-486, hsa-miR-21-5p, hsa-miR-423-5p, and hsa-let-7b-5p. Muscle-specific/enriched miRNAs were not significantly different, except for miR-206 and miR-486. According to the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, a large proportion of the differentially expressed miRNAs are targeted in cancer-related pathways, including proteoglycans in cancer and miRNAs in cancer and basal cell carcinoma. As the levels of circulating miRNAs (ci-miRNAs) are commonly known to be significantly deregulated in cancer patients, we further compared the levels of some well-studied miRNAs in different types of cancer patients with those in the S group and found that long-term exercise regulates the level of ci-miRNAs in an opposite direction to those in cancer patients. Conclusion: Long-term exercise significantly alters the profiles of plasma miRNAs in healthy young women. It may reduce the risk of certain types of cancers by regulating plasma miRNA levels.

A few studies have shown that plasma ci-miRNAs can also be affected by a single bout of acute exercise (Baggish et al., 2011;Aoi et al., 2013;Nielsen et al., 2014;Cui et al., 2016;Shah et al., 2017;Barber et al., 2019) or sustained training from 28 to 140 days (Baggish et al., 2011;Aoi et al., 2013;Nielsen et al., 2014;Li et al., 2018;Barber et al., 2019). Most of them focused on some specific ci-miRNAs, such as ci-miRNAs involved in inflammation, skeletal and cardiac muscle contractility, and hypoxia/ischemia adaptation (Baggish et al., 2011), muscle-enriched miRNAs (Aoi et al., 2013), miRNAs related to cardiovascular disease (Barber et al., 2019), and those involved in angiogenesis and inflammation and are enriched in muscle and/or cardiac tissues (Li et al., 2018). Only one study  quantified the expression levels of 188 known ci-miRNAs in the plasma of 32 healthy, trained men before and after 12 weeks of endurance exercise and found nine miRNAs expressed differently. By far, the response of whole ci-miRNAs in plasma to long-term exercise (>5 years), especially among long-term trained professional athletes, remains unknown. In this study, we chose athletes of a Synchronized Swimming Group as our subjects since the movement difference between players is small and the team members usually have years of regular exercise experience. We compared the total plasma miRNAs of 10 players from a Synchronized Swimming Group and 15 female college students without regular exercise experience, using second-generation sequencing technology, to study the impact of long-term exercise on the plasma miRNA content.

Subjects
The participants were 10 Synchronized Swimming Group members and 15 female college students of Guangdong Pharmaceutical University. All participants were informed of the methods, procedures, and risks, and then they signed an informed consent document, which was approved by the Ethics Committee of South China Normal University School of Sports Science (approval number 2017102001).

MiRNA Isolation From Plasma
Blood samples were collected around 7:00 a.m., and participants did not have breakfast or vigorous exercise for 24 h. Whole blood (10 ml) from subjects was collected via a direct venous puncture into tubes with ethylenediaminetetraacetic acid (EDTA) as an anticoagulant. Blood was put on ice for half an hour and then separated by 2,000 rpm for 10 min at 4 • C. Then, the plasma was transferred to an RNase-free tube and stored at −80 • C. The RNA samples were extracted with a TRIzol reagent (Life Technologies, cat. no. 15596018). The purity and concentration of the RNA samples were assessed with Nanodrop and Qubit 2.0. RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, United States) to ensure the use of qualified samples for sequencing.

MiRNA Sequencing
RNA samples undergo a series of strict quality control. Qualified samples are used for library construction. Library construction is strictly in accordance with NEBNext Ultra Small RNA Sample Library Prep Kit for Illumina, qualified libraries for highthroughput sequencing. The sequencing platform was Illumina HiSeq X Ten, and read length was single-end (SE) 50 nt.

Quality Control
Raw data (raw reads) of fastq format were firstly processed through in-house perl scripts. In this step, clean data (clean reads) were obtained by removing reads containing adapter or ploy-N and low-quality reads from raw data. The reads were trimmed and cleaned by removing the sequences smaller than 15 nt or longer than 35 nt. At the same time, Q20, Q30, GC content, and the sequence duplication level of the clean data were calculated. All the downstream analyses were based on clean data with high quality.

Comparative Analysis
Using Bowtie software, analysis of the clean reads was done, respectively, with the Silva database, GtRNAdb database, Rfam database, and Repbase database for sequence alignment, filter ribosomal RNA (rRNA), transfer RNA (tRNA), small nuclear RNA (snRNA), small nucleolar RNA (snoRNA), and other ncRNAs and repeats. The remaining reads were used to detect known miRNAs and novel miRNAs predicted by comparing with Genome and known miRNAs from miRBase. Randfold software was used for novel miRNA secondary structure prediction.

Differential Expression Analysis
The expression of miRNAs in each sample was counted and normalized by the Trusted Platform Module (TPM) algorithm (Fahlgren et al., 2007). The normalization formula for TPM is as follows: TPM = Actual miRNA read count/Total miRNA read count × 1,000,000. Differential expression analysis was performed using the DESeq R package (1.10.1). The resulting P values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate (FDR). MiRNAs with |log 2 (FC, fold change)| ≥ 1 and FDR < 0.05 found by DESeq were assigned as differentially expressed miRNAs.

Target Gene Functional Annotation
Target gene prediction was done using miRanda (Doron Betel et al., 2008) and RNAhybrid (Rehmsmeier et al., 2004). Gene function was annotated based on the following databases: Nr (NCBI non-redundant protein sequences), Pfam (Protein family), KOG/COG (Clusters of Orthologous Groups of proteins), Swiss-Prot (a manually annotated and reviewed protein sequence database), KEGG (KEGG Ortholog database), and GO (Gene Ontology).

Enrichment Analysis
A GO term enrichment analysis 1 was used to identify the biological processes, molecular functions, and cellular components associated with the target genes of differential expression miRNAs. KEGG pathway enrichment analysis 2 was used to identify the miRNA targets associated with signaling pathways in the S vs. the C group.

MiRNAs as Biomarkers of Certain Diseases (Study Selection Criteria)
We screened ci-miRNA biomarkers related to certain types of cancers through the Human MicroRNA Disease Database (HMDD 3 ). Studies were included if they (i) were original research and (ii) evaluated the ci-miRNA levels in a specific disease. Studies were excluded if they were published in a language other than English and are review articles, book chapters, conference abstracts, editorials/commentaries/expert opinion, theses, or dissertations.

The Expression Patterns of Ci-MiRNAs Were Highly Similar Within the Originally Divided Participants
The participants included 10 Synchronized Swimming Group members (S group) and 15 female college students from Guangdong Pharmaceutical University in China (C group). The average age of these 10 players was 19.1 years compared with 19.4 years of the control, and there was no significant age difference between these two groups (19.1 ± 2.08 vs. 19.4 ± 0.83 years, P = 0.67). The subjects were all healthy and without known diseases. C group members had no exercise habit and did not enroll in any long-term sports or dancing in the last 5 years. For details, see Supplementary Table 1. Pearson's correlation coefficient was used to test the similarity of these samples at a global level (Figure 1). The heat map of the correlations between each sample was generated and showed distinct expression pattern between the two groups. The clustering of samples showed that the same group of samples is very similar to the groups of samples originally divided.

Differentially Expressed MiRNAs Between the S and C Groups
Twenty-five specimens were successfully sequenced by using Illumina Genome Analyzer II. The number of raw reads obtained per specimen ranged from 15,872,436 to 32,103,564 (average = 23,148,465). From the raw reads, an average of 22,508,866 clean reads (ranging from 15,569,728 to 31,145,335) were filtered and mapped to human genome. After analyzing and normalizing all the mapped reads in both groups, we identified the existence of 3,531 miRNAs with 1,733 known miRNAs and 1,798 novel miRNAs. Then, we applied a stringent filtering criterion to differentiate S and C (FDR < 0.05, | log 2 (FC)| ≥ 1) and identified 380 miRNAs that were differentially expressed (for details, see Supplementary Table 2). Specifically, 238 miRNAs were upregulated and 142 were downregulated (Figure 2A). In addition, 23 miRNAs have an extremely significant difference between the two groups (| log 2 (FC)| ≥ 5) ( Figure 2B). Notably, all 23 miRNAs were upregulated in the S group, and the highest differentially expressed miRNA was novel-miR-1150 (| log 2 (FC)| = 8.09). Among the 380 differentially expressed miRNAs, the five most abundant in the S group were miR-451a, miR-486, miR-423-5p, let-7b-5p, and miR-21-5p, with | log 2 (FC)| = 4.77, 4.76, 2.13, 2.65, and 1. 53 and TPM = 224,122,119,250,85,314,43,867,and 26,667, respectively. On the other hand, the five most abundant miRNAs in the C group were miR-21-5p, let-7f-5p, miR-148a-3p, miR-423-5p, and miR-146a-5p. Taken together, there are consistent and significant differences in the plasma miRNA profiles between the athlete group and the normal college student group.

Target Gene Prediction of Differentially Expressed MiRNAs
To explore the functions of chronic exercise induced differentially expressed miRNAs, the target genes of these miRNAs were predicted by combining miRanda and RNAhybrid to reduce the probability of false positives. From this analysis, 20,449 target genes corresponding to 3,531 miRNAs were identified. For details, see Supplementary Table 3.

Functional Annotation Analysis
Gene ontology enrichment analysis of the 380 differentially expressed genes (DEGs) was implemented by the GOseq R packages based on Wallenius' non-central hyper-geometric distribution. In the S group, 1,162 GO categories of biological processes, 142 of cellular component, and 203 of molecular function were found to be significantly affected (Kolmogorov-Smirnov, KS < 0.05). The 10 most significantly enriched GO categories of biological processes, cellular component, and molecular function are listed in Table 1.

Pathway Annotation Analysis
Pathway annotation analysis of the differentially expressed miRNA target genes was done using the KEGG database see text footnote 2). The results are shown in Figure 3.  GO.ID, ID of the GO term; Term, GO function description; Annotated, the number of genes annotated with this function; Significant, the number of differentially expressed miRNA target genes annotated to this function; Expected, the expected number of differentially expressed miRNA target genes annotated to this function; KS statistical significance of the enrichment term (the smaller the KS value, the more significant the enrichment).

Circulating MiRNAs as Biomarkers of Different Diseases
We screened ci-miRNA biomarkers related to certain types of diseases through the HMDD. The results are shown in Table 2.

Differentially Expressed MiRNAs in Athletes Are Likely Due to Long-Term Exercise
As shown in Figure 1, the expression profiles of the ci-miRNAs in athletes were significantly different from those of college students who had not received any regular sports training. Notably, all the participants from both groups in this study come from different parts of China, eating different food. In addition, there are no immediate training seasons prior to the blood withdrawal in the athlete group. In addition to synchronized swimming, the athletes also underwent strength training and aerobic and anaerobic training. Thus, the differentially expressed miRNAs in the athlete group may be a result of combined exercise training. Altogether, the ci-miRNAs' profile changes in the athlete group are highly likely due to long-term exercise.

Ci-MiRNA Profile Changes in Athletes Are Correlated With Muscle-Enriched MiRNAs
A set of miRNAs, called myo-miRs, have been identified in skeletal muscle and/or the myocardium, including miR-1, miR-133a, miR-133b, miR-206, miR-208a, miR-208b, and miR-499 (McCarthy and Esser, 2007;Callis et al., 2008;van Rooij et al., 2008;Small et al., 2010). Their expression is at least 20-fold FIGURE 3 | KEGG pathway enrichment bubble plot of the differentially expressed miRNAs. X-axis label represents rich factor (rich factor = amount of differentially expressed genes enriched in the pathway/amount of all genes in the background gene set) and Y -axis label represents pathway. The size and color of the bubble represent the amount of differentially expressed genes enriched in the pathway and the enrichment significance (Q_ value is the P value corrected after multiple hypothesis testing), respectively.
higher than the mean expression value of the other tissues (Lee et al., 2008). Among them, miR-206 is skeletal musclespecific (McCarthy, 2008) and miR-208a is cardiac musclespecific (Sempere et al., 2004). MiR-486 is also a skeletal muscle and cardiac-enriched miRNA, which was included in myo-miRs in later studies (Horak et al., 2016). In our study, miR-1 was not detected in all samples. the TPMs of the remaining myo-miRs are shown in Table 3.
To test whether the level changes of the other differently expressed miRNAs were associated with musclespecific/enriched miRNAs, we first calculated a Pearson's correlation coefficient between miR-486 and the other differently expressed miRNAs. As a myo-miR, miR-486 mainly functions as a regulator of myoblast proliferation and migration (Alexander et al., 2011) and skeletal muscle size (Hitachi et al., 2014). The results showed that a large portion of the differently expressed miRNAs (>30%, R 2 > 0.6) were highly correlated with miR-486-5p and miR-486-3p (for details, see Supplementary Table 4), which provide another evidence that the miRNA profile changes in the athlete group were likely to be induced by exercise.
Muscle injures are common among athletes. In the plasma of healthy young men, the expression level of miR-206 was significantly higher after high-intensity intermittent exercise (HIIE) or vigorous-intensity continuous exercise (VICE) (Cui et al., 2016). miR-206 is also increased in the serum of mdx and CXMDJ mice (Mizuno et al., 2011) or in response to cardiotoxin (CTX)-induced injury, but markedly decreased from day 3 to day 5 after CTX injury (Matsuzaka et al., 2014). This suggests that high levels of miR-206 in plasma may be associated with skeletal muscle injury. In addition, miR-133 has also been demonstrated to be associated with skeletal muscle injury (Laterza et al., 2009). In the present study, miR-206 was slightly increased in the S group (|log 2 FC| = 1.99, TPM = 0-42.44), and miR-133 did not show significant difference between the two groups. To examine whether the level changes of miRNAs in the athlete group were also associated with muscle injury, we also calculated Pearson's correlation coefficient between miR-206 and the other differently expressed miRNAs and found that there were less than 3% of the differentially expressed miRNAs correlated with muscle injury (<3%, R 2 > 0.6; for details, see Supplementary Table 4).
vigorous-intensity physical activity can decrease the risk of breast cancer, and physically active individuals have a lower risk of lung cancer (Lee, 2003). The overall cancer incidence was lower in athletes than in the general population (Robsahm et al., 2010;Sormunen et al., 2014). And former college athletes had a significantly lower risk of breast cancer than do the nonathletes (Wyshak and Frisch, 2000). However, the underlying mechanism regarding exercise reducing cancer risks remains to  Table 2) can be used as biomarkers for cancer or other diseases (Komatsu et al., 2014;Zhang et al., 2016;Hara et al., 2017;Khalighfard et al., 2018). Strikingly, long-term exercise regulates the levels of these ci-miRNAs in an opposite direction to those in cancer patients, which suggests that long-term exercise may reduce the risk of cancers through regulating ci-miRNAs. Interestingly, some of the elevated miRNAs, such as miR-126 and miR-423, in heart failure and coronary artery disease (CAD) were also upregulated in the athlete group. Coronary artery calcification (CAC) is a strong predictor of incident CAD (Detrano et al., 2008). In a prospective study, 27% higher odds of CAC have been seen among participants who exceeded physical activity guidelines vs. those below physical activity guidelines (Laddu et al., 2017). CAD is also a leading cause of sudden cardiac death in athletes over 35 years of age (Morrison et al., 2018). But the incidence rates of sports-related sudden cardiac deaths in noncompetitive and competitive athletes are not different (Risgaard et al., 2014). Does high-intensity or volume exercise increase the risk of CAD? Is miR-126 not a solid biomarker of CAD? Further research is needed to address these questions.

CONSLUSION
Long-term exercise significantly alters the plasma miRNA profiles in healthy young women, which may reduce the risk of lung, breast, pancreas, melanoma, cholangiocarcinoma, esophageal neoplasms, and colorectal carcinoma. Our current study described an overall effect of long-term exercise on circulating miRNA profiles in the plasma of young women. Further studies are needed to determine whether different types of exercises have different effects on ci-miRNAs and whether there are gender and age differences.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Ethics Committee of South China Normal University School of Sports Science (Approval Number 2017102001). Written informed consent to participate in this study was provided by the participants.

AUTHOR CONTRIBUTIONS
RD and CL initiated and supported the project. FL and MB performed the experiments. FL and RD wrote the manuscript. All authors analyzed the data and commented on the manuscript.

FUNDING
This work was supported by the NSFC (31971096 and 31771256) and Guangdong Science and Technology Collaborative Innovation Center for Sports Science (2019B110210004).