Characteristics and Expression Pattern of MYC Genes in Triticum aestivum, Oryza sativa, and Brachypodium distachyon

Myelocytomatosis oncogenes (MYC) transcription factors (TFs) belong to basic helix-loop-helix (bHLH) TF family and have a special bHLH_MYC_N domain in the N-terminal region. Presently, there is no detailed and systematic analysis of MYC TFs in wheat, rice, and Brachypodium distachyon. In this study, 26 TaMYC, 7 OsMYC, and 7 BdMYC TFs were identified and their features were characterized. Firstly, they contain a JAZ interaction domain (JID) and a putative transcriptional activation domain (TAD) in the bHLH_MYC_N region and a BhlH region in the C-terminal region. In some cases, the bHLH region is followed by a leucine zipper region; secondly, they display tissue-specific expression patterns: wheat MYC genes are mainly expressed in leaves, rice MYC genes are highly expressed in stems, and B. distachyon MYC genes are mainly expressed in inflorescences. In addition, three types of cis-elements, including plant development/growth-related, hormone-related, and abiotic stresses-related were identified in different MYC gene promoters. In combination with the previous studies, these results indicate that MYC TFs mainly function in growth and development, as well as in response to stresses. This study laid a foundation for the further functional elucidation of MYC genes.


Introduction
bHLH (basic helix-loop-helix) family is the second largest plant transcription factor (TF) family. They are characterized by a bHLH domain, which consists of 50-60 amino acids that form two distinctive regions, the basic region and HLH (helix-loop-helix) region [1]. The basic region functions as a DNA-binding motif, while the HLH region forms two amphipathic α helices with a linking loop [1][2][3]. In the bHLH family, there are some special members, named myelocytomatosis oncogenes (MYC). These members are characterized by a so-called bHLH_MYC_N region in the N-terminal and a bHLH region in the C-terminal [1,2,4,5].
MYC genes function in various physiological and molecular processes, especially in growth and development. In Arabidopsis thaliana, MYC2 functions synergistically with MYC3 and MYC4 in regulating leaf senescence [6], root elongation [7], stamen development [8], seed production, and seed storage protein accumulation [8,9], as well as chlorophyll degradation [10]; other Arabidopsis MYC genes, like ALCATRAZ functions in cell separation in fruit dehiscence [11], SPATULA controls development of carpel margin tissues [12], and ABORTED MICROSPORES (AMS) plays a crucial role in tapetum cell development and pollen wall formation [13]. In rice (Oryza sativa), the orthologous AMS gene Tapetum Degeneration Retardation (TDR) is necessary for tapetum degradation and anther development [14]; OsMYC2 is expressed in all tissues, and highly expressed in the spikelets and floral
Except for TaMYC1 and OsMYC2, other 25 TaMYC genes were named as TaMYC2A to TaMYC12B according to their distribution on chromosomes and genomic homology; other 6 OsMYC and 7 BdMYC genes were named as OsMYC1, OsMYC3 to OsMYC7, and BdMYC1 to BdMYC7 based on their chromosomal localization (Figure 1).
The physical features of these MYC TFs were predicted. In wheat, the protein length varies from 456 (TaMYC7A) to 699 (TaMYC4D) amino acids; the PI (Isoelectric Point) varies from 4.96 (TaMYC6D) to 8.73 (TaMYC9D); and the molecular weight varies from 49.17 kDa (TaMYC9A) to 63.69 kDa (TaMYC7D). In rice, the protein length varies from 473 (OsMYC7) to 904 (OsMYC4) amino acids; the PI varies from 4.68 (OsMYC5) to 8.36 (OsMYC3); and the molecular weight varies from 50.67 kDa (OsMYC3) to 97.39 kDa (OsMYC4). In B. distachyon, the protein length varies from 470 (BdMYC5) to 706 (BdMYC7); the PI varies from 4.53 (BdMYC6) to 7.22 (BdMYC5); and the molecular weight varies from 50.43 kDa (BdMYC5) to 76.24 kDa (BdMYC7). The protein subcellular localization was predicted by CELLO web server [33] and results showed that 26 MYC proteins locate in the nucleus, while other 14 MYC proteins locate in the chloroplast, cytoplasm, mitochondria. The detailed information is listed in Table S1.   [33] and results showed that 26 MYC proteins locate in the nucleus, while other 14 MYC proteins locate in the chloroplast, cytoplasm, mitochondria. The detailed information is listed in Table S1.

Sequence Alignment and Phylogenetic Tree of MYC TFs
Based on previous studies [13,18,[34][35][36][37][38], protein sequences of 13 known animal MYCs and 8 known Arabidopsis MYCs obtained from GenBank [39], and sequences of 40 MYCs in this study were used to perform sequence alignment. Multiple sequence alignment identified a distinct N-terminal bHLH_MYC_N region ( Figure S1). As shown in Table S2, the consensus ratio of 56 conserved amino acid residues is more than 50% in the bHLH_MYC_N domains. Leu-100, Trp-113, Tyr-115, Trp-119, and Leu-147 are identical in these 40 MYC TFs, indicating that these amino acids are key component of bHLH_MYC_N domain. More importantly, one JAZ interaction domain (JID) and one putative transcriptional activation domain (TAD) were identified in the bHLH_MYC_N domain. Furthermore, a bHLH domain and a leucine zipper region were found in the C-terminal region ( Figure 2 and Figure  S1).

Sequence Alignment and Phylogenetic Tree of MYC TFs
Based on previous studies [13,18,[34][35][36][37][38], protein sequences of 13 known animal MYCs and 8 known Arabidopsis MYCs obtained from GenBank [39], and sequences of 40 MYCs in this study were used to perform sequence alignment. Multiple sequence alignment identified a distinct N-terminal bHLH_MYC_N region ( Figure S1). As shown in Table S2, the consensus ratio of 56 conserved amino acid residues is more than 50% in the bHLH_MYC_N domains. Leu-100, Trp-113, Tyr-115, Trp-119, and Leu-147 are identical in these 40 MYC TFs, indicating that these amino acids are key component of bHLH_MYC_N domain. More importantly, one JAZ interaction domain (JID) and one putative transcriptional activation domain (TAD) were identified in the bHLH_MYC_N domain. Furthermore, a bHLH domain and a leucine zipper region were found in the C-terminal region ( Figure 2 and Figure S1).
In animals, MYC proteins are key regulators of mammalian cell proliferation that activate genes as part of a heterodimeric complex with the protein Max [40]. As shown in Figure S2 According to the alignment of bHLH domains in rice and Arabidopsis bHLH proteins [41], some amino acids residues, such as Val-960 and Ala-962 in the basic region, Asn-972 and Val-981 in the first Helix, Ser-986, Lys-987, Met-1006, and Asp-1008 in the Loop, and Asp-1017 in the second Helix, are only found in plant MYC proteins ( Figure 2 and Table S3). In addition, although the bHLH-ZIP domain is commonly present in some human c-Myc and Max proteins [42], such as myc_H. sapiens, myc_B. taurus, c-myc_M. musculus, n-Myc_S. scrofa in Figure 2, they were lowly conserved in the leucine zipper region.
To understand the evolutionary relationships of MYC TFs, a neighbor-joining (NJ) phylogenetic tree was constructed based on the full-length alignment of 40 putative MYCs, 13 known animal MYCs, and 8 known Arabidopsis MYC. As shown in Figure 3, these MYC TFs are classified into six sub-groups with highly bootstrap values. The class V is the largest while the class II was the smallest. There are 9 TaMYCs, 2 OsMYCs, and 2 AtMYCs in class I; 1 OsMYC, 1 BdMYC, and 1 AtMYC in class II; 3 TaMYCs, 1 OsMYC, and 1 BdMYC in class III; 2 TaMYCs, 2 BdMYCs, 1 AtMYC in class IV; 12 TaMYCs, 3 OsMYCs, 3 BdMYCs, and 4 AtMYCs in class V; and 13 known animal MYCs in class "animal". In animals, MYC proteins are key regulators of mammalian cell proliferation that activate genes as part of a heterodimeric complex with the protein Max [40]. As shown in Figure  According to the alignment of bHLH domains in rice and Arabidopsis bHLH proteins [41], some amino acids residues, such as Val-960 and Ala-962 in the basic region, Asn-972 and Val-981 in the first Helix, Ser-986, Lys-987, Met-1006, and Asp-1008 in the Loop, and Asp-1017 in the second Helix, are only found in plant MYC proteins ( Figure 2 and Table S3). In addition, although the bHLH-ZIP domain is commonly present in some human c-Myc and Max proteins [42], such as myc_H. sapiens, myc_B. taurus, c-myc_M. musculus, n-Myc_S. scrofa in Figure 2, they were lowly conserved in the leucine zipper region.
To understand the evolutionary relationships of MYC TFs, a neighbor-joining (NJ) phylogenetic tree was constructed based on the full-length alignment of 40 putative MYCs, 13 known animal MYCs, and 8 known Arabidopsis MYC. As shown in Figure 3, these MYC TFs are classified into six subgroups with highly bootstrap values. The class V is the largest while the class II was the smallest. There are 9 TaMYCs, 2 OsMYCs, and 2 AtMYCs in class I; 1 OsMYC, 1 BdMYC, and 1 AtMYC in class II; 3 TaMYCs, 1 OsMYC, and 1 BdMYC in class III; 2 TaMYCs, 2 BdMYCs, 1 AtMYC in class IV; 12 TaMYCs, 3 OsMYCs, 3 BdMYCs, and 4 AtMYCs in class V; and 13 known animal MYCs in class "animal".

Gene Structures and Conserved Motifs
We used genomic DNA sequences and CDS (coding sequences) to analyze the gene structure ( Figure 4B). Gene structures of TaMYC, OsMYC, and BdMYC genes are similar within the same subgroup. In wheat, the exon number ranges from 1 to 11. In rice and B. distachyon, the exon number ranges from 1 to 10. Notably, most class V members only have one exon.

Gene Structures and Conserved Motifs
We used genomic DNA sequences and CDS (coding sequences) to analyze the gene structure ( Figure 4B). Gene structures of TaMYC, OsMYC, and BdMYC genes are similar within the same subgroup. In wheat, the exon number ranges from 1 to 11. In rice and B. distachyon, the exon number ranges from 1 to 10. Notably, most class V members only have one exon.
Conserved motifs are helpful to understand the functions of MYC TFs. In this study, 8 conserved motifs were identified ( Figure 4C). Motifs 2, 4, 6, 7, and 8 constitute the bHLH_MYC_N domain ( Figure 4D). Among these five motifs, motifs 6 and 8 form JID and motif 2 composes TAD. JID consists of approximately 90 amino acids, including a specific motif (W- TAD includes approximately 70 amino acids and a specific motif In addition, as MYCs belong to bHLH TF families [3,41,43], motifs 1 and 5 constitute bHLH domain ( Figure 4D), and motif 3 constitutes the leucine zipper in class I, II, and V members.  The tree was constructed with 1000 bootstrap replications using MEGA7 based on the full-length protein sequence. The exon-intron structure of these genes was graphically displayed by the Gene Structure Display Server using the coding sequences (CDS) and DNA sequences of MYC genes. The protein sequence of MYC proteins was used to predict the conserved motifs/region by using the MEME Suite web server.

Synteny and Homologous Gene Pairs
Gene duplication events include tandem duplication and segmental duplication. We analyzed the gene duplication by using the MCScanX software [44]. As shown in Table S4, only one segmental duplication pair (TaMYC6D-TaMYC7D) was found.
In addition, orthologous between MYC genes in wheat, B. distachyon, A. thaliana, rice, H. vulgare, The tree was constructed with 1000 bootstrap replications using MEGA7 based on the full-length protein sequence. The exon-intron structure of these genes was graphically displayed by the Gene Structure Display Server using the coding sequences (CDS) and DNA sequences of MYC genes. The protein sequence of MYC proteins was used to predict the conserved motifs/region by using the MEME Suite web server.

Synteny and Homologous Gene Pairs
Gene duplication events include tandem duplication and segmental duplication. We analyzed the gene duplication by using the MCScanX software [44]. As shown in Table S4, only one segmental duplication pair (TaMYC6D-TaMYC7D) was found.
Additionally, we also analyzed the expression of these genes in two-week seedlings with different treatments (Figure 7). Although the results showed no regularity, there are some meaningful findings. In wheat, the expression of TaMYC1/5B/5D is induced by GA (Gibberellin), indicating the involvement in GA signal. TaMYC7A, TaMYC7B, and TaMYC7D displayed different expression patterns, implying that their functions are diversified. In rice, the expression of six expressed genes is upregulated by different abiotic stresses, suggesting their important roles. In B. distachyon, BdMYC3 is drastically induced by heat.
Additionally, we also analyzed the expression of these genes in two-week seedlings with different treatments (Figure 7). Although the results showed no regularity, there are some meaningful   proteins in plants [43,65], bHLH-ZIP domain only appeared on group B which MYC and Max exist. Compared to bHLH genes, MYC genes have three special kinds of cis-elements in their gene promoters. These cis-elements associate with their special functions. Moreover, MYC genes show tissue-specific expression patterns. Additionally, different from other non-MYC-bHLH genes which have 2 or more exons [41,43,65], most class V MYC members, which is the largest class only, have one exon.

Functions of MYC TFs
Previous studies showed that MYC TFs function in plant development and growth [2,9,37]. For example, MYC proteins function as regulators in regulating plant seed production, root elongation, leaf senescence, and stamen development [6][7][8][9]; they also regulate plant secondary metabolism and are actively involved in hormone-mediated plant growth [21,66,67]. In this study, most MYC genes are expressed in roots, stems, leaves, and inflorescences ( Figure 6). In combination with the identification of many cis-elements related to plant growth and development and hormone stress, these results further suggest their functions in growth and development.
The expression profiles of many MYC genes indicate their probable functions in response to abiotic and biotic stresses. For example, Arabidopsis MYC2, MYC3, MYC4, and MYC5 are induced by jasmonate [17,27,29,30]; MYC3 and MYC4 act additively with MYC2 to regulate defense against insect herbivory [30]; and MYC5 regulates JA-mediated plant defenses against herbivores and is involved in JA-regulated plant resistance to pathogens [37]. Overexpression of OsMYC2 results in bacterial blight resistance in rice [16]. In maize, the expression of two MYC genes, ZmbHLH103 and ZmbHLH104, is significantly upregulated under drought stress [65]. In this study, the expression of 6 OsMYC genes is induced by different abiotic stresses (Figure 7), further indicating that MYC genes

The Characteristics of MYCs
bHLH TFs are the second largest class of plant TFs [61]. Plant MYC TFs belong to the bHLH superfamily. MYC proteins have one bHLH_MYC_N domain in the N-terminal region [1]. In this domain, Leu-100, Trp-113, Tyr-115, Trp-119, and Leu-147 amino acids residues are complete conserved in all TaMYC, BdMYC, and OsMYC proteins. The bHLH_MYC_N include two domains: JID and TAD. The former is necessary for interacting with JAZ proteins, while the latter is a putative transcriptional activation domain [2]. For example, Arabidopsis MYC2, MYC3, and MYC4 interact with the C-terminal JAS domain of JAZ proteins through JID, and AtMYC2 recruits the mediator complex required for transcription initiation through its TAD, which specifically interacts with the activator interaction domain [62,63]. We also identified one specific motif (

W-[TN]-Y-[AG]-[IVL]-[FYL]-W-X(6,19)-L-[GT]-W-[GK]-[DE]-G) in JID and one specific motif ([VL]-[TG]-[DEG]-[TA]-E-[WML]-[FY]-[FY]-X(2)-[SC]-[MA]-X(3)-F-X(4)-G-[LAG]-P-G-X(9)-W) in TAD.
The C-terminal region of TaMYCs, BdMYCs, and OsMYCs is conserved and includes the typical basic region and HLH domain, and most MYCs contain leucine zipper (bZIP) domain (Figure 2 and Figure S2). Previous reports showed that the basic region of MYC proteins mainly recognizes CACGTG (G-box) sequence and CATGTG sequence, both E-box DNA-binding sites (5 -CANNTG-3 ), and the bHLH domain are required for the heterodimer to bind to the G-box in target genes [2]. The bHLH-ZIP domain is present in human c-Myc and Max proteins [42], which have DNA-binding activity and has been predicted to mediate protein-protein interactions. The interaction of Max and c-Myc depends on the integrity of the c-Myc HLH-Zip domain, but not the basic region or other sequences outside the domain [42]. Different from animal MYCs, plants MYCs contain a longer (about 70 amino acids) leucine zipper. Furthermore, the C-terminal leucine zipper domain can form the dimer that affects the specificity of interaction with other TFs [64]. Compared to other bHLH proteins in plants [43,65], bHLH-ZIP domain only appeared on group B which MYC and Max exist.
Compared to bHLH genes, MYC genes have three special kinds of cis-elements in their gene promoters. These cis-elements associate with their special functions. Moreover, MYC genes show tissue-specific expression patterns. Additionally, different from other non-MYC-bHLH genes which have 2 or more exons [41,43,65], most class V MYC members, which is the largest class only, have one exon.

Functions of MYC TFs
Previous studies showed that MYC TFs function in plant development and growth [2,9,37]. For example, MYC proteins function as regulators in regulating plant seed production, root elongation, leaf senescence, and stamen development [6][7][8][9]; they also regulate plant secondary metabolism and are actively involved in hormone-mediated plant growth [21,66,67]. In this study, most MYC genes are expressed in roots, stems, leaves, and inflorescences ( Figure 6). In combination with the identification of many cis-elements related to plant growth and development and hormone stress, these results further suggest their functions in growth and development.
The expression profiles of many MYC genes indicate their probable functions in response to abiotic and biotic stresses. For example, Arabidopsis MYC2, MYC3, MYC4, and MYC5 are induced by jasmonate [17,27,29,30]; MYC3 and MYC4 act additively with MYC2 to regulate defense against insect herbivory [30]; and MYC5 regulates JA-mediated plant defenses against herbivores and is involved in JA-regulated plant resistance to pathogens [37]. Overexpression of OsMYC2 results in bacterial blight resistance in rice [16]. In maize, the expression of two MYC genes, ZmbHLH103 and ZmbHLH104, is significantly upregulated under drought stress [65]. In this study, the expression of 6 OsMYC genes is induced by different abiotic stresses (Figure 7), further indicating that MYC genes participate in response to abiotic stresses. Consistently, many cis-elements related to abiotic stresses were identified.
Taken together, we can draw a conclusion that MYC TFs mainly function in growth/development, and in response to environmental stresses.

Identification of MYC TFs in Wheat, Rice, and B. distachyon
The wheat, rice, and B. distachyon genome sequences, protein sequences, coding sequences (CDS), and upstream 2-kb genomic DNA sequences were downloaded from Ensembl plants [68]. The chromosomal distribution of MYC genes was obtained from the wheat, rice, and B. distachyon genome annotations in Ensembl plants [68]. To identify the MYC TFs in wheat, rice, and B. distachyon, 7 known MYC proteins, including wheat TaMYC1, barley HvMyc1, rice OsMYC2, Arabidopsis MYC2, MYC3, MYC4, and AMS proteins sequences, were downloaded from Ensembl plants [68] and were used to build a hidden Markov model (HMM), then searched against the genome protein sequences of wheat, rice, B. distachyon, respectively, with a threshold of E< 1e-5 and the length of amino acid >200 aa (amino acid). Then, the HMM profile of the bHLH-MYC_N domain (PF14215) was downloaded from the Pfam database [69] to search against the protein sequences of wheat, rice, and B. distachyon with a threshold of E < 1e-5. After manual correction to remove the redundancy and alternative splice, the NCBI-CDD database (NCBI Conserved Domains Database) and SMART database (Simple Modular Architecture Research Tool) were used to confirm the putative MYC proteins. The ExPASy webserver [70] was used to predict the theoretical isoelectric point and molecular weight of TaMYCs, OsMYCs, and BdMYCs. To further verify the existence of MYC genes in wheat, rice, and B. distachyon, we performed BLASTN [71] to search for ESTs using the CDS of MYC genes. The CELLO web server [33] was used to predict the subcellular localization of MYC proteins.

Multiple Sequence Alignment and Phylogenetic Analysis of MYCs
An unrooted neighbor-joining (NJ) tree was constructed by MEGA 7.0 software with 1000 bootstrap replications and Jones-Taylor-Thornton model based on full-length protein sequence alignments was done by using the ClustalX 2.0 method [72,73]. The multiple sequence alignment was performed using the ClustalX 2.0 and visualized by Jalview [74], and the phylogenetic tree was visualized by Evolview [75].

Analysis of Gene Structures and Conserved Motifs
The exon-intron structure of MYC genes was graphically displayed by the Gene Structure Display Server [76] by using the CDS and DNA sequences of wheat, rice, and B. distachyon MYC genes. The protein sequences of MYCs were used to predict the conserved motifs by using the MEME online program [77] with the following parameters: optimum width of motifs set from 5 to 200 amino acids and maximum number of motifs set at 8. The gene structure and motifs were visualized by Evolview [75].

Analysis of Cis-Elements, Gene Duplication, and Synteny
The 2-kb upstream genomic DNA sequences of MYC genes were submitted to the PlantCARE online tool to identify the cis-elements [44]. Gene duplication and systeny analysis of MYCs in different species were done by using the MCScanX software [43] and visualized by Dual Systeny Plotter software by CJ-Chen [78].

Expression Profile Analysis
The cultivar of T. aestivum 'Chinese Spring', rice cultivar 'Dongjin', and B. distachyon Bd-21 were planted in an artificial climate chamber at 26/22 • C (day/night) with a photoperiod of 16/8 h (day/night). For tissues analysis, roots, stems, leaves, and inflorescences were collected at the heading stage. For different abiotic stresses, 2-week-old seedling plants were subjected to H 2 O (CK), heat (42 • C), cold (4 • C), drought (20% PEG6000), salt (200 mM NaCl), ABA (100 µM), and GA (100 µM) for 2 h, then whole plants were collected for RNA isolation. Total RNA was extracted using RNAiso Reagent (TaKaRa, Dalian, China) according to the manufacturer's protocol and treated with DNaseI. The Transcriptor First Strand cDNA Synthesis Kit (Roche) was used to synthesize cDNA based on the manual. The QuantStudio TM Real-Time PCR Software (ThermoFisher Scientific) was used to carry out qRT-PCR, data acquisition, and analysis. The volume of each reaction was 15 µL containing 7.5 µL of SYBR Premix Ex Taq (TaKaRa, Dalian, China), 0.75 µL (10 pmol µL −1 ) each of forward and reverse primers, 0.5 µL of cDNA (5.0 ng µL −1 ), and 5.5 µL of ddH 2 O. The reference genes (Taactin, Osactin, Bdactin) [79][80][81] were used to normalize the expression of the TaMYC, OsMYC, and BdMYC genes and the 2 (−∆∆Ct) analysis method was used to determine the relative expression level [82]. The following program was used for qRT-PCR: 50 • C for 2 min, 95 • C for 10 min, followed by 40 cycles of 95 • C for 15 s, 60 • C for 1 min in the PCR stage and 95 • C for 1 min, 95 • C for 15 s in the melt curve stage. Primers for 6 sets of TaMYC genes (TaMYC3A/B/D, TaMYC4A/B/D, TaMYC1/5B/5D, TaMYC8A/B/D, TaMYC9A/B/D, and TaMYC10A/B/D) were universal in each set, because of the highly conserved sequences in A, B, and D sub-genomes, thus the detected expression is the combination of three copies of homologous genes. The primers are listed in Table S6.

Conclusions
In this study, 26 TaMYC, 7 OsMYC, and 7 BdMYCTFs were identified and divided into five groups. The MYC proteins contained a JAZ interaction domain (JID) and a putative transcriptional activation domain (TAD) in the bHLH_MYC_N domain of the N-terminal region, the bHLH domain was found in the C-terminal region, some MYCs followed by a bHLH-ZIP domain. The expression profiles and identified three kinds of cis-elements indicate that MYC TFs function in plant development and in response to environmental stresses. Taken together, our results provide a solid foundation for further structural and functional investigations on MYCs.
Supplementary Materials: The following data are available online at http://www.mdpi.com/2223-7747/8/8/274/s1, Figure S1: The JAZ interaction domain (JID), putative transcriptional activation domain (TAD) were identified in the MYC proteins; Figure S2: The bHLH domain and leucine zipper domain were identified in the MYC proteins; Table S1: The detailed information of TaMYCs, OsMYCs, BdMYCs; Table S2: The consensus ratio of the conserved amino acid residues in the bHLH_MYC_N domains, Table S3: The consensus ratio of the conserved amino acid residues in the bHLH domains of TaMYCs, OsMYCs, and BdMYCs, Table S4: Wheat, rice, and B. distachyon MYC gene duplication events, Table S5: The homologous genes between plant MYC genes, Table S6: Primer used in this study.
Author Contributions: H.L. conceived and designed the project; S.C., H.Z., T.L., and X.N. performed main experiments; H.L. and S.C. analyzed the data. H.L. and S.C. wrote the manuscript. All authors have read and approved the manuscript.
Funding: This work was supported by the National Natural Science Foundation of China (31571657). The funding body did not exert influence on the design of the study, and collection, analysis, and interpretation of data, or in writing of the manuscript.

Conflicts of Interest:
The authors have declared that they have no conflicts of interest.