A comprehensive evaluation of skin aging-related circular RNA expression profiles.

BACKGROUND: Circular RNAs (circRNAs) have been shown to play important regulatory roles in a range of both pathological and physiological contexts, but their functions in the context of skin aging remain to be clarified. In the present study, we therefore, profiled circRNA expression profiles in four pairs of aged and non-aged skin samples to identify identifying differentially expressed circRNAs that may offer clinical value as biomarkers of the skin aging process.
METHODS: We utilized an RNA-seq to profile the levels of circRNAs in eyelid tissue samples, with qRT-PCR being used to confirm these RNA-seq results, and with bioinformatics approaches being used to predict downstream target miRNAs for differentially expressed circRNAs.
RESULTS: In total, we identified 571 circRNAs with 348 and 223 circRNAs being up and downregulated that were differentially expressed in aged skin samples compared to young skin samples. The top 10 upregulated circRNAs in aged skin sample were hsa_circ_0123543, hsa_circ_0057742, hsa_circ_0088179, hsa_circ_0132428, hsa_circ_0094423, hsa_circ_0008166, hsa_circ_0138184, hsa_circ_0135743, hsa_circ_0114119, and hsa_circ_0131421. The top 10 reduced circRNAs were hsa_circ_0101479, hsa_circ_0003650, hsa_circ_0004249, hsa_circ_0030345, hsa_circ_0047367, hsa_circ_0055629, hsa_circ_0062955, hsa_circ_0005305, hsa_circ_0001627, and hsa_circ_0008531. Functional enrichment analyses revealed the potential functionality of these differentially expressed circRNAs. The top 3 enriched gene ontology (GO) terms of the host genes of differentially expressed circRNAs are regulation of GTPase activity, positive regulation of GTPase activity and autophagy. The top 3 enriched KEGG pathway ID are Lysine degradation, Fatty acid degradation and Inositol phosphate metabolism. The top 3 enriched reactome pathway ID are RAB GEFs exchange GTP for GDP on RABs, Regulation of TP53 Degradation and Regulation of TP53 Expression and Degradation. Six circRNAs were selected for qRT-PCR verification, of which 5 verification results were consistent with the sequencing results. Moreover, targeted miRNAs, such as hsa-miR-588, hsa-miR-612, hsa-miR-4487, hsa-miR-149-5p, hsa-miR-494-5p were predicted for circRna-miRna interaction networks.
CONCLUSION: Overall, these results offer new insights into circRNA expression profiles, potentially highlighting future avenues for research regarding the roles of these circRNAs in the context of skin aging.

INTRODUCTION

Aging is a complex process wherein cells and organisms undergo progressive changes at the molecular, tissue, and organ levels as a result of either normal physiology or pathological conditions. The skin is the largest organ in the body, serving as a barrier against external threats. Skin aging results in continuous alterations in the functionality and appearance of the skin and can occur as either a result of internal factors or due to extrinsic photoaging. , Photoaging is believed to account for approximately 80% of skin aging. Even so, the mechanisms that govern the skin aging process remain poorly clarified, and few aging‐related biomarkers have been identified to date. Transcriptomic analyses offer a powerful approach to identifying molecular biomarkers of skin aging. , Recent breakthroughs in the development of high‐throughput transcriptome sequencing technologies have led to the discovery of a diverse array of biomarkers of different human diseases, with non‐coding RNAs being commonly studied in this context. ,

Research conducted in the 1970s identified novel non‐coding circular RNAs (circRNAs) in the context of plant viral infections, although at the time these RNAs were believed to lack functional relevance and were instead thought to be a result of unusual splicing reactions. More recent work, however, suggests that circRNAs are key regulators of gene expression at the transcriptional and post‐transcriptional levels. , Indeed, some circRNAs are believed to function as competing endogenous RNAs (ceRNAs) or as de facto molecular sponges capable of sequestering and altering the functionality or expression of specific miRNAs or proteins. , While they generally lack coding potential, there is also some evidence that certain circRNAs may be translated in some contexts, thereby further modulating biological functionality. , , , , , , Indeed, circRNA dysregulation is a hallmark of conditions such as cancer, neurological disease, and cardiovascular diseases. , The functional importance of circRNAs in the context of skin aging, however, remains to be clarified.

Herein, we evaluated circRNA expression profiles in samples of aged and non‐aged skin tissue via high‐throughput sequencing to identify differentially expressed circRNAs (DECs). Appropriate bioinformatics analyses were then used to further predict the functional roles of these DECs in the aging process.

MATERIALS AND METHODS

Sample collection

The Institutional Review Board of The First Hospital of China Medical University approved the present study, with all participants having provided written informed consent. In total, we collected eyelid tissue samples from 28 females undergoing double eyelid surgery at The First Hospital of China Medical University between 2018 and 2019. Fourteen young patient samples were collected from individuals aged 17–23 years, whereas 14 aged patient samples were from those 55–70 years old. We ultimately used eight samples for RNA‐seq analyses (4 young, 4 aged), while 40 samples were used for downstream qRT‐PCR verification of our results.

High‐throughput sequencing.

A Hipure Total RNA Mini Kit (Magen) was used to extract RNA from these samples to the protocol. A Qubit 3.0 Fluorometer (Invitrogen), and Agilent 2100 Bioanalyzer (Applied Biosystems) were then used to assess RNA concentrations and integrity. Only samples yielding a RIN value of ≥7.0 were used for downstream RNA‐sequencing.

A total of 1 μg of RNA per sample was used together with a KAPA RNA HyperPrep Kit with RiboErase (HMR) for Illumina® (Kapa Biosystems, Inc.) to eliminate rRNA prior to library preparation. Samples were then treated for 30 minutes at 37℃ with 10U RNase R (Geneseed).

We next fragmented the remaining RNA, after which first‐ and second‐strand synthesis reactions were conducted. Tails and adapters were then ligated to purified cDNA samples, and amplification of the adapter‐ligated purified DNA was performed. A DNA 1000 chip was then used to evaluate library quality with an Agilent 2100 Bioanalyzer. A qRT‐PCR‐based KAPA Biosystems Library Quantification kit (Kapa Biosystems, Inc.) was used to accurately quantify prepared samples, after which libraries were diluted to a 10 nM concentration and pooled in equimolar amounts. We then conducted 150 bp paired‐end (PE150) sequencing of all samples.

Bioinformatics analysis

Initially, reads were mapped to the latest UCSC transcript set with Bowtie2 v2.1.0, after which RSEM v1.2.15 was used to estimate gene expression levels. Gene expression was normalized via a TMM (trimmed mean of M‐values) approach, with edgeR being used to identify differentially expressed genes. Genes were considered to be differentially expressed if they met the following criteria: p < 0.05 and >1.5 fold change.

To evaluate circRNA expression, STAR was used to map reads to the genome after which DCC was employed to evaluate circRNA expression levels. DECs were identified using edgeR, with resultant figures being generated using appropriate R packages.

In order to explore the potential functional relevance of identified DECs, we analyzed DEC target genes via gene ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Reactome functional enrichment analyses.

circRNA‐miRNA interaction network construction

To predict interactions between DECs and target miRNAs. We obtained miRNA sequences from the miRBase database with corresponding annotations, after which Miranda 3.3a was used to calculate binding interactions between miRNAs and circRNAs. We then generated a visualized version of the resultant DEC‐miRNA interaction network using Cytoscape (v3.7.2; Institute of Systems Biology).

Quantitative real‐time PCR

A Hipure Total RNA Mini Kit (Magen) was used to isolate RNA samples as above, after which a NanoDrop™ One (Thermo Fisher Scientific) instrument was used to assess RNA quality and quantity. In addition, 1.5% agarose gel electrophoresis was used to evaluate RNA integrity and to assess for the presence of any contaminating gDNA, while spectrophotometry at 260–280 nm was employed to assess the purity of these RNA samples. Next, 1 µg of total RNA per sample was used with a PrimeScript RT reagent kit (Takara Bio) to prepare cDNA samples. Primer3 (v. 0.4.0) was used to design all primers in the present study (Table 1), with GAPDH being used for normalization purposes. A 2−△△CT approach was used to confirm differential circRNA expression in analyzed samples, which were assessed in triplicate.

TABLE 1

The primers used for RT‐qPCR

Gene name	Primer sequences (5′ to 3′)	Product length (bp)
GAPDH	F: AGAAGGCTGGGGCTCATTTG R: GCAGGAGGCATTGCTGATGAT	140
hsa_circ_0137613	F: GCTGATGTCATTCTCCACAAGG R: AGCAGCAGCTGACACAGGAT	212
hsa_circ_0077605	F: CTGGAGCACATGAGCCTGCA R: TGAAAGGTCGCTCCCCTGTGT	208
hsa_circ_0000205	F: GTTGCTGACACTAGTCTTATTG R: CGAGCTGTTAGTTCTTCGTA	148
hsa_circ_0003803	F: CCTGGGCGCACAGAAAATCC R: CACCTCTCGGAGTTTCCTCTG	160
hsa_circ_0113488	F: GCTCATCAAAGACATTTATATGATA R: GAAGAAATTGTAGGCTGTTC	173
hsa_circ_0112861	F: ACATCAGTGGAGAACCTCAGT R: GACAGTGTTGGTCTTCCATTCA	166

Statistical analysis

Data are means ± standard deviation (SD). GraphPad Prism v7.0 was used for all statistical testing, and data were compared via Student's t‐tests with p < 0.05 as the significance threshold.

RESULTS

Detection of circRNA dysregulation in the context of skin aging

In order to identify circRNAs associated with the skin aging processes, we began by employing an RNA‐seq analysis approach to analyze four eyelid skin samples from young adults and four eyelid skin samples from aged adults, with young samples being treated as controls. In total, we identified 14915 circRNAs, of which 571 were found to be differentially expressed (Fold change >1.5 and p < 0.05) between these groups (Figure 1A). Of these DECs, 348 and 223 were found to be up and downregulated, respectively (Figure 1B). These DECs were distributed among all chromosomes, shown as a circos plots (Figure 1C). The details regarding the top 10 upregulated and 10 downregulated circRNAs are presented in Table 2

FIGURE 1

Profiling of skin aging‐related circRNAs. A, Those circRNAs that were up or downregulated (red and green, respectively) in aged skin tissue samples relative to young tissues samples were arranged in a heat map, with rows corresponding to individual circRNAs and columns corresponding to individual samples. B, Differentially regulated circRNAs were arranged in a Volcano plot, with blue dots corresponding to a lack of statistical significance, whereas red and green correspond to up and downregulated circRNAs, respectively (Fold change >1.5 and p < 0.05). C, Locations of DECs on human chromosomes are represented by circos plots, with chromosomes being represented by the outermost circle, whereas circRNAs are indicated in the middle circle, and DECs are shown in the innermost circle. Up and downregulated circRNAs are represented by red and blue lines, respectively

TABLE 2

Top 10 upregulated and 10 downregulated circRNAs in aging skin samples ranked by fold changes

circRNA ID	Host gene	Expression	Log 2 FC	p value	Chr
hsa_circ_0123543	UBP1	Up	6.61068	0.00098	3
hsa_circ_0057742	FAM126B	Up	6.57177	0.00149	2
hsa_circ_0088179	AKNA	Up	6.52851	0.00129	9
hsa_circ_0132428	DOPEY1	Up	6.52851	0.00112	6
hsa_circ_0094423	BTAF1	Up	6.36737	0.00121	10
hsa_circ_0008166	TDG	Up	6.14096	0.00052	12
hsa_circ_0138184	SETX	Up	6.09941	0.00345	9
hsa_circ_0135743	ASAP1	Up	5.99863	0.00387	8
hsa_circ_0114119	ACADM	Up	5.81751	0.00279	1
hsa_circ_0131421	KDM1B	Up	5.68353	0.00276	6
hsa_circ_0101479	ACIN1	Down	−6.99402	0.00017	14
hsa_circ_0003650	KPNB1	Down	−6.88085	0.00082	17
hsa_circ_0004249	PKM	Down	−6.51587	0.00079	15
hsa_circ_0030345	VPS36	Down	−6.47233	0.00121	13
hsa_circ_0047367	TRAPPC8	Down	−6.29717	0.00234	18
hsa_circ_0055629	STARD7	Down	−6.05903	0.00166	2
hsa_circ_0062955	DEPDC5	Down	−5.85501	0.00365	22
hsa_circ_0005305	VPS33A	Down	−5.67488	0.00231	12
hsa_circ_0001627	BACH2	Down	−5.16169	0.00329	6
hsa_circ_0008531	ARID4B	Down	−5.02842	0.00325	1

Functional enrichment analyses

We next conducted GO, KEGG, and Reactome analyses to explore the potential functional roles of identified DECs in the skin aging process. GO analyses suggested that these circRNAs are most closely associated with specific cellular components, molecular functions, and biological processes (Figure 2). Top enriched KEGG pathways associated with these DECs included the lysine degradation, inositol phosphate metabolism, and purine metabolism (Figure 3A). Top Reactome pathways for these DECs included SUMO E3 ligase SUMOylation of target proteins, Rab regulation of trafficking, and the Rho GTPase cycle (Figure 3B).

FIGURE 2

GO analysis of DEC host genes associated with aging

FIGURE 3

KEGG and Reactome pathway analyses in aged and young skin sample groups. A, The top 15 enriched KEGG pathways. B, The top 15 enriched Reactome pathways

Confirmation of differential circRNA expression

In order to confirm the validity of our RNA‐seq results, three random upregulated circRNAs (hsa_circ_0137613, hsa_circ_0077605, hsa_circ_0000205) and three random downregulated circRNAs (hsa_circ_0003803, hsa_circ_0113488, hsa_circ_0112861) were selected for qRT‐PCR analysis. For these analyses, we assessed circRNA expression levels in 40 total eyelid skin tissue samples (20 age, 20 young) using appropriate divergent primers. For full details regarding selected circRNAs, see Table 3.

TABLE 3

Six circRNAs were selected to perform further PCR validation

Circbase ID	Host gene	Fold change	p value	Expression	Chr	Length (nt)
hsa_circ_0137613	ABCA1	4.90	0.000686	Up	9	654
hsa_circ_0077605	ZBTB24	5.85	0.004417	Up	6	336
hsa_circ_0000205	WDR37	4.45	0.019171	Up	10	766
hsa_circ_0003803	UTRN	−5.08	0.006988	Down	6	459
hsa_circ_0113488	SPATA6	−3.65	0.015928	Down	1	708
hsa_circ_0112861	SCCPDH	−8.53	0.007961	Down	1	324

We found that all the six DECs exhibited identical trends upon qRT‐PCR analysis to those observed in our RNA‐seq results (Figure 4). Of these, five circRNAs (hsa_circ_0137613, hsa_circ_0077605, hsa_circ_0003803, hsa_circ_0113488, hsa_circ_0112861) exhibited significant differential expression between control and aging skin group (p < 0.05). We additionally utilized Sanger sequencing to verify the qRT‐PCR products, confirming them to be consistent with the sequences in the circbase (http://www.circbase.org/) database (Supplemental Figure S1).

FIGURE 4

qRT‐PCR validation of 6 differentially expressed circRNAs in 20 young and 20 aged skin samples

circRNA‐miRNA interaction network generation

To further evaluate the functional roles of these identified circRNAs, we assessed the potential miRNA binding of six circRNAs using miRnada (v 3.3a), with the resultant interaction network being visualized using the Cytoscape (V3.6.0) tool (Figure 5). In total, this network incorporated 6 circRNAs and 151 target miRNAs. Among these 151 miRNAs, we found that 6 miRNAs, hsa‐miR‐3064‐5p, hsa‐miR‐6762‐3p, hsa‐miR‐3194‐3p, hsa‐miR‐6731‐5p, hsa‐miR‐30c‐1‐3p and hsa‐miR‐6760‐5p not only bind to 1 circRNA, but they may be also regulated by multiple circRNAs at the same time.

FIGURE 5

Predicted circRNA‐miRNA interaction networks. Yellow and red rectangles correspond to circRNAs and target miRNAs, respectively

DISCUSSION

Skin aging is an ongoing process influenced by both intrinsic and extrinsic factors that ultimately leads to impaired tissue integrity and functionality. As skin ages, it loses its elasticity and undergoes thickening, drying, and wrinkling. Increased ultraviolet radiation exposure as a result of ozone layer degradation has further expedited the skin aging process in some areas of the world. Effective approaches to preventing skin aging, however, remain to be established and are of great interest to the pharmaceutical and cosmetic industries. At present, the molecular mechanisms governing skin aging remain poorly understood, and few relevant biomarkers of this progress have been identified.

Owing to their closed‐loop structures, circRNAs are more resistant to exonuclease‐mediated degradation than are linear RNAs, making them ideal as biomarkers of certain diseases. , , , However, relatively few studies have assessed the relationship between circRNA expression patterns and the skin aging process. Amaresh et al. evaluated circRNA expression profiles in WI‐38 cells from early and late passages via RNA‐seq and identified circPVT1 as a suppressor of cellular senescence that sequesters let7. Peng et al. similarly identified 29 circRNAs that were differentially expressed in human dermal fibroblasts following UVA irradiation and determined that circCOL3A1 serves as a regulator of the expression of type I collagen during the photoaging process. Si et al. also generated a UVB‐induced model of fibroblast senescence in which they identified 472 DECs via microarray. Of these, they identified circRNA_100797 acts as a molecular sponge capable of sequestering miR‐23a‐5p and inhibiting the UVB‐induced photoaging of these cells.

Unlike these prior studies analyzing cell models, in the present study we assessed differential circRNA expression in primary eyelid tissue samples from young and aged human donors. We ultimately identified 571 DECs and 5 circRNA expression levels in the two groups of samples were confirmed to have significant differences via qPCR. Hsa_circ_0137613 and has_circ_0077605 are significantly upregulated in the aging skin group, suggesting that they may have the biological functions of promoting skin cell apoptosis, reducing cell proliferation or arresting the cell cycle in G0/G1 phase. On the contrary, the expression of hsa_circ_0003803, hsa_circ_0113488 and hsa_circ_0112861 in aging skin tissues was significantly reduced, suggesting that they may have the function of accelerating skin cell proliferation and inhibiting skin cell apoptosis. We will use gain‐of‐function and loss‐of‐function experiments to verify the functions of these molecules in skin cell models. Moreover, we conducted GO, KEGG, and Reactome functional enrichment analyses to gain insight into the molecular mechanisms whereby these circRNAs functioned in this context. The interactions between circRNAs and miRNAs in the skin aging process are also poorly understood, and as such we utilized the miRnada program to construct a putative circRNA‐miRNA interaction network.

There are multiple limitations to the present study. For one, our sample size was relatively limited and should be expanded in future analyses. Furthermore, both in vivo and in vitro validation of the functional relevance of these DECs is still required. The identified circRNAs have the potential to function through a range of mechanisms including as ceRNAs or as translated peptides/proteins to modulate the aging process, and as such our analyses are too cursory to offer meaningful insights into their roles in this context. Despite these limitations, our results represent a novel dataset pertaining to the identification of circRNAs associated with the skin aging process, with hsa_circ_0137613, hsa_circ_0077605, hsa_circ_0003803, hsa_circ_0113488 and hsa_circ_0112861 being potential skin aging‐related biomarkers.

CONFLICT OF INTEREST

The authors declare they have no competing interests.

Fig S1

Click here for additional data file.