MiR-4319 Suppress the Malignancy of Triple-Negative Breast Cancer by Regulating Self-Renewal and Tumorigenesis of Stem Cells

Triple-negative breast cancer (TNBC) constitutes 10% to 20% of all breast cancers [1, 2] and is characterized by tumors with a lack of estrogen receptor (ER), progesterone receptor (PR), human epidermal growth factor receptor 2 (HER2), and higher rates of metastasis and poorer prognosis [3]. The incidence rate of TNBC has increased continuously over the past decades [4], especially in young women. A persistent problem among the young cohort is that some women ignore their regular physical examination, which makes early mammographic screening more difficult [5, 6]. Women in whom TNBC is diagnosed rarely benefit from endocrine therapy or chemotherapy, and as a result, TNBC has poor prognosis and more than 70% mortality in metastatic TNBC [1]. Furthermore, TNBC appeared with higher tumor heterogeneity than non-TNBC. Recently, TNBC has been divided into six subtypes based on the Lehmann/Pietenpol classification [7], and current chemotherapy treatments already evoke diverse responses and clinical outcomes in patients with TNBC. Consequently, there is still no optimized standard form of chemotherapy for TNBC, leading to an urgent need for the development of effective treatments for this malignancy. Several studies have suggested that cancer stem cells (CSCs) correlated positively with the malignant degree of breast cancer, and might be responsible for metastasis and drug resistance of patients with TNBC [8, 9].

CSCs, a subpopulation of cancer cells, are capable of initiating tumorigenesis by undergoing self-renewal and differentiation [10]. Breast cancer stem cells (BCSCs) have been defined and prospectively isolated by the cell surface markers CD44+, CD24-, and ALDH+ [11]. BCSCs play important roles in breast cancer, such as promoting tumor initiation, stimulating metastasis, and inducing chemoresistance [12]. In addition, promising evidence shows the existence of a positive association between the proportion of stem cells in tumors and disease progression, and that high levels of CSCs correlate with risk of metastasis and poor clinical outcome [13]. TNBC has been identified to comprise more CD44+/CD24-cells (CSCs) than non-TNBC in solid tumors [14, 15], which might contribute to chemoresistance and metastasis in patients with highly malignant TNBC.

It has been confirmed that microRNAs play a role in regulating biological processes through downstream target genes. Specifically, anti-microRNAs agents might be a new strategy for the development of clinical antitumor therapy [16, 17]. Over the past few decades, it has been demonstrated that appropriate combinations of microRNAs were closely related to the malignant progression of breast cancer [18]. Several microRNA clusters, involving miR-374-5p, miR-126-3p, and miR-27b-3p, are linked to development, recurrence, and response to therapeutic treatments in breast cancer [19, 20]. Previous studies have also suggested that aberrant microRNA expression contributes to the dysfunction of self-renewal and proliferation in cancer stem cells [21, 22], such as let-7 and miR-30, which could reduce the ability of self-renewal and mammosphere formation in BCSCs [22, 23]. Over the past few decades, protein markers such as ER, PR, and HER2, which are well known to stimulate the development of non-TNBC, have been used to differentiate TNBC and non-TNBC, and therapeutic agents, which targeted ER or HER2, inhibited non-TNBC efficiently. There is little progression in the development of treatment for TNBC. It has been proposed that CSCs might contribute to the malignancy of TNBC, and the correlation between CSC proportion and microRNAs might be helpful to find the potential therapeutic targets for TNBC by controlling CSC function.

Currently, there is little information about the TNBC-specific microRNAs in regulating CSC ability. We detected the TNBC-derived, CSCs-related microRNAs expression panel, and found an unknown miR-4319 especially dysregulated in TNBC CSCs. In this study, we aimed to identify the effect of miR-4319 on BCSCs and determine the undergoing mechanism.

The human breast cancers BT474, MCF-7, SKBR3, MDA-MB-468, MDA-MB-231 and MDA-MB-453 were all obtained from American Type Culture Collection (ATCC, Manassas, VA). BT474 and SKBR3 were cultured in RPMI-1640 medium (Gibco); MDA-MB-453 was cultured in L15 medium (Gibco); MCF-7, MDA-MB-231 and MDA-MB-468 were cultured in DMEM medium (Gibco); supplemented with 10% (v/v) fetal bovine serum (Gibco), 100 U/ml penicillin (Gibco), and 100 U/ml streptomycin (Gibco) at 37°C in a humidified atmosphere containing 5% CO₂ [24].

Anti-Human/Mouse CD44 FITC (eBioscience) and Anti-Human CD24 PE (eBioscience) were used for flow cytometric analysis [25]. Cells were harvested by trypsinization and washed once in PBS. Next, cells were counted and resuspended in cold PBS at a concentration of 1.0 * 10⁷ cells per 100ul. Staining was performed by incubating 100ul cells with 20ul antibody in 4°C for 25 min.

The cancer stem cells were plated into ultra-low attachment six well plates (corning 3471) at a density of 2000 cells per well [26]. The cancer stem cells were cultured with MammoCult^TM Human Medium Kit (Stem cell technologies, Vancouver, BC, Canada) according to manufacturer’s guidelines [27]. Mammospheres were resuspended in trypsin/EDTA and passaged by centrifuging at 300g for 10 min in 7-14 days. Secondary spheres were used to analyze the efficiency of mammosphere-formation by assessing the number of spheres formed per cells seeded [28]. Then formed spheres were counted manually and representative images were acquired [26] using Leica DMI3000B microscope with Leica Microsystems Framework.

Total RNA was isolated from non-TNBC stem cells (BT474, MCF-7, SKBR3) and TNBC stem cells (MDA-MB-468, MDA-MB-231, MDA-MB-453) using Trizol (Ambion, Carlsbad, CA), and then assessed by Agilent 2100 BioAnalyzer (Agilent Technologies, Palo Alto, CA). Two small RNA libraries of BT474 CSC, and MDA-MB-468 CSC were established by following the described methods [29]. We filtered the low quality reads according to base quality value. Small RNAs (18–30 nt), acquired from the total RNA, with 5’ and 3’ adaptors ligated to them, were subsequently transcribed into cDNA using RT Master Mix (Takara, Dalian, China) followed by PCR amplification. The PCR products were purified and submitted to Solexa sequencing (Illumina, CA) Beijing Genomics Institute (BGI, Shenzhen, China). Reads (counts < 2) were discarded when calculating the microRNA expression. Fold-change (≥2.0 or ≤0.5) of the reads of sequenced microRNAs was employed to study the differential expression profiles of the CSC microRNAs in BT474 and MDA-MB-468 groups based on TagCount_ALL_Isoform data.

The initial sequence reads were generated by Illumina Genome Analyzer. To identify conserved microRNAs, the filtered sequences were originally applied to miRBase 21.0 (http://www.mirbase.org/) with BLASTN. The criterion was implemented based on a reported microRNA protocol [30]. All unannotated RNAs were subjected to novel microRNA prediction using Mireap software (https://sourceforge.net/projects/mireap/). The expression of known microRNAs was compared between the two groups to screen differentially expressed microRNAs.

The possible targets of differentially expressed microRNAs were predicted via Mirbase, Target-Scan, and miRanda programs with the default parameters. MicroRNA target intersections of the three programs were regarded as potential microRNA target genes. The MAS3 program (http://bioinfo.capitalbio.com/mas3/) was used to annotate the functions of potential microRNA target genes via the Gene Ontology (GO) (http://www.geneontology.org/) project and the relative pathways were analyzed via the Kyoto Encyclopedia of Genes and Genomes (KEGG) program (http://www.genome.jp/kegg/pathway.html).

A stem-loop primer based PCR assay was used to validate the differentially expressed microRNAs. Total RNA was extracted from TNBC CSCs (MDA-MB-468, MDA-MB-231, MDA-MB-453) and non-TNBC (BT474, MCF7, SKBR3) using Trizol (Life Technologies, Shanghai, China), then was subjected to reverse transcription and qRT-PCR using the Bulge-LoopTM microRNA qRT-PCR Starter Kit (Ribobio, Guangzhou, China). qRT-PCR was performed on ABI 7900 PCR System (Applied Biosystems, USA). After denaturation, the reactions were conducted as the following: 2s at 95°C, 20s at 60°C and 10s at 70°C for 40 cycles. Results were expressed as ΔCt values and the expression of snRNA U6 was served as normalization control [31]. The RNA expression of the specific stem cell properties (CD24, CD44, Oct4, ALDH, and Nanog [32]) was analyzed by qRT-PCR (SYBR Green Master, Takara, Dalian, China). qRT-PCR was performed with a Roche LightCycler 96 Real-Time PCR System. Data was analyzed by the 2^-ΔΔCT method with normalization to GAPDH. Partial primer sequences are listed at Table 1.

The CD44+CD24- human BCSCs were isolated by flow cytometry and then were infected by 20 moi of miR-4319-expressing lentivirus, miR-4319-inhibiting lentivirus or control lentivirus by spin infection for 2 hours followed by incubation at 37°C to generate stable expression cell lines (MDA231-miR-4319-OE, MDA231-miR-4319-KD and MDA231-NC). The expression of miR-4319 in MDA231-miR-4319-OE, MDA231-miR-4319-KD and MDA231-NC was confirmed via real-time PCR.

Infected CD24-CD44+ cells (MDA231-miR-4319-OE, MDA231-miR-4319-KD and MDA231-NC) were injected into mammary fat pad of female immunodeficient mice mixed with Matrigel (BD Biosciences) at the concentration of 5*10⁴/150ul. After 45 days, all the mice were sacrificed to harvest the xenograft. All experiments were carried out under the approval of the Administrative Panel on Laboratory Animal Care of Stanford University. The liver and lung tissues were fixed in 4% paraformaldehyde in continuing dehydration in ascending grades of Ethanol (Merck, Germany), cleared and embedded in xylol and paraffin, respectively. Sections of 5 µm were taken, stained with Hematoxylin-Eosin (H-E), and investigated under the light microscope (Leica DMI3000B microscope).

MDA231 cells were co-transfected with lentivirus (miR-4319-expressing, negative control) and 0.5µg of dual luciferase reporter plasmid containing wild-type or mutated E2F2 3’UTRs (Synbio technologies, Suzhou) using Lipofectamine 2000 (Life Technologies Corporation). 48h post-transfection luciferase activity was quantified following the instruction of Dual Luciferase Reporter kit (Promega Inc).

MicroRNAs expressions in the treated and the control cells were identified by applying one way ANOVA (analysis of variance) method followed by post-hoc Tukey’s analysis. A p-value of less than 0.05 was considered statistically significant (p value< 0.05) and each experiment was conducted in triplicate. Computer-based calculations were performed using SPSS version 20.0 (SPSS Inc, Chicago, Illinois, USA).

BCSC-enriched subpopulations were sorted from MDA-MB-468 (ER-, PR-, HER2-), MDA-MB-231 (ER-, PR-, HER2-), MDA-MB-453 (ER-, PR-, HER2-), BT474 (ER+), MCF-7 (ER+), and SKBR3 (HER2+) with CD24 and CD44 cell-surface markers by flow cytometric analysis (Fig. 1A, data showed partially). The proportion of CD44+CD24- cells in TNBC was higher than non-TNBC, which is consistent with previous reported results [33]. The ability of tumorsphere formation was also examined and the data showed that the isolated single cell (both TNBC- and non-TNBC-CSCs) could form dominant tumorspheres in 7 to 14 days (Fig. 1B). However, TNBC-derived CSCs (MDA-MB-468, MDA-MB-231, MDA-MB-453) appeared to have a stronger ability to form tumorsphere in less than 7 days. The number of tumorspheres was counted after 7 days. The results showed that the tumorsphere formation ability of MDA-MB-468 (54.4%±10.96%) was stronger than that of BT474, and the tumorsphere formation ratio of MDA-MB-468 (54.4%±10.96%) was higher than BT474 (36.3%±7.02%) (p > 0.05). The same situation was also observed in MDA-MB-231, MDA-MB-453, MCF-7, and SKBR3. The expression of the specific surface markers (CD44, ALDH, Oct4, and Nanog) was significantly increased, whereas the expression of CD24 was significantly decreased in cancer stem cells (Fig. 1B, partial data shown).

Total RNA of each sample was used to prepare the microRNA sequencing library. The Solexa [34] data showed that microRNAs were the major components of small RNAs (< 30 bp) in the BCSCs. The 193 microRNAs were expressed differentially between the two small RNA libraries (P< 0.05), in which 112 microRNAs were upregulated and 81 were downregulated in MDA-MB-468 CSC comparing to BT474 CSC (Fig. 2A, 2B).

Stem-loop real-time PCR quantification was implemented to verify 10 microRNAs that were selected from the Solexa sequencing microRNA expression results in TNBC (MDAMB-468, MDA-MB-231, MDA-MB-453) and non-TNBC (BT474, MCF7, SKBR3) CSCs. Four of 10 (hsa-miR-205, hsa-miR-181-2-3p, hsa-miR-222-3p, and hsa-miR-135) were upregulated and 2 microRNAs (hsa-miR-30a and hsa-miR-4319) were stably downregulated in TNBC CSCs (MDA-MB-468, MDA-MB-231, MDA-MB-453) (shown in Fig. 3).

To further understand the biological functions of the microRNAs, target genes were applied to MiRbase, MiRanda, and TargetScan (http://www.targetscan.org); there were 236 target genes down-regulated and 523 target genes upregulated with a value of P< 0.05. The microRNA-target genes were associated with the following biological processes: stem cell development, embryonic morphogenesis, single-organism development, multicellular organismal development, cell differentiation, system development, organ development, regulation of metabolic process, regulation of stem cell differentiation, and response to chemicals (Fig. 4A). The microRNA-target genes were also involved in cellular components (Fig. 4B) and potentially regulated almost 201 gene ontology (GO) terms of molecular function of which the top three were mRNA binding, protein binding, and enzyme binding (Fig. 4C). The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis showed that a total of 55 pathways were identified in which most were associated with the function of CSCs, such as signaling pathways regulating pluripotency of stem cells, Notch signaling pathway, mitogen-activated protein kinase (MAPK) signaling pathway, estrogen signaling pathway, Janus kinase/signal transducers and activators of transcription signaling pathway, forkhead transcription factors of the O class signaling pathway, mechanistic target of rapamycin signaling pathway, and transforming growth factor beta signaling pathway (Fig. 4D). The brief mechanism diagram was listed in Fig. 4F, and the orange- labeled proteins were predicted to be dysregulated between TNBC-CSC and non-TNBC-CSC, which indicated that these proteins probably participated in maintaining the high pluripotency of TNBC-CSC. Most importantly, our analysis suggested that the microRNAs played important roles in tumor initiation and development.

We performed function analysis of miR-4319 in MDA-MB-231 CSCs. The expression of miR-4319 in stable cell lines (MDA231-miR-4319-OE, MDA231-miR-4319-KD, and MDA231-NC) was confirmed by quantitative real-time polymerase chain reaction (qRT-PCR) (Fig. 5B). Compared to the miR-4319-NC group, a decrease in tumorsphere-forming efficiency was observed in the miR-4319-OE group and an increase was observed in the miR-4319-KD group. Tumorspheres in the miR-4319-KD group grew faster than that in the miR-4319-NC group and the results in the miR-4319-OE group showed the reverse. The number of tumorspheres in suspension was expanded in the miR-4319-KD group compared to the control group (Fig. 5A, 5C). qRT-PCR demonstrated that inhibition of miR-4319 markedly increased the transcript abundance of Sox2, Oct4, Nanog, and CD133, which were specific for isolated CSCs (Fig. 5D).

The xenograft model of TNBC CSCs was adopted to assess the role of miR-4319 in tumor initiation as xenotransplantation remains the gold standard for the balance between BCSC self-renewal and differentiation in vivo. A higher number of macroscopic hepatic metastases could be observed in the MDA231-miR-4319-KD group, and few macroscopic hepatic metastases could be seen in the MDA231-miR-4319-OE group (Fig. 6A). Compared to control, the overexpression of miR-4319 inhibited tumor formation, and knockdown of miR-4319 promoted tumor initiation and extensive metastases, which indicated that miR-4319 was a regulator of CSC functions. MiR-4319 knockdown resulted in a higher number of growing tumors than in the control group, which showed the strong ability of miR-4319 on regulation of the tumorigenicity expected from TNBC CSCs (Fig. 6B). Hematoxylin and eosin staining showed metastatic focus in liver (Fig. 6C) and lung (Fig. 6D); the same results were seen via visual inspection. These results indicate that miR-4319 knockdown could significantly enhance the tumorigenesis of stem cells and consequently be responsible for extensive metastasis in TNBC.

Gene expression profiles of MDA231-miR-4319-OE group and MDA231-NC group were determined using Illumina Hiseq 4000. Among the various genes expressed, E2F2 has been known to maintain the pluripotency of human embryonic stem cells and three miR-4319 binding sites (M1, M2, M3) on E2F2 3′ untranslated region (UTR) were predicted by TargetScan (Fig. 7A). Both qRT-PCR and western blot results showed that overexpression of miR-4319 in MDA231-CSCs significantly decreased the expression of E2F2 (Fig. 7B, 7C). The 3′UTR of E2F2 was cloned into pGL3-control vector and the specific binding sites were confirmed by luciferase assays via mutation of the three sites. Overexpression of miR-4319 greatly reduced luciferase activity of the reporter gene by almost 50% and when the mismatch mutation was introduced into the region of M1, M2, and M3 at the same time, the inhibition of luciferase activity was almost abolished (Fig. 7C). The results demonstrated that miR-4319 could bind to the three sites (M1, M2, and M3) on the E2F2 3′UTR to regulate its biological function in CSCs.

A rescue experiment has been performed to confirm the function of E2F2 in the MDA231-miR-4319-OE/KD group. MDA231-miR-4319-OE was transfected with E2F2 vector or empty vector (E2F2+ or control-1); MDA231-miR-4319-KD was transfected with E2F2-small interfering RNA or scrambled control (E2F2 siRNA or control-2). The expression levels of E2F2 were determined by qRT-PCR and western blot analysis (Fig. 8A). E2F2 overexpression accelerated the formation of tumorspheres and increased the expression of the specific stem cell surface markers (Sox2, Oct4, Nanog and CD133), and both functions were suppressed by miR-4319-OE in TNBC CSCs; in contrast, E2F2 depletion reduced the CSC functions in miR-4319-KD CSCs (Fig. 8B). MiR-4319 inhibited the formation of tumorspheres in TNBC CSCs, which was partially overcome by overexpressing E2F2.

These results showed that miR-4319 could inhibit the stem cell characteristics of TNBC CSCs through regulating E2F2 expression levels.

TNBC is a highly aggressive cancer accompanied by poor prognosis, with no effective treatments currently available. Additionally, women with TNBC are more likely to experience chemotherapy drug resistance. Previous studies have suggested that CSCs were positively related to the malignant degree of breast cancer and might be responsible for metastasis and drug resistance in TNBC [8, 9]. In this study, we isolated CD44+CD24- cells to enrich CSCs from TNBC (MDA-MB-468, MDA-MB-231, MDA-MB-453) and non-TNBC (BT474, MCF-7, SKBR3) by flow cytometric analysis. The results of microsphere formation showed that TNBC-derived CSCs showed a stronger ability to form microspheres and expression of specific stem cell markers was higher in TNBC-derived CSCs than non-TNBC derived CSCs (Fig. 1), which is consistent with previous studies [33].

Emerging evidence showed that microRNAs could regulate key features of BCSCs as oncogenes or tumor suppressor genes. Dysregulated expression levels of microRNAs have been related to poor prognosis and poor therapeutic efficiency in patients with breast cancer [35]. However, there is little information about the TNBC-specific microRNAs in regulating CSC ability. In this study, we detected the TNBC-derived, CSCs-related microRNAs expression panel utilizing Solexa sequencing, and the data showed that 112 microRNAs were upregulated and 81 were downregulated in MDA-MB-468–derived CSCs compared to BT474 CSCs (Fig. 2), including 13 novel microRNAs and 28 microRNAs that have been identified in regulating TNBC malignancies, which might play important roles in the control of CSC hallmark functions, such as self-renewal, tumorsphere formation, and metastasis. We performed bioinformatic analysis for further study and the results of GO and KEGG pathway analysis showed that these microRNAs were mainly involved in pluripotency maintenance of CSCs. Biological process of the microRNA-target genes most contained stem cell development, embryonic morphogenesis, and regulation of stem cell differentiation. The Notch pathway and MAPK signaling pathway have been reported to regulate the self-renewal of stem cells and dysregulation [36]. The results of GO and KEGG pathway analysis indicated that target genes of these differently expressed microRNAs probably were involved in maintenance of stemness and contributed to the differential properties (such as self-renewal, microsphere formation, chemoresistance, and invasion) between TNBC-CSCs and non-TNBC-CSCs (Fig. 4A-4F).

The expression levels of microRNAs (cutoff level of fold-change ≥2.0 or ≤0.5), were verified by qRT-PCR and 6 of 10 showed stable expression in TNBC CSCs (Fig. 3). MiR-4319 is an unknown miRNA, and showed a stable low expression level in TNBC CSCs, compared to non-TNBC CSCs. It has been reported that miR-4319 was related to acute myeloid leukemia [37], whereas the brief function and mechanism are not clear. In this study, we investigated the effect of miR-4319 in TNBC CSCs as well as the underlying mechanism. Overexpression of miR-4319 could significantly suppress tumorsphere formation and growth in suspension (Fig. 5A, 5C). The specific stem cell markers also were downregulated simultaneously (Fig. 5D). The results showed the reverse in the miR-4319-KD group. Nude mouse transplantation tumor experimentation showed the same results as in vitro experiments. Compared to control, overexpression of miR-4319 inhibited the capacity of self-renewal, tumor formation, and extensive metastases of TNBC CSCs whereas miR-4319 knockdown could significantly enhance the tumorigenicitys of stem cells and consequently be responsible for extensive metastasis in TNBC (Fig. 6). Together, these findings supported the hypothesis that miR-4319 was involved in regulating biological function of stem cells and inhibited the tumorigenesis, progression, and metastasis of TNBC CSCs.

MicroRNAs participate in a series of cellular processes by binding to 3′ UTR of target mRNA, and then inducing targeted messenger RNA (mRNA) slicing by Dicer complex in the cytoplasm [38]. Bioinformatics analysis results (TargetScan) and microarray data indicated that E2F2 was the potential target gene of miR-4319. Furthermore, there were three miR-4319 binding sites (1446,bp ∼1452,bp, 1481,bp∼1488,bp, and 1863,bp∼1870,bp) in 3′ UTR of E2F2 mRNA, and the luciferase assay results showed that miR-4319 could bind to each site to reduce the expression of E2F2. Moreover, overexpression of E2F2 abolished the effect of miR-4319 on the self-renewal and tumorigenicity in TNBC CSCs. E2F2 is a member of the E2F family of transcription factors. It plays a central role in regulating G1/S transition and progression through S phase, promoting cellular transformation [39]. Previous studies have reported that E2F2 is a transcription factor with important functions in the control of embryonic development and differentiation [40]. E2F2 is also essential for self-renewal of stem cells by inducing the transcription of oncogenes typically associated with poor prognosis, such as BMI1 and MELK [41, 42]. Hence, miR-4319 might reduce the stemness of TNBC CSCs mainly through inhibition of E2F2 expression.

Our data elucidated the differences in miRNA expression patterns between TNBC and non-TNBC CSC subpopulation, and further investigated the effect of miR-4319 on BCSCs as well as the underlying mechanism: we identified miR-4319 as a tumor suppressor gene in the maintenance of CSC stemness and confirmed that miR-4319 could inhibit the tumor initiation and metastasis of CSCs mainly through downregulation of E2F2, which afforded the database to help understand the aggressive phenotype of TNBC. However, the retrospective analysis indicated that the stemness of TNBC-CSCs may be regulated by a series of miRNAs, and the correlation between miRNAs network and TNBC-CSCs requires further confirmation. In conclusion, miR-4319 plays an important role in regulating the self-renewal and tumorigenicity of stem cells and might be a remarkable prognostic factor or therapeutic target for patients with TNBC.

This study was financially supported by the National Natural Science Foundation of China (81172503, 81302304, 81572603, and 81570804), Key projects of Jiangsu Health Department (RC201153), the Second Level of 333 Project of Jiangsu Province (BRA2017534), the Third Level of 333 Project of Jiangsu Province (BRA2015470), the Project of Jiangsu Provincial Medical Talent (ZDRCA2016023), the Project of Jiangsu Provincial Medical youth Talent (QNRC2016095), General Project of Nanjing Health Bureau (YKK17179), and the Project of Wu Jieping Medical Foundation (320.6750.17006).

The miRNA-seq data in this paper are deposited to GEO (http://www.ncbi.nlm.nih.gov/geo/) with accession number GEO: GSE97765.

The authors declare that there is no conflict of interests.

J. Chu, Y. Li and X. Fan contributed equally to this work.