Comprehensive RNA Sequencing in Adenoma-Cancer Transition Identified Predictive Biomarkers and Therapeutic Targets of Human CRC.

Specific molecular biomarkers for predicting the transition from colorectal adenoma to cancer have been identified, however, circular RNA (circRNA)-related signatures remain to be clarified. We carried out high-throughput RNA sequencing to determine the expression profiles of circRNAs, microRNAs (miRNAs), and mRNAs in human colorectal cancer (CRC), adenoma, and adjacent normal tissues. We identified 84 circRNAs, 41 miRNAs, and 398 mRNAs that were commonly differentially expressed in CRC and adenoma tissues compared with normal tissues. Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, and protein-protein interaction (PPI) analyses identified numerous cancer-related hub genes that might serve as potential therapeutic targets in CRC. Competing endogenous RNA (ceRNA) networks, including three circRNAs, three miRNAs, and 28 mRNAs were constructed, suggesting their potential role in cancer progression. Representative differentially expressed RNAs were validated by the Cancer Genome Atlas (TCGA) database and real-time PCR experiments. Receiver operating characteristic (ROC) curve analysis identified three circRNAs (hsa_circ_0049487, hsa_circ_0066875, and hsa_circ_0007444) as possible novel biomarkers predicting the transition from colonic adenoma to cancer. Overall, our findings may provide novel perspectives to clarify the mechanisms of the transition from premalignant adenoma to cancer and identify specific circRNA-related signatures with possible applications for the early diagnosis of and as potential therapeutic targets in CRC.

Introduction

Colorectal cancer (CRC) is the third most common cancer worldwide, with increasing incidence and mortality in China., The transition from adenoma to carcinoma has been proposed as a critical step in colorectal carcinogenesis, including multiple genetic and epigenetic risk factors., Up to 90% of CRC cases follow this transition over a period of many years, offering opportunities for the early screening, prediction, and prevention of CRC.

Genetic and epigenetic alterations are critical in driving adenoma-carcinoma progression, and numerous molecular signatures have been identified. For example, EZH2 and COX-2 gene expression levels were increased in both colorectal adenoma and carcinoma compared with normal tissues, contributing to CRC progression and with potential uses in the early diagnosis of CRC., Whole-genome expression profiling of colorectal adenoma and CRC samples revealed a set of downregulated transcripts (e.g., MAL, SFRP1, SULT1A1, etc.) along the adenoma-carcinoma transition. Researchers also identified specific microRNA (miRNA) profiles in the adenoma-to-carcinoma progression using microarray analysis.

Circular RNAs (circRNAs) are a novel class of endogenous noncoding RNAs that form covalently closed continuous loops without 3′ end poly(A) tails and 5′ end caps. circRNAs regulate gene expression at the transcriptional or post-transcriptional level by acting as miRNA sponges (competing endogenous RNAs, ceRNAs) or binding to other molecules. Accumulating evidence has suggested that circRNAs play an important role in cancer progression. circLMTK2 expression was reduced in gastric cancer, suggesting that it might be a novel prognostic biomarker and therapeutic target. Two circRNAs (hsa_circRNA_103809 and hsa_circRNA_104700) have been implicated in CRC development and as potential biomarkers. RNA sequencing and in vitro experiments showed that circDDX17 was decreased in CRC tissues and might act as a tumor suppressor and novel biomarker. However, previous studies were limited to circRNAs that were differentially expressed between CRC and normal tissues, and their expression changes in adenomas were not clarified. Little is therefore known about circRNA-related ceRNAs predicting the adenoma to cancer transition.

We aimed to clarify the pathogenesis of CRC and identify potential biomarkers along the colorectal adenoma to CRC transition by comprehensive expression profiling of circRNAs, miRNAs, and mRNAs in colorectal adenoma, CRC, and adjacent normal tissues. We focused on mining common molecular events in colorectal adenoma and CRC, constructed circRNA-related ceRNA networks, and identified predictive biomarkers and potential therapeutic targets of CRC. Our findings may provide novel perspectives in CRC pathogenesis and suggest possible candidate biomarkers for the early diagnosis and treatment of CRC.

Results

Expression profiles of Candidate circRNAs, miRNAs, and mRNAs

A total of 575 differentially expressed circRNAs (DECs), 243 differentially expressed miRNAs (DEMs), and 3,950 differentially expressed mRNAs (DEGs) were identified in colorectal adenoma compared with adjacent normal tissues (adenoma versus normal), and 171 DECs, 53 DEMs, and 506 DEMs in CRC compared with normal tissues (cancer versus normal) (Figures 1A–1F). The overlapped DECs (84), DEMs (41), and DEGs (398) between adenoma versus normal and cancer versus normal were filtered out using Venn diagrams (see Figure S1). The log-transformed fold changes and p values of representative overlapped DECs, DEMs, and DEGs are listed in Table 1, and the details are described in Tables S1–S3, respectively.

Figure 1

Identification and Hierarchical Cluster Analysis of Candidate RNAs

Significantly changed RNAs were visualized in volcano plots. Red and green dots indicate upregulated and downregulated genes, respectively. RNA expression patterns were visualized by hierarchical cluster heatmap. Colors indicate expression values, with brighter red indicating higher values and brighter green lower values. (A and B) Differentially expressed circRNAs in adenoma compared with normal tissues (A) and in CRC compared with normal tissues (B). (C and D) Differentially expressed miRNAs in adenoma compared with normal tissues (C) and in CRC compared with normal tissues (D). (E and F) Differentially expressed mRNAs in adenoma compared with normal tissues (E) and in CRC compared with normal tissues (F). (G–I) Heatmaps of overlapped circRNAs (G), miRNAs (H), and mRNAs (I).

Table 1

Representative Overlapped RNAs between Adenoma and CRC

RNA Name	Adenoma versus Normal log2 (Fold Change)	Adenoma versus Normal p Value	Cancer versus Normal log2 (Fold Change)	Cancer versus Normal p Value
hsa_circ_0003915	3.93	0.00	3.46	0.00
hsa_circ_0008309	3.06	0.00	3.28	0.00
hsa_circ_0000660	3.34	0.00	2.68	0.02
hsa_circ_0000566	2.97	0.00	2.00	0.02
hsa_circ_0071411	−5.17	0.00	−2.13	0.04
hsa_circ_0049487	−6.75	0.00	−2.19	0.00
hsa_circ_0009130	−2.72	0.03	−2.30	0.02
hsa_circ_0001279	−5.46	0.00	−2.30	0.03
hsa_circ_0066875	−3.26	0.00	−2.30	0.01
hsa_circ_0007444	−1.90	0.00	−2.53	0.00
hsa-miR-449a	−4.96	0.02	−4.92	0.00
hsa-miR-187-3p	−4.53	0.00	−2.50	0.02
hsa-miR-126-5p	−4.27	0.00	−2.99	0.00
hsa-miR-126-3p	−3.49	0.00	−1.88	0.02
hsa-miR-449b-5p	−2.53	0.04	−2.30	0.03
hsa-miR-20a-5p	1.36	0.02	1.49	0.03
hsa-miR-17-3p	1.45	0.02	1.58	0.03
hsa-miR-17-5p	1.63	0.00	1.56	0.02
hsa-miR-135b-3p	3.01	0.00	3.09	0.01
hsa-miR-135b-5p	4.65	0.00	4.97	0.00
CDC6	1.42	0.00	1.37	0.05
CDK1	1.92	0.00	1.58	0.03
DDR2	−4.71	0.00	−1.69	0.03
EFS	−2.48	0.00	−1.15	0.03
FOXN3	−1.42	0.00	−1.23	0.00
IL6R	−1.83	0.00	−1.16	0.05
KLF9	−3.29	0.00	−1.81	0.05
MID2	−1.28	0.00	−1.17	0.04
SDK1	−1.39	0.00	−1.05	0.04
TCF21	−1.48	0.00	−1.16	0.04

We explored the expression patterns of the overlapped RNAs among the groups by hierarchical cluster analysis. Most of the overlapped RNAs exhibited similar expression patterns in adenoma and CRC tissues compared with adjacent normal tissues (Figures 1G–1I). For example, has-miR-135b-5p, has-miR-20a-5p, and has-miR-17-5p levels were elevated in both adenoma and CRC tissues, whereas hsa_circ_0066875, hsa_circ_0007444, hsa_circ_0049487, TNS1, and SDK1 levels were decreased in both adenoma and CRC samples.

Functional Enrichment and Protein-Protein Interaction (PPI) Analysis

GO and KEGG pathway enrichment analyses were performed to explore the main biological functions of the overlapped DEGs. Among the top 30 enriched GO terms and KEGG pathways (Figure 2; Tables S4 and S5), regulation of gene silencing, DNA damage checkpoint, cell cycle, p53 signaling pathway, etc. were related to tumorigenesis.

Figure 2

Functional Enrichment of Candidate DEGs

(A) Top 30 enriched GO terms; y axis represents GO terms, and x axis represents rich factor. Size and color of the bubble represent number of DEGs enriched in GO terms and enrichment significance, respectively. (B) Top 30 enriched KEGG pathways; y axis represents pathway names, and x axis represents rich factor. Size and color of the bubble represent number of DEGs enriched in the pathway and enrichment significance, respectively.

PPI analysis was performed to identify important hub genes among the overlapped DEGs, and a network with 301 nodes was generated (Figure 3A). A total of 180 genes (e.g., CASP8, TNS1, CXCL1, etc.) interacted with more than five other genes. The degree of each node was calculated using CytoNCA, and a sub-network with the top 50 genes (e.g., CDC6, CDK1, PTTG1, etc.) was presented in Figure 3B.

Figure 3

PPI Networks of Candidate DEGs

Nodes and edges represent genes and interactions, respectively. (A) PPI network with 301 nodes. (B) Sub-network with top 50 hub genes.

miRNA-mRNA Interaction and circRNA-Related ceRNA Network Analysis

To construct circRNA-related ceRNA networks, negatively regulated miRNA-mRNA pairs were first screened and visualized. A network with 24 downregulated miRNAs and 158 upregulated mRNAs was constructed (Figure 4A), e.g., decreased has-126-5p was correlated with two increased mRNAs (DTL and NDC1). Another network with 17 upregulated miRNAs and 86 downregulated mRNAs was also generated (Figure 4B), e.g., upregulated hsa-135b-5p interacted with eight downregulated mRNAs (TNS1, SDK1, and ALPK3, etc.). The interaction between miRNAs and circRNAs was predicted using starBase, and a ceRNA network (Figure 4C) was constructed including three downregulated circRNAs (hsa_circ_0049487, hsa_circ_0066875, and hsa_circ_0007444), three upregulated miRNAs (hsa-miR-135b-5p, hsa-miR-20a-5p, and hsa-miR-17-5p), and 25 downregulated mRNAs (ALPK3, TNS1, SDK1, etc.).

Figure 4

miRNA-mRNA Regulatory Networks and circRNA-Related ceRNA Network

Rectangles represent circRNAs, triangles represent miRNAs, ellipses represent mRNAs. Red nodes represent upregulated and green nodes represent downregulated. (A) Downregulated miRNA with upregulated mRNA regulatory network. (B) Upregulated miRNA with downregulated mRNA regulatory network. (C) circRNA-related ceRNA network.

The Cancer Genome Atlas (TCGA) Database and Real-Time PCR Validation

To validate our findings, the RNA expression levels in the ceRNA network were compared with the TCGA database and confirmed by real-time PCR. The expression patterns of six downregulated mRNAs and two upregulated miRNAs were consistent with the TCGA data (Figures 5A–5H). Another 10 differentially expressed RNAs (three downregulated circRNAs, one upregulated miRNA, and six downregulated mRNAs) were validated by real-time PCR in both colorectal adenoma and cancer tissues (Figures 5I and 5J). The RNA sequencing and real-time PCR results were in good accordance.

Figure 5

Representative RNAs Expression Validation

(A–H) RNA expression levels in cancer and normal samples according to TCGA database, (A) DDR2, (B) EFS, (C) FOXN3, (D) IL6R, (E) KLP9, (F) TCF21, (G) hsa-miR-17-5p, (H) hsa-miR-20a-5p. (I and J) Real-time PCR validation. The x axis represents RNA names, and the y axis represents log2 (fold change) based on the ratio of adenoma and normal average expression values (I) or the ratio of CRC and normal average expression values (J). Brown bars represent real-time PCR data, and blue points represent RNA sequencing data.

Receiver Operating Characteristic (ROC) Curve Analysis

The expression levels of three circRNAs (hsa_circ_0049487, hsa_circ_0066875, and hsa_circ_0007444) in the ceRNA network were confirmed by real-time PCR. Literature retrieval revealed that circRNA-related miRNAs and mRNAs in the ceRNA network might play a role in cancer progression, suggesting potential roles for the above three circRNAs in cancer. Therefore, ROC curve analysis was employed to evaluate the diagnostic values of the candidate circRNAs in colorectal adenoma and CRC. The area under the curve (AUC) of has_circ_0066875 was 0.806, 0.907 in cancer and adenoma, respectively (Figure 6A). The AUC of hsa_circ_0007444 was 0.771, 0.682 in cancer and adenoma, respectively (Figure 6B). The AUC of hsa_circ_0049487 was 0.809, 1 in cancer and adenoma, respectively (Figure 6C). These results suggested that the above three circRNAs possessed good diagnostic values in colorectal adenoma and cancer and might thus serve as early biomarkers to predict CRC.

Figure 6

ROC Curve Analysis of Three circRNAs

Blue dotted lines represent diagnostic values for discriminating cancer from normal tissues; green dotted lines represent diagnostic values for discriminating adenoma from normal tissues. (A–C) hsa_circ_0066875 (A), hsa_circ_0007444 (B), and hsa_circ_0049487 (C).

Discussion

The pathogenesis of CRC is complicated, and an adenoma-carcinoma sequence hypothesis has been proposed. The duration of the adenoma-to-carcinoma transition offers an opportunity to reduce the incidence and mortality of CRC if specific biomarkers can be employed to predict this transition. Identification of specific molecular signatures is urgently needed, and little is currently known about circRNA-related signatures.

In the present study, we employed comprehensive RNA sequencing to obtain circRNA, miRNA, and mRNA expression profiles for adjacent normal, colorectal adenoma, and CRC tissues. We identified several candidates by focusing on the RNAs that are commonly changed in colorectal adenoma and CRC. Integrated analysis of functional enrichment and PPI obtained a batch of hub genes (e.g., CDC6, CDK1, PTTG1, etc.) related to cancer progression. Cell division cycle 6 (CDC6) is an important regulator of cell cycle, and its aberrant expression in human cells poses a serious risk of carcinogenesis. CDC6 has been reported as a potential prognostic marker of lung adenocarcinoma and prostate cancer., In the present study, CDC6 expression levels were elevated in both colorectal adenoma and CRC tissues, implying a potential role in the adenoma-to-carcinoma transition. Cyclin-dependent kinase1 (CDK1) is a member of the CDK family associated with cancer development. Previous bioinformatics data mining identified CDK1 as an important hub gene with elevated expression in all CRC stages compared with adjacent normal cells. Our findings were partly consistent with the above previous studies. We noted that the hub gene CDK1 was upregulated in both human colorectal adenoma and CRC compared with adjacent normal tissues. Pituitary tumor transforming gene-1 (PTTG1) was also overexpressed in colorectcal carcinoma tissues and may serve as an oncogene in CRC. In the current study, PTTG1 expression was elevated in both colorectal adenoma and CRC tissues. Further investigation of the hub genes in the present study may obtain novel predictive biomarkers and therapeutic targets to treat CRC.

Long noncoding RNA (lncRNA)-related ceRNA networks have been extensively investigated in human cancers such as hepatocellular carcinoma, breast cancer, as well as CRC.22, 23, 24 However, the involvement of circRNA-related ceRNA networks in CRC, especially their effects in colorectal adenoma and its transition to CRC, remains unclear. In the present study, we constructed three circRNA-related ceRNA networks using candidate RNAs that overlapped between adenoma and CRC. We noted that two circRNAs (hsa_circ_0007444 and hsa_circ_0066875) might be ceRNAs of eight mRNAs (ADCY5, ALPK3, SDK1, MID2, TNS1, RBPMS, GNAZ, and RRAS) by binding to hsa-miR-135b-5p. Although little is known about the biological functions of the above two circRNAs, their related miRNAs and mRNAs have revealed potential roles in cancer. For example, hsa-miR-135b-5p might modulate the APC gene and influence gastric carcinogenesis, while ADCY5 has been associated with several human carcinomas, including breast cancer, prostate cancer, and CRC., Regarding the ALPK3 gene, an analysis of Gene Expression Omnibus (GEO) datasets indicated its relationship with the metastasis of osteosarcoma, which was the most common primary solid bone malignancy. The role of TNS1 in different types of cancer is controversial. TNS1 was downregulated in human kidney cancer and breast cancer but upregulated in CRC. We noted that TNS1 was downregulated in both colorectal adenoma and CRC tissues in the current study. However, further large-scale investigations are needed to confirm the role of TNS1 in CRC progression. The RBPMS gene has been implicated in the migration, invasion, and apoptosis of human CRC cells, and the remining four genes, MID2, SDK1, GNAZ, and RRAS, have also been associated with various human cancers.33, 34, 35, 36

In addition, we observed that hsa_circ_00049487 might be a ceRNA of 18 mRNAs (ALPK3, FOXN3, TCF21, etc.) by sponging two miRNAs (hsa-miR-17-5p and hsa-miR-20a-5p). The biological functions of hsa_circ_00049487 have not yet been determined but may be inferred from the related miRNAs and mRNAs. A large-scale pan-cancer analysis identified hsa-miR-17-5p as an oncogenic miRNA, and another study reported that hsa-miR-17-5p and hsa-miR-20a-5p were elevated in human CRC tissues, which might serve as important discriminators for CRC. In the present study, the expression levels of hsa-miR-17-5p and hsa-miR-20a-5p were upregulated in both colorectal adenoma and CRC tissues, which were partly consistent with these previous data. Of interest, the binding targets of the above two miRNAs have been demonstrated in cancer progression. For instance, FOXN3 was downregulated in colon cancer tissues, and its loss of function promoted colon cancer cell growth, invasion, and metastasis, indicating its suppressive role in colon cancer. TCF21 is also a tumor suppressor gene in both human CRC tissues and cells, and its downregulation indicated a poor prognosis in CRC patients., We consistently found that FOXN3 and TCF21 expression levels were reduced in both colorectal adenoma and CRC tissues.

Overall, our findings suggested that hsa_circ_0007444, hsa_circ_0066875, and hsa_circ_00049487 might play roles in both colorectal adenoma and CRC. Finally, ROC curve analysis revealed that the above three circRNAs could be used to discriminate colorectal adenoma/CRC from adjacent normal tissues and might thus act as early predictive biomarkers to prevent CRC. Although novel circRNAs have been strongly associated with the transition from colorectal adenoma to cancer, further studies would be worthy to be investigated. First, our data were based on RNA sequencing expression profiles, and the underlying functions and mechanisms of the circRNA-related ceRNAs are to be clarified. Second, the numbers of clinical samples were limited in the present study, and larger-scale validation is needed to verify these findings.

In conclusion, we comprehensively analyzed the expression profiles of circRNAs, miRNAs, and mRNAs from colorectal adenoma, CRC, and adjacent normal tissues. We mined RNAs that are commonly changed in colorectal adenoma and CRC, constructed circRNA-related ceRNA networks, and identified three novel circRNAs as potential predictive biomarkers and possible therapeutic targets of CRC. Our findings may offer new perspectives for understanding the pathogenesis of the transition from colorectal adenoma to CRC and provide novel biomarkers and candidate targets for the early diagnosis and treatment of CRC.

Materials and Methods

Participants and Samples

A total of 90 samples were collected (30 colorectal adenoma, 30 CRC tissues, and 30 adjacent normal tissues) from patients who underwent surgical treatment in Longhua Hospital, Shanghai, China. All resected specimens were snap frozen in liquid nitrogen and stored at −− 80 °C. This study was approved by the Ethics Committee of Longhua Hospital, and informed consent was obtained from all participants.

RNA Extraction and RNA Sequencing Data Acquisition and Processing

Total RNA was extracted from five colorectal adenoma, five CRC, and five adjacent normal tissues using TRIzol reagent (Life Technologies, CA, USA) following the manufacturer’s instructions. The RNA quality was assessed by spectrophotometry and denaturing agarose gel electrophoresis. The RNA sequencing data were acquired as described previously. In brief, the sequencing libraries for circRNAs, miRNAs, and mRNAs were prepared, purified, qualified and sequenced using Illumina Hiseq 4000. Clean reads were obtained and aligned to the human genome. Significantly differentially expressed RNAs were identified using DEseq2 software with a threshold of 2-fold change and p value less than 0.05.

DECs, DEMs, and DEGs were filtered out in pairwise groups (colorectal adenoma versus adjacent normal and CRC versus adjacent normal). The differentially expressed RNAs were visualized on volcano plots constructed by plotting-log10 (p value) on the y axis and log2 (fold change) on the x axis. Overlapped differentially expressed RNAs between colorectal adenoma and CRC tissues (compared with adjacent normal tissues) were obtained by Venn diagrams and subjected to further analysis. Hierarchical cluster analysis was performed to exhibit the expression patterns in pairwise groups.

Functional Enrichment and PPI Analysis

We further explored the main biological functions of the identified DEGs by functional enrichment analysis based on the GO and KEGG databases., PPI networks were screened based on the String database and visualized using Cytoscape software (version 7.0). The Cytoscape plug-in app CytoNCA was used to calculate the node degree, and the top 50 representative hub nodes with higher degree were filtered out.

miRNA-mRNA Interaction and circRNA-Related ceRNA Network Construction

miRNA-mRNA interaction networks were constructed as described previously. To construct the circRNA-related ceRNA networks, interactions between circRNAs and miRNAs were predicted using starBase. According to the ceRNA hypothesis, circRNAs act as miRNA sponges and negatively regulate miRNA-mediated gene silencing. Negatively interacting circRNA-miRNA pairs were filtered out, and circRNA-related ceRNA networks were constructed and visualized using Cytoscape software (version 7.0).

RNA Expression Validation by TCGA Database and Real-Time PCR

We validated our findings by examining the expression levels of eight RNAs in ceRNA networks by large-scale TCGA data (471 CRC and 41 normal tissues) using the online tool starBase and validated a further 10 RNAs by real-time PCR as described previously. The primer sequences are listed in Table 2.

Table 2

Primer Sequences in Real-Time PCR Experiments

RNAs	Primer (5′ to 3′)
hsa_circ_0049487	forward, CATACACAGGTGCAGTCC
	reverse, GGCTTCACCCCATACTTG
hsa_circ_0066875	forward, AGGAGCTGTCACGGGAAGT
	reverse, GAATGAAGCCTCGTGTGG
hsa_circ_0007444	forward, AAGTTGAAAGATTCTGGGGATG
	reverse, TGTGACGCTTCAGCCTTT
hsa-miR-17-5p	forward, CAAAGTGCTTACAGTGCAGGTAG
	reverse, ATCCAGTGCAGGGTCCGAGG
hsa-miR-20a-5p	forward, TAAAGTGCTTATAGTGCAGGTAG
	reverse, ATCCAGTGCAGGGTCCGAGG
hsa-miR-135b-5p	forward, TATGGCTTTTCATTCCTATGTGA
	reverse, ATCCAGTGCAGGGTCCGAGG
U6	forward, AGAGAAGATTAGCATGGCCCCTG
	reverse, ATCCAGTGCAGGGTCCGAGG
MID2	forward, CCATGACCACCTCCCAAT
	reverse, AAACAAACAGCCAACCCT
TNS1	forward, TGGAGAAGTCGGGCAGAG
	reverse, AGAAGCGAAGGATGTCAGTG
RBPMS	forward, CGGAGTCCAGGGTAAATTAGTAGCA
	reverse, CCAGTTGTGAATAAGCCATAGGTAGAAG
ALPK3	forward, TGAGGCAGAAGTCGGTGG
	reverse, GATTCGGGTGGAGCAGTT
ADCY5	forward, ATCGCCCAGGCTGTAGTT
	reverse, AAGATCGTGCGGTCCAAG
SDK1	forward, ACGATCAGGCTGGACAGAG
	reverse, GTGAATGAGGCGGGCTAC
ACTB	forward, GAAGAGCTACGAGCTGCCTGA
	reverse, CAGACAGCACTGTGTTGGCG

Statistical Analysis

Statistical analysis was performed using MedCalc statistical software (version 15.8). p values less than 0.05 was considered statistically significant. ROC curves were generated, and the AUC was calculated to assess the diagnostic values of circRNA expression for discriminating colorectal adenoma/CRC from normal tissues.

Author Contributions

G.J. and W.Z. conceived, designed, and supervised the study. M.Z., Y.D., Z.Y., Y.L., L.Z., and Y.X. collected samples. M.Z. and Y.D. performed the experiments and analyzed the data. M.Z., W.Z., and G.J. wrote the paper. All authors reviewed and approved the final manuscript.

Conflicts of Interest

The authors declare no competing interests.