Hu Tao
Juan Yang
Pengpeng Zhang
Nian Zhang1
Mingxin Chen- 1Hubei Key Laboratory of Animal Embryo Engineering and Molecular Breeding, Institute of Animal Husbandry and Veterinary, Hubei Academy of Agricultural Sciences, Wuhan, China
- 2Department of Biotechnology, College of Life Sciences, Xinyang Normal University, Xinyang, China
Long noncoding RNAs (lncRNAs) were identified recently as a large class of noncoding RNAs (ncRNAs) with a length ≥200 base pairs (bp). The function and mechanism of lncRNAs have been reported in a growing number of species and tissues. In contrast, the regulatory mechanism of lncRNAs in the goat reproductive system has rarely been reported. In the present study, we sequenced and analyzed the lncRNAs using bioinformatics to identify their expression profiles. As a result, 895 lncRNAs were predicted in the pre-ovulatory ovarian follicles of goats. Eighty-eight lncRNAs were differentially expressed in the Macheng black goat when compared with Boer goat. In addition, the lncRNA XR_311113.2 acted as a sponge of chi-miR-424-5p, as assessed via a luciferase activity assay. Taken together, our findings demonstrate that lncRNAs have potential effects in the ovarian follicles of goats and may represent a promising new research field to understand follicular development.
Introduction
Ovulation rate is an important goat reproductive trait that determines the upper limit of the female goat litter size; in turn, the growth and development of follicles in the ovary determine the ovulation rate (Cui et al., 2009). Follicle development is a very complex biological process and its regulatory mechanism warrants further study. During mammalian ovarian folliculogenesis, which involves multiple transcription factors, primordial follicles develop into pre-ovulatory follicles, followed by ovulation, which releases mature oocytes (Choi and Rajkovic, 2006; Pan et al., 2012). The Macheng black goat is an excellent goat breed that is unique to China and is characterized by pure black hair color and its breed function is higher than Boer goat (Tao et al., 2018). In the breeding production process, PMSG-hCG is usually used to treat the ewe, and then artificial fertilization. We found that after hormone treatment, the litter size of Macheng black goat remains above Boer goat. To investigate the phenotypic differences, we analyze the pre-ovulatory follicles of ewes after hormone treatment via High-throughput sequencing and analyze the development of follicles in a highly fertile local goat.
Noncoding RNAs (ncRNAs) have always been regarded as “DNA junk” without the capability to encode proteins. More recently ncRNAs, which include mainly microRNAs (miRNAs), piwi-interacting RNAs (piRNAs), and circular RNAs (circRNAs), were demonstrated to have important biological functions. They are basically derived from the transcription activity of the genome (Jin et al., 2018). Long noncoding RNAs (lncRNAs), which were recently identified as a large class of ncRNAs, are defined as ncRNAs of more than 200 bp in length with no protein-coding capacity (Mercer et al., 2009; Hung and Chang, 2010). A recent study showed that lncRNAs play multifunctional roles in a wide variety of important biological processes, including cell proliferation (Wang et al., 2017), cell apoptosis (Lu et al., 2016), genomic splicing (Parikshak et al., 2016), chromatin remodeling (Han and Chang, 2015), and transcription (McHugh et al., 2015). In addition, lncRNAs regulate the expression of miRNA by acting as molecular sponges and reducing the inhibitory effect of miRNAs on their target genes (Paraskevopoulou and Hatzigeorgiou, 2016).
During ovarian folliculogenesis, functional lncRNAs, such as growth arrest specific 5 (GAS5) (Wang et al., 2018), lncRNA HCG26 (Liu et al., 2017), lncRNA Neat1 (Nakagawa et al., 2014; Li et al., 2017), and lncRNA MALAT1 (Li et al., 2018a), have been found to regulate follicle development and regeneration through diverse mechanisms. For instance, GAS5, which is expressed in female germ-line stem cells (FGSCs) and oocytes, promotes FGSC reproduction and survival in vitro and is highly expressed in neonatal mouse ovaries (Wang et al., 2018). High levels of lncRNA SRA stimulate the growth of mouse follicle granulosa cells, increase the levels of estrogen and progesterone, and upregulate the expression of key enzymes (i.e., cytochrome P450 family 19, subfamily A, member 1 (CYP19A1) and cytochrome P450 family 11, subfamily A, member 1 (CYP11A1)) (Li et al., 2018b). Although miRNAs and piRNAs have been extensively studied in ovarian folliculogenesis, the role of other valuable noncoding transcripts remains to be studied. Moreover, researchers are paying increasing attention to the role of lncRNAs in follicular development, and have used high-throughput sequencing to analyze the expression pattern of lncRNAs during follicular development. Finally, the identification of the key lncRNAs for follicular development is expected.
In the present study, differentially expressed lncRNAs were detected through deep RNA sequencing (RNA-seq) in samples of pre-ovulatory follicles from Macheng black goats and Boer goats. Eight hundred and ninety-five new lncRNAs were identified in the ovarian follicles of goats. Eighty-eight lncRNAs were differentially expressed between Macheng black and Boer goats. The lncRNA expression patterns were validated using quantitative real-time polymerase chain reaction (qRT-PCR). The greatly differentially expressed lncRNA, XR_311113.2 was demonstrated to sponge chi-miR-424-5p via a luciferase activity assay. Therefore, in this study we examined the expression of lncRNAs by RNA-seq during goat follicle development and explored their expression profiles and functions in goat follicles.
Materials and Methods
Ethics Statement
All studies involving animals were conducted according to the regulation (No. 5 proclaim of the Standing Committee of Hubei People’s Congress) approved by the Standing Committee of Hubei People’s Congress, P. R. China. Sample collection was approved by the ethics committee of Hubei Academy of Agricultural Sciences. Animals were humanely sacrificed as necessary to ameliorate suffering.
Animals and Tissues
Aged adult ewes weighing 40 kg were obtained from the goat stud farm of the Institute of Animal Husbandry and Veterinary Medicine, Hubei Academy of Agricultural Sciences, feeding with the same pattern and environment. Macheng black goats (Kidding rate and fecundity rate: 219 and 346%) have higher reproductive performance than Boer goats (Kidding and fecundity rates: 189 and 210%). Three Macheng black goats and three Boer goats exhibiting normal estrous cycles were treated with 1000 IU PMSG (SanSheng, Zhejiang, China) and 500 IU hCG (SanSheng, Zhejiang, China) as previously described (Saharrea et al., 1998; Torner et al., 1998). PMSG-hCG was used to initiate and synchronize the follicular phase. After 36 h, the PMSG-hCG stimulated ewes were slaughtered, and the ovaries of each goat in the pre-ovulatory phase were immediately removed, stored in liquid nitrogen until RNA extraction. The total RNAs were extracted from fourteen tissues (i.e., liver, spleen, heart, kidneys, small intestine, fat, ovarian folliclelungs, lungs, uterus, abomassum, muscle, reticulum, rumen and omasum) of the Boer goats.
RNA Extraction and Qualification
A total of 5 μg of RNA was isolated from each individual sample using the TRIzol reagent (Invitrogen, MD, United States) according to the manufacturer’s protocol. The purity and quantity of the total RNA were measured using a NanoDrop instrument and agarose gel electrophoresis. RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, United States), and only the samples showing RNA integrity number (RIN) scores higher than 6 were used in this study.
Library Preparation for lncRNA and mRNA Sequencing
Three micrograms (μg) of RNA was used as input material. Ribosomal RNA was completely removed using the Epicentre Ribo-zero™ rRNA Removal Kit (Epicentre, WI, United States), and then sequencing libraries were generated with the NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina® (NEB, Ispawich, United States) following manufacturer’s recommendations. First-strand cDNA was synthesized using M-MuLV Reverse transcriptase and random hexamer primers, and second-strand cDNA was synthesized using RNase H and DNA polymerase I. After adenylation of the 3′ ends of DNA fragments, NEB Next Adaptors were ligated to prepare for hybridization. The library fragments were purified for selecting cDNA fragments (150–200 bp in length). Finally, the products were purified with an AMPure XP system (Beckman Coulter, Beverly, United States). The libraries were sequenced on an Illumina Hiseq 4000 platform, and 150 bp paired-end reads were generated. Sequencing data were deposited with the NCBI Sequence Read Archive (SRA) under accession number PRJNA648013.
Quality Control and Transcriptome Assembly
Firstly, Raw reads were processed through in-house Perl scripts developed by the Novogene Bioinformatics Institute (Beijing, China). The clean data were obtained by removing reads that contain adapter or ploy-N and low-quality reads from raw data. Meanwhile, The Phred score (Q20, Q30) and GC content of the clean data were calculated. Reference genome and gene model annotation files were downloaded from National Center for Biotechnology Information (ARS1, https://www.ncbi.nlm.nih.gov/genome/?term=goat). Index of the reference genome was built by Bowtie v2.4.2 and paired-end clean reads were aligned to the reference genome using TopHat v2.1.1. The mapped reads of each sample were assembled by both Scripture (beta2) and Cufflinks (v2.2.1).
Coding Potential and Conservative Analysis
We used a Coding Potential Calculator (CPC, 0.9-r2) (Kong et al., 2007), Coding-Non-Coding-Index (CNCI, v2) (Sun et al., 2013) and PfamScan (v1.3) (Punta et al., 2012) to assess the protein-coding potential of each novel transcript. CPC examined sequences in NCBI eukaryotes’ protein database and set the e-value “1e-10.” Default parameters were used in CNCI profiles. Pfam searches use default parameters of -E 0.001 -domE 0.001. Candidate sets of lncRNAs that possessed no coding potential were chosen and filtered out from the predicted transcript results from the three tools listed.
Analysis of Differential Expression
Cuffdiff (v2.2.1) was used to estimate the fragment per kilobase of exon per million fragments mapped (FPKMs) of the transcripts in each sample (Trapnell et al., 2012). The unit of measurement is FPKM. Cuffdiff calculated for the level of expression in each transcript. Gene and transcript expressions based on FPKM values were calculated by Cufflinks transcript quantification engine. Transcripts with a p < 0.05 were assigned as differentially expressed.
Target Gene and miRNA Binding Site Prediction of lncRNAs
To explore the function of candidate lncRNAs, the target genes of candidate lncRNAs were predicted in the cis way. Coding genes 100 k upstream and downstream of lncRNAs were examined as the target genes in cis (co-location genes). The co-expression of coding genes with lncRNAs in different chromosomes were determined with the pearson correlation coefficient (PCC); with PCC >0.95 or < −0.95, the lncRNA-mRNA pair was considered to represent the target genes in trans (co-expressed genes) (Schober et al., 2018). MiRNA binding sites were predicted for the goats with miRanda software, whose principles are based on the miRanda prediction algorithm (Enright et al., 2003). The miRNAs and target lncRNAs were considered when the miRanda score was 140 or higher, and the energy threshold was set to -1.
Analysis of lncRNA-miRNA Network Interactions
The lncRNA-miRNA interaction network was built according to the prediction of miRNA binding sites. The lncRNA-miRNA interaction analysis was conducted with Cytoscape software (Version 3.6.1). In the network diagram, the connections indicate possible regulatory relationships. The square represents lncRNAs, the circle represents miRNAs, the red represents up-regulated expression and the green represents down-regulated expression.
Quantitative real‐time PCR Analysis
Total cDNA was synthesized using reverse transcriptase Kit (TaKaRa, Dalian, China). QRT-PCR were performed using CFX96 Touch™ Real-Time PCR Detection System and SYBR® Green PCR Supermixe (Bio-Rad, CA, United States). Each PCR reaction (in 20 μL) involved 10 μL SYBR® Green PCR Supermixe (Bio-Rad, CA, United States), 0.25 μL of each primer, 1 μL cDNA and 8.5 μL H2O. The cycling conditions included an initial single cycle (94°C for 3 min), and followed by 40 cycles (94°C for 30 s; 60°C for 30 s; 72°C for 20 s). The primers were designed using Primer 5 (shown in Table 1). The expression level was defined based on the threshold cycle (Ct), and relative expression levels were calculated via the 2−ΔΔCt method. The correlation between the RNA-seq and qRT-PCR results was calculated using the SPSS (Version 18.0.0). β-actin was served as internal standard control, and all reactions were performed in triplicate.
TABLE 1. The nucleotide sequence used in this study.
Plasmid Construction
Two fragments of XR_311113.2 sequence were amplified using I-5™ 2 × High-Fidelity Master Mix (Tsingke, Wuhan, China). All the PCR products were inserted into pmirGLO vector (Promega, WI, United States), respectively. The primers are listed in Table 1. Then two recombinant plasmids were digested with PmeI and NheI (Thermo Scientific, WLM, United States). All constructs were sequenced by Sangon Biotech Co., Ltd. (Shanghai, China). The mutants of binding sites were generated using a MutanBEST Kit (TaKaRa, Dalian, China) and mutagenic primers (shown in Table 1).
Cell Culture, Cell Transfection, and Dual Luciferase Reporter Assays
The goat kidney epithelial cells (GK cells) were obtained from the Kunming cell bank of the Chinese academy of sciences. Cells were seeded at a density of 1.5 × 105/ml using Dulbecco minimum essential medium (DMEM) (Hyclone, UT, United States) supplied with 10% fetal bovine serum (FBS) medium (Hyclone, UT, United States). The miRNA mimics and negative control (NC) were synthesized by Ribobio (Ribobio, Guangzhou, China). The sequences used are shown in Table 1. Cells were seeded in 24-well plates at a density of 1.5 × 105 cells per well 24 h before transfection. The cells were co-transfected with a mixture of 250 ng pmirGLO recombinant vector plasmid and 3 μL miRNA mimics or NC per wells using lipofectamine 3000 (Life Technologies, MD, United States). The luciferase activity in the transfected cells was then measured after 24 h according to the manufacturer’s instructions with a VICTOR™ X3 PerkinElmer 2030 Multilabel Plate Reader (PerkinElmer, MA, United States).
Statistical Analyses
Descriptive results of the study are expressed as means ± SD. All data analyses were performed by GraphPad Prism 5 (San Diego, CA). The means of the groups were compared by Student’s t‐test. p < 0.05 was considered to indicate statistically.
Results
Identification of lncRNAs in Goat Pre-ovulatory Follicles
To identify the expression of lncRNAs in goat pre-ovulatory follicles, six goats average from two goat breeds, Boer and Macheng goats, were sequenced using an Illumina HiSeq 4000 platform. A total of 41.3 Gb of lncRNA clean data were obtained and mapped to the goat reference genome using TopHat2. We carried out a series of rigorous screening analyses using three analytical software programs (CPC, PFAM, and CNCI) (Figures 1A,B), and a total of 895 lncRNAs from 99,106 assembled transcripts were identified (Figure 1C and Supplementary Table S1). These lncRNAs consisted of 88.0% long intergenic noncoding RNAs (lincRNAs) and 12.0% antisense lncRNAs, with virtually no intronic lncRNAs being detected (Figure 1D). We analyzed the characteristics and expression levels of the lncRNA and mRNA transcripts. As shown in Figure 1E, we found that the predicted lncRNAs were shorter in length than the mRNAs and had fewer exons than the average mRNA. Furthermore, most of the lncRNAs tended to exhibit a relatively shorter open reading frame (ORF) length than the mRNAs. The mRNA levels were significantly higher than those of lncRNAs or transcripts of uncertain coding potential (TUCP) in the pre-ovulatory follicles of goat samples (Figures 1F–H).
FIGURE 1. Characterization of pre-ovulatory follicles for goat lncRNAs, mRNAs, and TUCPs. (A,B). Steps of the novel lncRNA screening. (C) Coding potency filter with three mainstream coding potential analytical methods (CPC, CNCI, and FFAM). (D) Classification of the lncRNAs identified here. (E) Transcript size distribution, ORF length, and exon number of lncRNA and mRNA transcripts. (F−H) Box plot, violin plot, and density distribution diagram showing the expression features of lncRNAs, mRNAs, and TUCPs in the pre-ovulatory follicles of goats.
Expression Profiles of lncRNAs and mRNAs in Goat Follicles
Here, we identified the differential expression profiles of lncRNAs and mRNAs in the pre-ovulatory follicles of Boer and Macheng Black goats. Specifically, a total of 88 lncRNAs (67 upregulated and 21 downregulated) and 1216 mRNAs (743 upregulated and 473 downregulated) were differentially expressed (Figures 2A,B and Supplementary Tables S2, S3). A hierarchical clustering analysis was performed on the differentially expressed lncRNAs and mRNAs, respectively. To gain insight into the similarities of the pre-ovulatory follicles between the Boer and Macheng goats, data from all of the differentially expressed lncRNAs and mRNAs were used in a systematic cluster analysis. The heat map clearly indicated self-segregated clusters in the Boer (BO) and Macheng (MC) groups (Figures 2C,D).
FIGURE 2. Differential expression of lncRNAs and mRNAs in the pre-ovulatory follicles of goats. (A,B) Volcano plots of differentially expressed lncRNA and mRNA transcripts. (C,D) Hierarchical clustering of the expression profiles of differentially expressed lncRNAs and mRNAs.
Overview of Differential lncRNA Expression Profiles
Five upregulated (XR_001297559.1, XR_311113.2, XR_001297560.1, LNC_000026, and XR_310768.2) and two downregulated (LNC_000155 and XR_001295597.1) lncRNAs were randomly selected and amplified via qRT-PCR to confirm RNA-seq accuracy. As shown in Figure 3A, log2-fold changes (MC/BO) were calculated based on the RNA-seq and qRT-PCR results. The observed expression trends indicated that the two methods produced consistent results. Furthermore, the fold change values measured by RNA-seq and qRT-PCR were significantly correlated (correlation coefficient, 0.944) (Figure 3B). Moreover, we investigated the relative expression levels of the greatly differentially expressed lncRNA XR_311113.2 in 14 goat tissues of Macheng black goat and Boer goat (i.e., liver, spleen, heart, kidneys, small intestine, fat, ovarian follicle, lungs, uterus, abomasum, muscle, reticulum, rumen, and omasum) (Figure 3C). XR_311113.2 tended to exhibit high expression levels in the liver, spleen, heart, and kidneys of Macheng black goat and Boer goat, followed by the small intestine, fat, ovarian follicle, and lungs. By contrast, its low expression was noted in the uterus, abomasum, muscle, reticulum, rumen, and omasum. Notably, the expression level of XR_311113.2 in Macheng ovarian follicle was far greater than Boer, and the fold change was nearly 6. These results demonstrated that there was a general consistency between the qRT-PCR and RNA-seq results, although the fold changes were not exactly the same between the two different technologies. In addition, XR_311113.2 was only significantly dimorphically expressed in ovarian follicle.
FIGURE 3. Validation of selected lncRNAs and mRNAs using qRT-PCR. (A) Fold changes in the relative expression of lncRNAs, as assessed by qRT-PCR. Results of the comparison of the seven lncRNAs using qRT-PCR and RNA‐seq. The vertical axis indicates the mean fold change (log2-fold change) of each lncRNA. (B) Correlation analysis of the fold changes between qRT-PCR and RNA-seq. (C) Relative expression of lncRNAs. Expression profile of lncRNA XR_311113.2 in 14 tissues of Macheng and Boer goat expressed as the mean ± SD. Data were normalized to the reference gene (β-actin).
lncRNA-miRNA Interaction Networks
The lncRNA-miRNA interaction networks were constructed using the RNA-seq data and potential lncRNA-miRNA connections were explored using the Cytoscape 3.6.1 software (http://www.cytoscape.org/). The lncRNA-miRNA interaction network was established based on the relationship between the differentially expressed (DE) lncRNAs and DE miRNAs. In the network depicted in Figure 4A, we identified the top 20 lncRNAs and top 20 miRNAs as differentially expressed profiles. The top differentially expressed lncRNA XR_311113.2 included binding sites for four miRNAs (chi-miR-135a, chi-miR-424-5p, chi-miR-544-5p, and chi-miR-3955-5p). The interaction network of XR_311113.2 and these four miRNAs was predicted by miRanda (Figure 4B). This result implies that These results indicate that the aforementioned miRNA and lncRNA might have a tight correlation and regulation relationship, which could be our main focus for further study.
FIGURE 4. Biomathematically predicted XR_311113.2 targeted in the lncRNA-miRNA network. The networks of lncRNA-miRNA were indicated by Cytoscape. Red: upregulated gene; green: downregulated gene; circle: lncRNA; triangle: miRNA. (A) The interaction network of the top 20 differentially expressed lncRNAs and mRNAs. (B) Predicted lncRNA-miRNA network of XR_311113.2 with four miRNAs.
The lncRNAs may affect the expression of miRNAs by sponging them.
XR_311113.2 Serves as a Sponge for Chi-miR-424-5p
The miRanda software was used to predict the putative miRNA binding sites of XR_311113.2. Figure 5A indicates that four miRNAs (chi-miR-135a, chi-miR-424-5p, chi-miR-544-5p, and chi-miR-3955-5p) bind to XR_311113.2 via putative binding sites. To verify that these four miRNAs bind to XR_311113.2, the miRNA binding site regions and the mutated regions of chi-miR-135a (719 bp/741 bp), chi-miR-424-5p (975 bp/995 and 3197 bp/3216 bp), chi-miR-544-5p (680 bp/700 bp), and chi-miR-3955-5p (3569 bp/3591 bp) were cloned into two pmirGLO vectors, named pmirGLO-1 (159 bp/1163 bp, 1005 bp) and pmirGLO-2 (2104 bp/3628 bp, 1525 bp), respectively. Three miRNA mutants (chi-miR-544-5p mut, chi-miR-135a mut, and chi-miR-424-5p mut) of pmirGLO-1 and two miRNA mutants (chi-miR-424-5p mut and chi-miR-3955-5p mut) of pmirGLO-2 were transfected into GK cells. We discovered that chi-miR-424-5p mut (p < 0.001) and chi-miR-3955-5p mut (p < 0.01) completely abolished the suppression of luciferase activity compared with the wild-type vector (Figure 5B). These results suggest that chi-miR-424-5p and chi-miR-3955-5p bind to the putative miRNA binding site of XR_311113.2. To clarify further the binding activity of miRNAs, the mimic of chi-miR-424-5p and chi-miR-3955-5p was cotransfected with the pmirGLO-2 luciferase reporter plasmid into GK cells. Compared with the NC group, chi-miR-424-5p significantly reduced luciferase activity (p < 0.01) (Figure 5C). These results suggest that XR_311113.2 serves as a sponge of chi-miR-424-5p.
FIGURE 5. XR_311113.2 serves as a sponge for chi-miR-424-5p. (A) Alignment of potential chi-miR-135a, chi-miR-424-5p, chi-miR-544-5p, and chi-miR-3955-5p base pairing with XR_311113.2, as identified by the miRanda software. (B) The luciferase activity of miRNA binding sites mutant pmirGLO vector by luciferase assay in GK cells. Site-directed mutagenesis in the chi-miR-135a, chi-miR-424-5p, chi-miR-544-5p, and chi-miR-3955-5p binding site of the XR_311113.2. The pmirGLO-1 and pmirGLO-2 were wild vector of XR_311113.2. (C) Luciferase reporter assay to determine the luciferase activity of pmirGLO-2 in GK cells transfected with chi-miR-424-5p and chi-miR-3955-5p mimics, to identify miRNAs that bound to the XR_311113.2 sequence. **p < 0.01; ***p < 0.001; NS, not significant.
Discussion
Goat is a domestic animal with a long history. The Macheng black goat is a unique breed indigenous to mountainous areas of Central China Region and characterized by systemic black and high prolificacy. Boer goat is one of the most popular breeds of meat goat in the world due to their excellent carcass qualities and high growth efficiency (Naing et al., 2010). The reproductive performance of ewes is an important part of goat productivity. The kidding and fecundity rates of Boer goats are approximately 189% and 210% (Malan, 2000). The Macheng black goat is highly fertile, the kidding rate and fecundity rate are approximately 219% and 346% (Tao et al., 2018). Recent studies indicate that a high ovulation rate can increase the probability of higher litter size (Wang et al., 2020). As a high propagation capability goat breed, how lncRNAs regulate its follicular development is still unclear. Therefore, we investigated the differential expression of lncRNAs in ewe follicular development between Macheng goat and Boer goat using RNA-seq.
Folliculogenesis is a complex process that includes primordial follicle recruitment, granulosa cell (GC) proliferation, oocyte maturation, steroidogenesis, and ovulation (McLaughlin and McIver, 2009; Field et al., 2014). GCs play a crucial role in each stage of folliculogenesis, provide nutritional support for the development of oocytes through the gap junction, and enhance the maturation of oocytes. The interaction between GCs and theca cells is an important condition for maintaining follicular function and promoting normal development. The growth and differentiation processes of GCs are key criteria for the normal initiation and growth of primordial follicles, and mediate the development and atresia process of growth-period follicles, thereby playing an important regulatory role in the development of follicles (Uyar et al., 2013; Johnson, 2015; Sugiyama et al., 2016; Nguyen et al., 2019; Kim et al., 2020).
LncRNAs are newly discovered ncRNAs that have important regulatory functions in a range of biological events (Peng et al., 2017), especially their role as therapeutic targets and diagnostic biomarkers of various diseases (Ji et al., 2003; Yan et al., 2015; Huang et al., 2017; Huang, 2018; Simion et al., 2019; Yan et al., 2019). A growing body of evidence has shown that lncRNAs are widely distributed in mammals and plants (Necsulea et al., 2014; Wang et al., 2014; Zhao et al., 2018; Sarropoulos et al., 2019), but little has been reported regarding the folliculogenesis process of goats. In the field of molecular biology, the genetic research of goats is subject to many restrictions compared with the research of other species, such as humans and mice. The database of lncRNAs mainly targets common research subjects (humans, mice, etc). The lack of various databases can hamper the analysis of lncRNA data.
In the current study, we used high-throughput sequencing to analyze the lncRNA profiles in goat ovarian follicle samples. In addition, we investigated the manner in which lncRNA and mRNA networks contribute to follicle development, and confirmed the regulatory function of lncRNA XR_311113.2. We predicted a total of 895 lncRNAs in samples of pre-ovulatory follicles from Macheng black goats and Boer goats, which was less than that detected in the developmental skeletal muscle of fetal goat (3981) (Zhan et al., 2016) and the fetal skin of goats (1336) (Ren et al., 2016), and more than that detected in the anagen-phase skin samples of cashmere goats (437) (Zheng et al., 2019). These results revealed that the expression of lncRNAs was tissue-specific in different goat breeds.
A recent study has revealed that lncRNAs can function as miRNA sponges to further affect the expression of miRNA target genes (Du et al., 2016). For example, lncRNA TRPM2-AS acts as the sponge of miR-612 to promote gastric cancer progression (Xiao et al., 2020), and lncRNA MDNCR binds to miR-133a, thus promoting cell differentiation in bovine primary myoblasts (Li et al., 2018). In the present study, we found that lncRNA XR_311113.2 was differentially expressed in the pre-ovulatory follicles of the Boer and Macheng goat groups. The bioinformatics analyses showed that four miRNAs potentially interacted with XR_311113.2 (as shown in Figure 4). In addition, we validated the interacting relationship between XR_311113.2 and chi-miR-424-5p using luciferase activity assays (as shown in Figure 5). Therefore, our results suggest that XR_311113.2 targets chi-miR-424-5p directly by functioning as a sponge in the pre-ovulatory follicles of goats. Moreover, the previous study results demonstrated that miR-424-5p plays an important role in diseases of the reproductive system by directly targeting several key functional genes (Xu et al., 2013; Liu et al., 2018). In addition, miR-424 suppresses the proliferation and promotes the apoptosis of human ovarian granulosa cells (Du et al., 2020). These observations suggest that lncRNAs regulate goat ovarian follicular development by binding to miRNAs that regulate target gene expression.
One result we find out is really interesting. In Figure 5B, when we mutated the binding site of miR-3955-5p in XR_311113.2 sequence, the luciferase activity of the mutant pmirGLO vector significantly increased compared with wild type pmirGLO vector. The results showed that the mutant site might prevent miR-3955-5p bind to pmirGLO, so that the luciferase activity of pmirGLO will not be inhibited by miR-3955-5p. To clarify further the binding activity of miR-3955-5p, the chi-miR-3955-5p mimic was synthesized and co-transfected with the wild pmirGLO plasmid. Compared with the NC group, chi-miR-3955-5p did not significantly reduce luciferase activity of wild pmirGLO. Why? The reason may be that the mutated binding site is not necessarily the factual site of miR-3955-5p, but may be the binding site of other unknown miRNAs. When this fake site is mutated, the unknown miRNA cannot inhibit the luciferase activity. This is why we continue to use the results of Figure 5C to further eliminate the fake results of Figure 5B and ensure the reliability of the results.
MiRNAs have emerged as key post-transcriptional regulators of target expression through complementary base pairing with the target (Filipowicz et al., 2008). A large number of computational prediction tools for the prediction of putative miRNA targets have been developed, and commonly used tools in mammals include miRWalk, PicTar, miRanda, Targetscan, RNAhybrid, PITA, miRmap, DIANA-microT and RNA22V2 (Riffo-Campos et al., 2016). However, the described tools are mainly developed for a few species such as humans, mice, pigs etc. Some of these tools can only set a few parameters and cannot enter sequences. For obscure livestock such as goats, it is impossible to accurately predict the targeting relationship between miRNA and target. The lack of various databases and tools can make it difficult to carry out the study of goats.
Conclusion
In conclusion, we analyzed the lncRNA and mRNA expression profiles of goat ewe pre-ovulatory follicles to predict the interactions between lncRNAs and miRNAs. Overall, 895 lncRNAs were identified, 88 of which were preliminarily determined to show a marked differential expression, indicating potentially substantial effects on goat ovarian follicles. The tissue expression profile of lncRNA showed tissue specificity. LncRNA XR_311113.2 may function as a sponge of chi-miR-424-5p. Our RNA-seq data contributed to the known types of lncRNA species. LncRNAs may represent a promising new field of research in the area of ovarian follicular development.
Data Availability Statement
Sequencing data were deposited in the NCBI Sequence Read Archive (SRA) under accession number PRJNA648013.
Ethics Statement
The animal study was reviewed and approved by The ethics committee of Hubei Academy of Agricultural Sciences. Written informed consent was obtained from the owners for the participation of their animals in this study.
Author Contributions
Conceived and designed the experiments: HT and MC. Wrote the manuscript: HT. Performed the experiments: HT and JY. Collection of the sample: NZ and PZ. RNA extraction: XS. Critically revised the manuscript: XL and YL. All authors read and approved the final manuscript.
Funding
This study was supported by the National Natural Science Foundation of China (NSFC, 31972995 and 31501932), Hubei Province Key R&D Program (2021BBA086) and Hubei Announcement System Technology Project (2020BED028).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Acknowledgments
We thank the reviewers for their critical reading and discussion of the manuscript. The authors were grateful to Online English (https://www.oleng.com.au/) for manuscript editing.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.760416/full#supplementary-material
References
Choi, Y., and Rajkovic, A. (2006). Genetics of Early Mammalian Folliculogenesis. Cell. Mol. Life Sci. 63 (5), 579–590. doi:10.1007/s00018-005-5394-7
Cui, H. X., Zhao, S. M., Cheng, M. L., Guo, L., Ye, R. Q., Liu, W. Q., et al. (2009). Cloning and Expression Levels of Genes Relating to the Ovulation Rate of the Yunling Black Goat1. Biol. Reprod. 80 (2), 219–226. doi:10.1095/biolreprod.108.069021
Du, J., Lin, X., Wu, R., Gao, Z., Du, Y., Liao, Y., et al. (2020). miR-424 Suppresses Proliferation and Promotes Apoptosis of Human Ovarian Granulosa Cells by Targeting Apelin and APJ Expression. Am. J. Transl. Res. 12 (7), 3660–3673.
Du, Z., Sun, T., Hacisuleyman, E., Fei, T., Wang, X., Brown, M., et al. (2016). Integrative Analyses Reveal a Long Noncoding RNA-Mediated Sponge Regulatory Network in Prostate Cancer. Nat. Commun. 7, 10982. doi:10.1038/ncomms10982
Enright, A. J., John, B., Gaul, U., Tuschl, T., Sander, C., and Marks, D. S. (2003). MicroRNA Targets in Drosophila. Genome Biol. 5 (1), R1. doi:10.1186/gb-2003-5-1-r1
Field, S. L., Dasgupta, T., Cummings, M., and Orsi, N. M. (2014). Cytokines in Ovarian Folliculogenesis, Oocyte Maturation and Luteinisation. Mol. Reprod. Dev. 81 (4), 284–314. doi:10.1002/mrd.22285
Filipowicz, W., Bhattacharyya, S. N., and Sonenberg, N. (2008). Mechanisms of post-transcriptional Regulation by microRNAs: Are the Answers in Sight? Nat. Rev. Genet. 9 (2), 102–114. doi:10.1038/nrg2290
Han, P., and Chang, C.-P. (2015). Long Non-coding RNA and Chromatin Remodeling. RNA Biol. 12 (10), 1094–1098. doi:10.1080/15476286.2015.1063770
Huang, Y. (2018). The Novel Regulatory Role of lncRNA-miRNA-mRNA axis in Cardiovascular Diseases. J. Cel. Mol. Med. 22 (12), 5768–5775. doi:10.1111/jcmm.13866
Huang, Y., Zhang, J., Hou, L., Wang, G., Liu, H., Zhang, R., et al. (2017). LncRNA AK023391 Promotes Tumorigenesis and Invasion of Gastric Cancer through Activation of the PI3K/Akt Signaling Pathway. J. Exp. Clin. Cancer Res. 36 (1), 194. doi:10.1186/s13046-017-0666-2
Hung, T., and Chang, H. Y. (2010). Long Noncoding RNA in Genome Regulation. RNA Biol. 7 (5), 582–585. doi:10.4161/rna.7.5.13216
Ji, P., Diederichs, S., Wang, W., Böing, S., Metzger, R., Schneider, P. M., et al. (2003). MALAT-1, a Novel Noncoding RNA, and Thymosin β4 Predict Metastasis and Survival in Early-Stage Non-small Cell Lung Cancer. Oncogene 22 (39), 8031–8041. doi:10.1038/sj.onc.1206928
Jin, J. J., Lv, W., Xia, P., Xu, Z. Y., Zheng, A. D., Wang, X. J., et al. (2018). Long Noncoding RNASYISLregulates Myogenesis by Interacting with Polycomb Repressive Complex 2. Proc. Natl. Acad. Sci. USA 115 (42), E9802–E9811. doi:10.1073/pnas.1801471115
Johnson, A. L. (2015). Ovarian Follicle Selection and Granulosa Cell Differentiation. Poult. Sci. 94 (4), 781–785. doi:10.3382/ps/peu008
Kim, B. H., Ju, W. S., Kim, J.-S., Kim, S.-U., Park, S. J., Ward, S. M., et al. (2020). Effects of Gangliosides on Spermatozoa, Oocytes, and Preimplantation Embryos. Int. J. Mol. Sci. 21 (1), 106. doi:10.3390/ijms21010106
Kong, L., Zhang, Y., Ye, Z.-Q., Liu, X.-Q., Zhao, S.-Q., Wei, L., et al. (2007). CPC: Assess the Protein-Coding Potential of Transcripts Using Sequence Features and Support Vector Machine. Nucleic Acids Res. 35, W345–W349. doi:10.1093/nar/gkm391
Li, H., Yang, J., Jiang, R., Wei, X., Song, C., Huang, Y., et al. (2018). Long Non-coding RNA Profiling Reveals an Abundant MDNCR that Promotes Differentiation of Myoblasts by Sponging miR-133a. Mol. Ther. - Nucleic Acids 12, 610–625. doi:10.1016/j.omtn.2018.07.003
Li, R., Harvey, A. R., Hodgetts, S. I., and Fox, A. H. (2017). Functional Dissection of NEAT1 Using Genome Editing Reveals Substantial Localization of the NEAT1_1 Isoform outside Paraspeckles. RNA 23 (6), 872–881. doi:10.1261/rna.059477.116
Li, Y., Liu, Y.-d., Chen, S.-l., Chen, X., Ye, D.-s., Zhou, X.-y., et al. (2018a). Down-regulation of Long Non-coding RNAMALAT1inhibits Granulosa Cell Proliferation in Endometriosis by Up-Regulating P21 via Activation of the ERK/MAPK Pathway. Mol. Hum. Reprod. 25 (1), 17–29. doi:10.1093/molehr/gay045
Li, Y., Wang, H., Zhou, D., Shuang, T., Zhao, H., and Chen, B. (2018b). Up-regulation of Long Noncoding RNA SRA Promotes Cell Growth, Inhibits Cell Apoptosis, and Induces Secretion of Estradiol and Progesterone in Ovarian Granular Cells of Mice. Med. Sci. Monit. 24, 2384–2390. doi:10.12659/MSM.907138
Liu, J., Gu, Z., Tang, Y., Hao, J., Zhang, C., and Yang, X. (2018). Tumour-suppressive microRNA-424-5p Directly Targets CCNE1 as Potential Prognostic Markers in Epithelial Ovarian Cancer. Cell Cycle 17 (3), 309–318. doi:10.1080/15384101.2017.1407894
Liu, Y.-d., Li, Y., Feng, S.-x., Ye, D.-s., Chen, X., Zhou, X.-y., et al. (2017). Long Noncoding RNAs: Potential Regulators Involved in the Pathogenesis of Polycystic Ovary Syndrome. Endocrinology 158 (11), 3890–3899. doi:10.1210/en.2017-00605
Lu, W., Huang, S. Y., Su, L., Zhao, B. X., and Miao, J. Y. (2016). Long Noncoding RNA LOC100129973 Suppresses Apoptosis by Targeting MIR-4707-5p and MIR-4767 in Vascular Endothelial Cells. Sci. Rep. 6 (1), 21620. doi:10.1038/srep21620
Malan, S. W. (2000). The Improved Boer Goat. Small Ruminant Res. 36 (2), 165–170. doi:10.1016/S0921-4488(99)00160-1
McHugh, C. A., Chen, C.-K., Chow, A., Surka, C. F., Tran, C., McDonel, P., et al. (2015). The Xist lncRNA Interacts Directly with SHARP to Silence Transcription through HDAC3. Nature 521 (7551), 232–236. doi:10.1038/nature14443
McLaughlin, E. A., and McIver, S. C. (2009). Awakening the Oocyte: Controlling Primordial Follicle Development. Reproduction 137 (1), 1–11. doi:10.1530/REP-08-0118
Mercer, T. R., Dinger, M. E., and Mattick, J. S. (2009). Long Non-coding RNAs: Insights into Functions. Nat. Rev. Genet. 10 (3), 155–159. doi:10.1038/nrg2521
Naing, S. W., Wahid, H., Mohd Azam, K., Rosnina, Y., Zuki, A. B., Kazhal, S., et al. (2010). Effect of Sugars on Characteristics of Boer Goat Semen after Cryopreservation. Anim. Reprod. Sci. 122 (1-2), 23–28. doi:10.1016/j.anireprosci.2010.06.006
Nakagawa, S., Shimada, M., Yanaka, K., Mito, M., Arai, T., Takahashi, E., et al. (2014). The lncRNA Neat1 Is Required for Corpus Luteum Formation and the Establishment of Pregnancy in a Subpopulation of Mice. Dev 141 (23), 4618–4627. doi:10.1242/dev.110544
Necsulea, A., Soumillon, M., Warnefors, M., Liechti, A., Daish, T., Zeller, U., et al. (2014). The Evolution of lncRNA Repertoires and Expression Patterns in Tetrapods. Nature 505 (7485), 635–640. doi:10.1038/nature12943
Nguyen, X. P., Nakamura, T., Osuka, S., Bayasula, B., Nakanishi, N., Kasahara, Y., et al. (2019). Effect of the Neuropeptide Phoenixin and its Receptor GPR173 during Folliculogenesis. Reproduction 158 (1), 25–34. doi:10.1530/REP-19-0025
Pan, Z., Zhang, J., Li, Q., Li, Y., Shi, F., Xie, Z., et al. (2012). Current Advances in Epigenetic Modification and Alteration during Mammalian Ovarian Folliculogenesis. J. Genet. Genomics 39 (3), 111–123. doi:10.1016/j.jgg.2012.02.004
Paraskevopoulou, M. D., and Hatzigeorgiou, A. G. (2016). Analyzing MiRNA-LncRNA Interactions. Methods Mol. Biol. 1402 (1), 271–286. doi:10.1007/978-1-4939-3378-5_21
Parikshak, N. N., Swarup, V., Belgard, T. G., Irimia, M., Ramaswami, G., Gandal, M. J., et al. (2016). Genome-wide Changes in lncRNA, Splicing, and Regional Gene Expression Patterns in Autism. Nature 540 (7633), 423–427. doi:10.1038/nature20612
Peng, W.-X., Koirala, P., and Mo, Y.-Y. (2017). LncRNA-mediated Regulation of Cell Signaling in Cancer. Oncogene 36 (41), 5661–5667. doi:10.1038/onc.2017.184
Punta, M., Coggill, P. C., Eberhardt, R. Y., Mistry, J., Tate, J., Boursnell, C., et al. (2012). The Pfam Protein Families Database. Nucleic Acids Res. 40 (1), D290–D301. doi:10.1093/nar/gkr1065
Ren, H., Wang, G., Chen, L., Jiang, J., Liu, L., Li, N., et al. (2016). Genome-wide Analysis of Long Non-coding RNAs at Early Stage of Skin Pigmentation in Goats (Capra hircus). BMC Genomics 17, 67. doi:10.1186/s12864-016-2365-3
Riffo-Campos, Á., Riquelme, I., and Brebi-Mieville, P. (2016). Tools for Sequence-Based miRNA Target Prediction: What to Choose? Int. J. Mol. Sci. 17 (12), 1987. doi:10.3390/ijms17121987
Saharrea, A., Valencia, J., Balcázar, A., Mejía, O., Cerbón, J. L., Caballero, V., et al. (1998). Premature Luteal Regression in Goats Superovulated with PMSG: Effect of hCG or GnRH Administration during the Early Luteal Phase. Theriogenology 50 (7), 1039–1052. doi:10.1016/S0093-691X(98)00206-4
Sarropoulos, I., Marin, R., Cardoso-Moreira, M., and Kaessmann, H. (2019). Developmental Dynamics of lncRNAs across Mammalian Organs and Species. Nature 571 (7766), 510–514. doi:10.1038/s41586-019-1341-x
Schober, P., Boer, C., and Schwarte, L. A. (2018). Correlation Coefficients. Anesth. Analgesia 126 (5), 1763–1768. doi:10.1213/ANE.0000000000002864
Simion, V., Haemmig, S., and Feinberg, M. W. (2019). LncRNAs in Vascular Biology and Disease. Vasc. Pharmacol. 114, 145–156. doi:10.1016/j.vph.2018.01.003
Sugiyama, M., Sumiya, M., Shirasuna, K., Kuwayama, T., and Iwata, H. (2016). Addition of Granulosa Cell Mass to the Culture Medium of Oocytes Derived from Early Antral Follicles Increases Oocyte Growth, ATP Content, and Acetylation of H4K12. Zygote 24 (6), 848–856. doi:10.1017/S0967199416000198
Sun, L., Luo, H., Bu, D., Zhao, G., Yu, K., Zhang, C., et al. (2013). Utilizing Sequence Intrinsic Composition to Classify Protein-Coding and Long Non-coding Transcripts. Nucleic Acids Res. 41 (17), e166. doi:10.1093/nar/gkt646
Tao, H., Xiong, Q., Zhang, F., Zhang, N., Liu, Y., Suo, X., et al. (2018). Circular RNA Profiling Reveals Chi_circ_0008219 Function as microRNA Sponges in Pre-ovulatory Ovarian Follicles of Goats (Capra hircus). Genomics 110 (4), 257–266. doi:10.1016/j.ygeno.2017.10.005
Torner, H., Brüssow, K.-P., Alm, H., Ràtky, J., and Kanitz, W. (1998). Morphology of Porcine Cumulus-Oocyte-Complexes Depends on the Stage of Preovulatory Maturation. Theriogenology 50 (1), 39–48. doi:10.1016/S0093-691X(98)00111-3
Trapnell, C., Roberts, A., Goff, L., Pertea, G., Kim, D., Kelley, D. R., et al. (2012). Differential Gene and Transcript Expression Analysis of RNA-Seq Experiments with TopHat and Cufflinks. Nat. Protoc. 7 (3), 562–578. doi:10.1038/nprot.2012.016
Uyar, A., Torrealday, S., and Seli, E. (2013). Cumulus and Granulosa Cell Markers of Oocyte and Embryo Quality. Fertil. Sterility 99 (4), 979–997. doi:10.1016/j.fertnstert.2013.01.129
Wang, J.-J., Zhang, T., Chen, Q.-M., Zhang, R.-Q., Li, L., Cheng, S.-F., et al. (2020). Genomic Signatures of Selection Associated with Litter Size Trait in Jining Gray Goat. Front. Genet. 11, 286. doi:10.3389/fgene.2020.00286
Wang, J., Gong, X., Tian, G. G., Hou, C., Zhu, X., Pei, X., et al. (2018). Long Noncoding RNA Growth Arrest-specific 5 Promotes Proliferation and Survival of Female Germline Stem Cells In Vitro. Gene 653 (4), 14–21. doi:10.1016/j.gene.2018.02.021
Wang, K., Jin, W., Song, Y., and Fei, X. (2017). LncRNA RP11-436H11.5, Functioning as a Competitive Endogenous RNA, Upregulates BCL-W Expression by Sponging miR-335-5p and Promotes Proliferation and Invasion in Renal Cell Carcinoma. Mol. Cancer 16 (1), 166. doi:10.1186/s12943-017-0735-3
Wang, Y., Fan, X., Lin, F., He, G., Terzaghi, W., Zhu, D., et al. (2014). Arabidopsis Noncoding RNA Mediates Control of Photomorphogenesis by Red Light. Proc. Natl. Acad. Sci. 111 (28), 10359–10364. doi:10.1073/pnas.1409457111
Xiao, J., Lin, L., Luo, D., Shi, L., Chen, W., Fan, H., et al. (2020). Long Noncoding RNA TRPM2-AS Acts as a microRNA Sponge of miR-612 to Promote Gastric Cancer Progression and Radioresistance. Oncogenesis 9 (3), 29. doi:10.1038/s41389-020-0215-2
Xu, J., Li, Y., Wang, F., Wang, X., Cheng, B., Ye, F., et al. (2013). Suppressed miR-424 Expression via Upregulation of Target Gene Chk1 Contributes to the Progression of Cervical Cancer. Oncogene 32 (8), 976–987. doi:10.1038/onc.2012.121
Yan, L., Zhou, J., Gao, Y., Ghazal, S., Lu, L., Bellone, S., et al. (2015). Regulation of Tumor Cell Migration and Invasion by the H19/let-7 axis Is Antagonized by Metformin-Induced DNA Methylation. Oncogene 34 (23), 3076–3084. doi:10.1038/onc.2014.236
Yan, S., Wang, P., Wang, J., Yang, J., Lu, H., Jin, C., et al. (2019). Long Non-coding RNA HIX003209 Promotes Inflammation by Sponging miR-6089 via TLR4/NF-Κb Signaling Pathway in Rheumatoid Arthritis. Front. Immunol. 10, 2218. doi:10.3389/fimmu.2019.02218
Zhan, S., Dong, Y., Zhao, W., Guo, J., Zhong, T., Wang, L., et al. (2016). Genome-wide Identification and Characterization of Long Non-coding RNAs in Developmental Skeletal Muscle of Fetal Goat. BMC Genomics 17 (1), 666. doi:10.1186/s12864-016-3009-3
Zhao, X., Li, J., Lian, B., Gu, H., Li, Y., and Qi, Y. (2018). Global Identification of Arabidopsis lncRNAs Reveals the Regulation of MAF4 by a Natural Antisense RNA. Nat. Commun. 9 (1), 5056. doi:10.1038/s41467-018-07500-7
Keywords: goat, ovarian follicle, lncRNA, transcriptome, miRNA sponge
Citation: Tao H, Yang J, Zhang P, Zhang N, Suo X, Li X, Liu Y and Chen M (2022) Characterization of XR_311113.2 as a MicroRNA Sponge for Pre-ovulatory Ovarian Follicles of Goats via Long Noncoding RNA Profile and Bioinformatics Analysis. Front. Genet. 12:760416. doi: 10.3389/fgene.2021.760416
Received: 18 August 2021; Accepted: 10 December 2021;
Published: 03 January 2022.
Edited by:
Marcos De Donato, Instituto de Tecnología y Educación Superior de Monterrey (ITESM), MexicoReviewed by:
Xin Qi, Ocean University of China, ChinaShahin Eghbalsaied, Islamic Azad University, Isfahan, Iran
Copyright © 2022 Tao, Yang, Zhang, Zhang, Suo, Li, Liu and Chen. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Mingxin Chen, cmx1963@163.com