Evolution of the ErbB gene family and analysis of regulators of Egfr expression during development of the rat spinal cord.

Egfr, a member of the ErbB gene family, plays a critical role in tissue development and homeostasis, wound healing, and disease. However, expression and regulators of Egfr during spinal cord development remain poorly understood. In this study, we investigated ErbB evolution and analyzed co-expression modules, miRNAs, and transcription factors that may regulate Egfr expression in rats. We found that ErbB family members formed via Egfr duplication in the ancient vertebrates but diverged after speciation of gnathostomes. We identified a module that was co-expressed with Egfr, which involved cell proliferation and blood vessel development. We predicted 25 miRNAs and nine transcription factors that may regulate Egfr expression. Dual-luciferase reporter assays showed six out of nine transcription factors significantly affected Egfr promoter reporter activity. Two of these transcription factors (KLF1 and STAT3) inhibited the Egfr promoter reporter, whereas four transcription factors (including FOXA2) activated the Egfr promoter reporter. Real-time PCR and immunofluorescence experiments showed high expression of FOXA2 during the embryonic period and FOXA2 was expressed in the floor plate of the spinal cord, suggesting the importance of FOXA2 during embryonic spinal cord development. Considering the importance of Egfr in embryonic spinal cord development, wound healing, and disease (specifically in cancer), regulatory elements identified in this study may provide candidate targets for nerve regeneration and disease treatment in the future.

Introduction

Epidermal growth factor receptor (EGFR, also known as ErbB1) is a member of the receptor tyrosine kinase (RTK) superfamily. The other three EGFR family members include ErbB2, ErbB3, and ErbB4 (Ayati et al., 2020). EGFR was the first receptor discovered in this family. Until now, at least seven primary ligands can activate EGFR, including epidermal growth factor (EGF), transforming growth factor-α, heparin-binding EGF, β-cellulin, amphiregulin, epiregulin, and epigen (Knudsen et al., 2014). The ErbB family plays a critical role in vital signaling pathways, such as the phosphatidylinositol 3-kinase/Akt, Ras/Raf/MEK/ERK1/2, and phospholipase C pathways, which are essential for cellular proliferation, differentiation, and survival (Schlessinger, 2002). Each member of the ErbB family contains an EGF-like domain (a conserved three-loop structure). EGFR consists of an extracellular ligand-binding domain, single transmembrane domain, and cytoplasmic domain. After ligand binding, EGFR can homo- or hetero-dimerize, and the EGFR-ErbB2 dimer forms the most active receptor (Klapper et al., 2000). After dimerization, specific tyrosine residues are phosphorylated and new binding sites are formed for adaptor proteins with Src homology 2 domains, which enable the activation of downstream signaling pathways.

EGFR plays a crucial role during vertebrate development. In mammals, EGFR improves embryo implantation and placenta development during embryogenesis (Chia et al., 1995). EGFR is involved in the development of many organs, including the brain, spinal cord, heart, bone, and liver (Chen et al., 2016). EGFR can be detected in cells of the central nervous system (CNS) during development, such as neurons, astrocytes, oligodendrocytes, and progenitor cells of the subventricular zone (Seroogy et al., 1994; Kornblum et al., 1997). However, the expression of EGFR is reduced in the adult brain, including in the nigrostriatal system, and the highest levels of expression have been detected during the neonatal stage (Seroogy et al., 1994). In the peripheral nervous system, EGFR is also found in the neurons of dorsal root ganglia (DRG), Schwann cells, and satellite glial cells (Morris et al., 1999). Sibilia et al. (2007) reported that EGFR-null mice displayed serious neural defects. Hence, EGFR participates in the development, differentiation, maintenance, and regeneration of various organs and the nervous system. Understanding the regulation of Egfr expression is of great significance for comprehending the development of the nervous system and exploring potential targeted therapies for related diseases and such spinal cord injury.

Gene expression is regulated by several factors, including chromatin modification, transcription factors (TFs), and non-coding RNAs (Kouzarides, 2007; Patil et al., 2014). In addition, accumulating studies showed that activation of G protein-coupled receptor (GPCR) by agonists (e.g. angiotensin II) could promote EGFR activity via transactivation (Forrester et al., 2016). TFs recognize cis-regulatory regions to correctly initiate gene transcription. MicroRNAs (miRNA, 20–30 NT) regulate gene expression post-transcriptionally by interacting with the 3’ untranslated region (3’-UTR) of target mRNAs (Filipowicz et al., 2008). TFs and miRNAs are important regulators of gene expression in eukaryotes. Several TFs (Sp1, EVI1, ETF) have been shown to activate Egfr (Kageyama et al., 1988; Mizuguchi et al., 2019; Tsai et al., 2019). Chromatin immunoprecipitation (ChIP)-seq technology has accelerated the discovery of TF binding sites (Mundade et al., 2014), and the majority have been studied in humans and mice. MiRNAs, such as miR-34a (Li et al., 2017) and miR-133a (Guo et al., 2018), and long non-coding RNAs, such as LINC00265 (Ge et al., 2019) and EGFR-AS1 (Wang et al., 2019a), have been validated as regulators of Egfr expression in cancers. In our previously sequenced transcription profile of rat spinal cord development (Yang et al., 2021), we demonstrated that Egfr was highly expressed in the late embryonic period (embryonic day 18, E18d) until postnatal week 1 (P1w) with low expression in adults (postnatal weeks 4–8), suggesting dynamic regulation of Egfr during spinal cord development. However, limited studies have explored potential regulators of Egfr expression in the spinal cord.

In this study, we performed an evolutionary analysis of the ErbB gene family in metazoans, screened a module containing genes that were co-expressed with Egfr, investigated GPCRs that may be involved in EGFR transactivation, and predicted miRNAs that targeted Egfr and TFs that regulated Egfr expression based on our previous transcription profiling of spinal cord development. We also validated the direct binding between TFs and Egfr using a dual luciferase reporter assay in vitro. Finally, we investigated the expression of candidate TFs during spinal cord development at the mRNA and protein level using qRT-PCR and immunohistochemistry, respectively, in rats ranging in age from embryonic to adult.

Materials and Methods

Animals and tissue preparation

The study was approved by the Animal Care and Use Committee of Nantong University and was reported in accordance with the ARRIVE 2.0 guidelines (Animal Research: Reporting of In Vivo Experiments) (Percie du Sert et al., 2020). Specific pathogen-free Sprague-Dawley rats, including pregnant rats (body weights: 300–330 g; n = 2), male pups (postnatal day 1 [P1d] and postnatal week 1 [P1w]; body weights: 7 g and 20 g, respectively; n = 4/group), and adult males (postnatal week 8 [P8w] and 12 [P12w]; average body weights: 200 g and 320 g, respectively; n = 4/group) were used for sample collection at different developmental stages. All rats were purchased from the Experimental Animal Center of Nantong University, Nantong, China (license No. SCXK (Su) 2019-0001). Four rats were housed per cage, and all rats were allowed access to food and water ad libitum. The animal room temperature was maintained at 23 ± 1°C. Spinal cord tissue samples were collected from embryos of pregnant rats at embryo day 11 (E11d), 13 (E13d), and 14 (E14d) and from pups/adults at P1d, P1w, and P8w. The sampling time was 6:30 a.m. The pregnant rats were deeply anesthetized by intraperitoneal injection of mixed narcotics (4.25 g chloral hydrate, 2.12 g magnesium sulfate, 886 mg sodium pentobarbital, 14.25 mL ethanol, 33.8 mL propylene glycol, and distilled water to a final volume of 100 mL) and the embryos were removed. The embryos were placed in a pre-prepared, oxygen-saturated artificial cerebrospinal fluid (ACSF, ice-water mixture; ACSF composition: 130 mM NaCl, 5 mM KCl, 2 mM KH2PO4, 1.5 mM CaCl2, 6 mM MgSO4, 10 mM glucose, 10 mM HEPES, pH 7.2, 305 mOsm osmotic pressure). The neural tubes of E11d rat embryos displayed a translucent, soft gel state. The strip outline of the neural tube and its dark boundary with surrounding tissues were observed under a microscope (ZEISS Type: Stemi 508, Jena, Germany). The neural tubes were carefully isolated using micro scissors and micro tweezers. For E13d, 14d, and E18d rat embryos, the cartilaginous spinal canals were carefully cut open with micro scissors, and the spinal cords were isolated using micro tweezers. Adult rats (P8w) were anesthetized by intraperitoneal injection of mixed narcotics, placed on an ice bag, and decapitated at the neck with a specialized decapitator. The spinal columns were separated immediately and washed twice with pre-prepared, oxygen-saturated ACSF. The spines were dissected in ACSF (ice-water mixture), and the muscles, bones, dura mater, arachnoid membranes, pia mater, and large blood vessels were carefully removed. All spinal cords were immediately quick-frozen and preserved in liquid nitrogen (–180°C).

Evolutionary analysis of the ErbB gene family in 32 metazoan species

We used BLAST (v2.2.26, https://ftp.ncbi.nlm.nih.gov/blast/) (Altschul et al., 1990) to identify potential ErbB gene family members from the proteomes of 32 species. The longest protein sequence was used if multiple isoforms existed. The phylogenetic relationship of 32 species was inferred by combining published literature. To understand the evolution of ErbB genes in the 32 species, we performed phylogenetic analysis as follows. Multiple sequence alignments were performed using MUSCLE (v3.8.31, https://www.ebi.ac.uk/Tools/msa/muscle/) (Edgar, 2004) and poor alignment was further trimmed using TrimAl (v1.2, http://trimal.cgenomics.org) (Capella-Gutierrez et al., 2009) with the parameter “-automated1”. Finally, IQ-TREE (v1.6.7, http://www.iqtree.org) (Nguyen et al., 2015) was used to reconstruct the phylogenetic tree using the maximum-likelihood algorithm. The iTOL tool (https://itol.embl.de/) (Letunic and Bork, 2019) was used for gene tree visualization. Synteny analysis was performed using JCVI (https://github.com/tanghaibao/jcvi) (Tang et al., 2015). MEGA-X (https://www.megasoftware.net/) (Kumar et al., 2018) was employed to analyze let-7 family evolution.

Prediction of regulatory miRNAs targeting Egfr

Four transcripts of the Egfr gene were present in rat reference annotation (RefSeq accession: GCF_000001895.5, rn6). The sequence of the longest 3’-UTR was extracted from the rat Egfr gene (accession No. XM_008770416.2). Mature miRNA sequences for rats were retrieved from the miRbase database (http://www.mirbase.org/) (Kozomara et al., 2019). The miRanda software (v3.3a, https://cbio.mskcc.org/miRNA2003/miranda.html) (Enright et al., 2003) was used to predict miRNAs targeting Egfr with the criterion of “-sc 140 -en -20”. To avoid multiple predictions, we also employed two other widely used programs: PITA (http://genie.weizmann.ac.il/pubs/mir07) (Kertesz et al., 2007) and TargetScan (http://www.targetscan.org/vert_72/) (Agarwal et al., 2015) using default parameters. MiRNAs with experimentally validated interactions in Tarbase (v8.0, http://www.microrna.gr/tarbase) (Karagkouni et al., 2018) were also collected. MiRNAs predicted by at least two programs were chosen for downstream analyses. Pearson’s correlation of expression between Egfr and miRNAs in rats was calculated using the function “cor.test” in the R program (https://www.r-project.org/) (Ripley, 2001) based on our previous RNA-seq data of rat spinal cord development. Significant negative correlations were used to further select candidate regulatory miRNAs (P < 0.05).

Prediction of potential TFs that regulate the expression of Egfr

We extracted a 5 kb sequence upstream of rat Egfr transcriptional start site and predicted potential promoter regions using the webserver Promoter 2.0 (http://www.cbs.dtu.dk/services/Promoter/) (Knudsen, 1999). Then, we predicted and screened potential regulatory TFs and binding sites by combining four methods: 1) JASPAR database (https://jaspar.genereg.net/) (Mathelier et al., 2016); 2) TRANSFAC database (https://genexplain.com/transfac/) (Matys et al., 2003) (v8.3); 3) animal TFDB database (http://bioinfo.life.hust.edu.cn/AnimalTFDB) (Zhang et al., 2015); and 4) MAST algorithms (https://meme-suite.org/meme/, from The MEME Suite) (Bailey and Gribskov, 1998). We focused on predicted TFs with Gene Ontology annotation by assigning “angiogenesis” and high expression in the early embryonal stage. ChIP-seq evidence of TFs targeting Egfr in humans and mice was also evaluated using ChIPBase (v2.0, https://rna.sysu.edu.cn/chipbase/) (Zhou et al., 2017).

Validation of direct interaction between TFs and the Egfr promoter in rats by dual-luciferase reporter assay

To further confirm nine predicted TFs, we employed a dual-luciferase reporter assay (Sherf et al., 1996). Sequences of nine TFs and Egfr promoter region were retrieved from the Genbank database and synthesized by GeneChem (Shanghai, China). GV141-TF (Foxa2, Sox17, Bhlha15, Klf4, Sp1, Ddit3, Arid3b, Klf1, and Stat3) over-expression vectors were constructed by GeneChem. A GV141-empty (control) vector was also provided by GeneChem. The plasmids individually containing full-length coding sequences of the nine TFs were used as the effectors. The plasmid containing the Egfr promoter was used as the reporter. Human embryonic kidney 293T cells (Cell bank of Shanghai Chinese Academy of Sciences, Shanghai, China, RRID:19375.09.3101HUMSCSP502) were transfected with these vectors (an effector and the reporter) using X-tremeGENE HP transfection reagent (Roche, Basel, Switzerland) in accordance with the manufacturer’s instructions. After 48 hours, the cells were harvested, and a dual-luciferase assay kit (Promega Biotech Co., Ltd.; Beijing, China) was used to detect LUC and REN luciferase activity via a Luminoskan Ascent Microplate Luminometer (Thermo Scientific, Rockford, IL, USA) in accordance with the manufacturer’s instructions. The transcriptional activity of each TF after binding to the rat Egfr promoter was estimated using the LUC/REN value. At least three independent experiments were performed for each TF-promoter pair.

Quantitative reverse transcription-polymerase chain reaction (qRT-PCR) validation

The mRNA expression of Egfr and six candidate TFs (Bhlha15, Ddit3, Foxa2, Klf1, Sox17, and Stat3) in the developing rat spinal cord was examined by quantitative reverse transcription-polymerase chain reaction (qRT-PCR). Total RNA was separately extracted from the spinal cord of E13d, E18d, P1w, and P8w rats using TRIzol (Gibco; Carlsbad, CA, USA) in accordance with the manufacturer’s protocol. Egfr expression was analyzed at the additional time points E11d, E14d, P1d, and P12w. Total RNA (1 μg) was used to synthesize cDNA using the PrimeScriptTM RT reagent kit (TaKaRa Biotechnology CO., Ltd.; Dalian, Liaoning Province, China). The qRT-PCR was performed using the TB Green® Premix Ex TaqTM (TaKaRa). The expression of target genes was normalized to Gapdh expression using the ΔCt method (2–ΔΔCt). The primers used in this study are listed in .

Immunofluorescence staining

Immunofluorescence staining was performed on tissues from E11d, E14d, E18d, P1d, P1w, and P8w rats to determine the expression and cellular localization EGFR. Immunofluorescence staining was performed on tissues from E13d, E18d, P1w, and P8w rats to evaluate FOXA2 expression and cellular localization during spinal cord development. For embryos (E11d–E18d) and neonatal rats (P1d, P1w), the spinal cord tissue was dissected in ice-cold phosphate-buffered saline (PBS) after a 75% ethanol wipe. Animal anesthesia method was the same with description above. To obtain the adult spinal cord (P8w), the animals were transcardially perfused with normal saline followed by 4% paraformaldehyde under deep anesthesia via intraperitoneal injections of mixed narcotics as described above. All collected spinal cord tissues were post-fixed in 4% paraformaldehyde at 4°C for 8 hours, dehydrated in sequential 10–20–30% sucrose solutions, and sectioned on a Leica cryostat (Leica, Wetzlar, Germany) into 12-μm-thick sections. Immunofluorescence procedures were performed as follows. The sections were placed in blocking buffer (Beyotime; Shanghai, China), incubated with the primary antibody at 4°C for 12 hours, and then treated with the secondary antibody at room temperature for one hour. Subsequently, the nuclei were counterstained with Hoechst 33342 in PBS (1:2000) for 10 minutes at room temperature. Images were observed using fluorescence microscopy (Leica). To quantify protein expression, the “integrated density” of the target channels (EGFR and FOXA2 only) with the same area at each developmental stage was calculated using ImageJ Fiji software (https://imagej.net/software/fiji/) (Schindelin et al., 2012) and then normalized to the first developmental stage (E11d for EGFR and E13d for FOXA2). Nestin was used as a marker of neural stem cells in the embryonic rat spinal cord, and NeuN was used as the marker of neuronal nuclei in the postnatal rat spinal cord. Primary and secondary antibodies used in this study are listed as follows: rabbit anti-EGFR (1:100, Proteintech, Rosemont, IL, USA, Cat# 18986-1-AP, RRID: AB_10596476), rabbit anti-FOXA2 (1:300, Abcam, Cambridge, UK, Cat# ab108422, RRID: AB_11157157), mouse anti-Nestin (1:100, Millipore, Bedford, MA, USA, Cat# MAB353, RRID: AB_94911), mouse anti-NeuN (1:100, MilliporeSigma, Burlington, MA, USA, Cat# MAB377A5, RRID: AB_2814948), Alexa Fluor 488 goat anti-mouse IgG (H+L) (1:2000, Thermo Fisher Scientific, Waltham, MA, USA, Cat# A-11001, RRID: AB_2534069), and Cy3-conjugated goat anti-rabbit IgG (H+L) (1:400, Proteintech, Cat# SA00009-2, RRID: AB_2890957).

Statistical analysis

No statistical methods were used to predetermine sample sizes; however, our sample sizes were similar to those reported in a previous publication (Zha et al., 2016). Histological analyses were blinded to experimental conditions. Data from the dual-luciferase reporter assays are presented as mean ± SD and were analyzed using one-way analysis of variance (ANOVA) followed by the Dunnett multiple comparison test. Quantitative data from qRT-PCR and immunofluorescence experiments are presented as means ± SEM after analysis using Student’s t-tests. All statistical analyses were performed using GraphPad Prism (v8.3, GraphPad Software Inc., San Diego, CA, USA). P-values < 0.05 were considered statistically significant.

Results

ErbB gene family members are duplicated in vertebrates but diverge in Gnathostomata

We analyzed ErbB family genes in 32 species from Eukaryota (Placozoa, Cnidaria, Protostomia, and Deuterostomia) to understand the evolution of this gene family. The results indicated that ErbB genes occurred in species from Metazoa except for the phyla Placozoa and Cnidaria (). Multiple copies were observed in most species from vertebrate lineages but only in some invertebrate species, such as those in Platyhelminthes and in the earthworm Eisenia andrei (). According to the annotation in the public database, the ErbB gene in invertebrates was Egfr, while Egfr, Erbb2, Erbb3, and Erbb4 were found in vertebrates except for species from Agnatha which showed the best hit with Erbb4 (termed Erbb4-like). Phylogenetic analysis of ErbB genes indicated that Egfr was duplicated in the last common ancestor of Schistosoma mansoni and Schmidtea mediterranea and subsequently underwent independent expansion (). Egfr in E. andrei underwent species-specific expansion with a total of 11 copies. Interestingly, although duplication of ErbB genes was observed in vertebrates, divergence of ErbB genes occurred in Gnathostomata. Five to six ErbB genes in three species belonging to Agnatha were not diverged, but four diverged genes were present in the closest species, i.e., cartilaginous fish (shark), and remaining gnathostomes, such as humans and mice (). ErbB genes in species from Agnatha were clustered together with Erbb3-Erbb4 and divided into two subgroups before speciation (). Six ErbB genes in lampreys were located in six different scaffolds. Synteny analysis between humans and lampreys showed that blocks containing four ErbB genes in humans had multiple orthologs and cross relationships in lampreys (). These data suggested that the ErbB gene family may have originated via small-scale duplication of ancient Egfr in vertebrates. Gene duplication is an important source of novelty; however, most duplicates are eventually lost, and the remains may be preserved by neofunctionalization and subfunctionalization (Gout and Lynch, 2015). Transcriptomic information provides more details about the potential fate of duplicated genes. Planarian flatworms (S. mediterranea) have powerful abilities of regeneration after amputation owing to abundant stem cells (neoblasts) throughout the body (Adler et al., 2014). RNA interference (RNAi) knockdown that inhibited the expression of Egfr homologues smed-egfr-1–6 in S. mediterranea showed expression divergence in tissues with different copies even in the same group and one copy in each group mainly contributed to the phenotype (Barberan et al., 2016). It was also shown that RNAi specific for smed-egfr-1 and smed-egfr-3 influenced regeneration (Barberan et al., 2016). Nine out of 11 copies of Egfr in E. andrei were consistently up-regulated at 6 hours after injury, suggesting that most expanded Egfr copies were possibly a gene dosage response to injury and regeneration (Shao et al., 2020) (). In vertebrates, time-series transcriptome analysis of sea lamprey embryogenesis showed that most ErbB genes were upregulated in the early (1–2.5 days) period and decreased later (3–5 days) post-fertilization (Bryant et al., 2016) (). Expression divergence of different ErbB genes was observed in our previous rat spinal cord development dataset (Yang et al., 2021) ().

Identification and evolution of the ErbB gene family in Metazoa.

(A) Number of genes from the ErbB gene family that were found in 32 metazoan species. (B) The phylogeny of ErbB gene family members in 32 species. Four ellipses were used to label four typical members of the ErbB family in Gnathostomata. Egfr from A. queenslandica was used as a root. Leaf shapes and colors are the same as in (A) that depict different species. (C) Syntenic blocks containing ErbB genes that compare humans and lampreys. Shaded lines depict syntenic blocks of different ErbB genes in humans compared with that from lampreys. (D) Gene expression patterns of Efgr during the regeneration process of earthworms post-injury. Expression profiles were retrieved from the study (Shao et al., 2020). Smooth fit was performed using polyfit with df = 3. (E) Gene expression patterns of ErbB genes during lamprey embryogenesis. The gene expression profile was reanalyzed from the data (Bryant et al., 2016). Low expression levels or undetected genes are not shown. (F) Gene expression patterns of ErbB genes during rat spinal cord development from the embryonic stage to adulthood. Gene expression was retrieved from our previous study (Yang et al., 2021). Smooth fit was performed using polyfit with df = 5. Colors depict different members of the ErbB gene family.

Spatiotemporal expression of EGFR during spinal cord development

Egfr and Erbb2 from the ErbB gene family have important roles in tumor etiology and progression that are characterized by expression or activation of Egfr and altered expression of Erbb2 (Hynes and Lane, 2005; Tebbutt et al., 2013). In addition to cancer research, Egfr has been reported to be substantially upregulated in proximal nerve segments on days 4 and 7 after sciatic nerve transection (Gu et al., 2015). In this study, we focused on Egfr and its regulatory elements in the spinal cord. Egfr expression patterns during spinal cord development have not been systematically reported. Our previous RNA-seq data showed that Egfr expression was low in the early embryonic stage (E9d, E11d) and high in the later embryonic and neonatal stages (E18d–P1w). Expression decreased in adolescence and adulthood (P2w–P12w). These data were consistent with the qRT-PCR results for the developmental stages ranging from E11d to P8w (). Next, we investigated EGFR protein expression and cellular localization by immunohistochemistry. Consistent with the observation at the transcriptional level, immunostaining results showed low levels of EGFR protein at E11d (Figure and ). At E14d, we detected EGFR-positive signals located approximately in the ependymal cells around the central canal of the rat spinal cord (). An obvious increase in EGFR protein was detected at E18d and was diffusely distributed within the spinal cord tissue. A relatively strong positive signal was observed in the anterior horn of gray matter (). In the P1d–P8w rats, EGFR protein levels were relatively stable and mostly confined to the ventral horn of the spinal cord with some co-localization with neurons (possibly motoneurons) ().

Spatiotemporal dynamics of EGFR expression during spinal cord development.

(A) qRT-PCR analysis of Egfr mRNA expression during spinal cord development. Data are expressed as means ± SEM (n = 3). The relative expression of Egfr from embryo day 11 (E11d) to postnatal week 12 (P12w). Expression levels were normalized to Gapdh. (B) Immunostaining for EGFR during spinal cord development. Colors depict EGFR-positive cells (red), neural stem cells (Nestin, green), or neuron (NeuN, green), and Hoechst 33342+ nuclei (blue). Individual stains (b–d) and merged (e) staining of markers within white boxed regions (a) shown on the left. Staining was performed at different time points during spinal cord development. E11d and E14d: embryo day 11 and 14, respectively; P1d: postnatal day 1; P1w and P8w: postnatal week 1 and 8, respectively. Left (a) and right (b-e) scale bars represent 100 μm. (C) Histogram depicts the relative fluorescence intensity of EGFR. All values are presented as the mean ± SEM (n = 3). *P < 0.05 vs. E11d (Student’s t-test).

Co-expression of Egfr-related genes inferred from rat spinal cord development data

Gene co-expression network analysis aids in identifying gene modules with potential biological functions and has been widely applied in biomedical research. Based on RNA-seq data from our previous rat spinal cord development study, we first constructed co-expression modules using weighted correlation network analysis. We identified a total of 27 modules of which Egfr and Erbb3 were contained in the blue module (). Integrating co-expression modules and interaction relationships retrieved from the STRING database, we extracted Egfr-neighboring subnetworks and labeled these genes as Egfr-related () with a total of 67 genes. Gene ontology enrichment analysis of the 67 Egfr-related genes demonstrated that terms related to cell proliferation, blood vessel development, and ErbB signaling pathway were greatly enriched (). Gene expression profiles consisted of 65 genes with fragments per kilobase of exon per million (FPKM) greater than 1.0. The gene profiles showed that the majority of Egfr-related genes were upregulated in later development of the spinal cord (). We also examined the correlation (Pearson’s correlation coefficient r) of these genes with Egfr using our previous time-series transcriptome dataset from post-injury rat spinal cord (Yu et al., 2019). Forty-six out of 65 (71%) genes showed strong negative or positive correlations with Egfr expression ( and ; P < 0.05); for example, the correlation for Erbb3 was r = 0.73 (P = 2.94 × 10–17).

A co-expression module with Egfr inferred by weighted gene co-expression network analysis (WGCNA) based on an RNA-seq dataset of rat spinal cord development.

(A) Co-expressed modules inferred by the WGCNA method based on an RNA-seq dataset of rat spinal cord development and a blue module containing Egfr. (B) The subnetwork of first neighbors of Egfr and protein-protein interactions (PPIs) predicted from genes in the blue module. PPIs in the blue module were retrieved from the STRING database. (C) Gene ontology analysis of biological processes shows enrichment of genes in the subnetwork of first neighbors of Egfr. (D) Gene expression profile of genes in the subnetwork of first neighbors of Egfr. * Indicates the Pearson’s correlation (r) calculated based on the spinal cord injury dataset from Yu et al. (2019).

GPCR-mediated EGFR transactivation during spinal cord development

GPCRs have been reported to mediate transactivation of EGFR (Elliott et al., 2013; Forrester et al., 2016). GPCR-mediated transactivation of EGFR has been widely studied in cancers; however, data is limited for its role in spinal cord development. We found that five genes (Gpr153, Gpr21, Dear, Crcp, and Galr2) assigned to G protein-coupled receptor activity (GO:0004930) were clustered together with Egfr in the blue module. Previous studies have reported that some orphan receptors are involved in EGFR transactivation (Wang et al., 2010; Li et al., 2011; Fujiwara et al., 2012; Girgert et al., 2012; Ochiai et al., 2013; Cho-Clark et al., 2014). We determined that the two orphan receptor genes Gpr153 and Gpr21 showed positive correlations with the expression of Egfr (Pearson’s correlations of 0.87 and 0.56, respectively). Gpr153 was mainly expressed in the CNS (Fathi et al., 1998; Ruiz-Opazo et al., 1998). Mechanisms of GPCR-mediated EGFR transactivation include two types: the triple membrane passing signal (TMPS) and ligand-independent pathways. Our results showed that some important genes involved in the TMPS pathway were positively related to the expression of Egfr, such as Adam 1a, Adam4l1, Mmp16, Mmp17, Mmp23, Akt1, Sos2, Pik3cd, Raf1, and Map2k2 (). In the ligand-independent pathway, three genes (Src, Fyn, and Lck) belonging to the Src family also showed a positive relationship with the expression of Egfr (). The roles of these GPCR-related genes that are co-expressed with Egfr in rat spinal cord development require further experimental investigation.

Prediction of miRNAs targeting Egfr

MiRNAs are a class of small, regulatory non-coding RNAs that normally bind to the 3’-UTR of target genes and influence gene expression at the post-transcriptional level (Bushati and Cohen, 2007). There are 59 human and 17 mouse miRNAs recorded in the TarBase database that are experimentally supported. Various miRNAs targeting humans and mouse Egfr have been widely investigated, but only miR-128-3p has been reported in rats. MiRNAs generally negatively regulate gene expression of their targets. Thus, we combined three widely applied prediction tools (PITA, miRanda, and TargetScan) and calculated correlations with Egfr expression to predict candidate miRNAs targeting Egfr in rats. In total, 120 miRNAs were supported by at least two programs of which 15 miRNAs were predicted by three of the programs (). To select miRNAs that may regulate Egfr during spinal cord development, we calculated Pearson’s correlation coefficients between 121 miRNAs, including miR-128-3p, and Egfr. Finally, 25 miRNAs showed a significant negative correlation with Egfr expression (Pearson’s correlation coefficients were –0.67 to –0.58; and ). Of the 25 candidates, three miRNAs (miR-34a-5p, miR-21-3p, and let-7b-5p) have been reported to target human EGFR and five miRNAs (miR-1843a-5p and let-7a, 7b, 7c, and 7e) have been reported to target mouse Egfr according to miRNA-target pairs deposited in Tarbase. Phylogenetic analysis of the let-7 family (six members: let-7a–f) in humans, mice, and rats demonstrated that five members (let-7b, 7c, 7d, 7e, and 7f) are conserved between these three species and, except for let-7e, they are conserved in their 5’ seed regions (nucleotides 2–8) (). These data suggested that let-7 in the rat may also target Egfr. Additionally, retrieval of publications deposited in PubMed showed that 13 out of 25 miRNAs have been reported to play a role in inhibiting cell proliferation (Wang et al., 2011, 2019b, c; Kim et al., 2012, 2015; Zhao et al., 2013; Xu et al., 2014; Yang et al., 2015, 2019; Lu et al., 2016, 2017; Fang et al., 2017; Wu et al., 2017; Sun et al., 2018) in humans or mice by targeting Egfr or other targets and include miR-330-5p, miR-29c-3p, and miR-34a-5p ().

Prediction of miRNAs that target the 3’-UTR of rat .

(A) Venn diagram of predicted miRNAs targeting the 3’-UTR of rat Egfr using three programs and supported by experimentally validated data from the Tarbase database (indicated by *). Red numbers indicate miRNAs that were predicted by at least two programs or supported by experimental evidence. Fifteen miRNAs supported by three programs are shown in the dashed rectangle. (B) Expression of 25 candidate miRNAs that target Egfr in silico that were predicted based on the spinal cord development dataset. Boxed text indicates that miRNAs were reported to inhibit cell proliferation in humans or mice. MiRNAs with significant negative correlations (P < 0.05) with Egfr expression are shown. (C) Phylogenetic analysis of the let-7 family in humans, mice, and rats. Seed sequences of let-7c, 7b, 7d, 7f, and 7e are highlighted with a red line.

Potential TFs regulating Egfr in silico and confirmed in vitro

TFs are the key cellular components that control gene expression via cis-regulatory activity. We predicted potential TFs in rats that regulate Egfr expression by combining multiple methods (). A total of nine TFs were identified as candidates (). A dual-luciferase reporter assay was used to determine the targeting relationship between the nine TFs and Egfr in vitro. Six out of nine TFs showed a significant change in Egfr reporter activity compared with that in the control (P < 0.05; ). KLF1 and STAT3 significantly repressed Egfr reporter activity and FOXA2, SOX17, BHLHA15, and DDIT3 activated the Egfr reporter (P < 0.05; ). These data suggested that these six TFs may be regulators of Egfr in rats. We also used the same methods to predict TFs that regulate mouse Egfr and human EGFR. Five TFs (FOXA2, KLF4, STAT3, SP1, and KLF1) were shared by the three species (). KLF4 and STAT3 were supported by ChIP-seq results for humans deposited in the ChIPBase database, and SP1 activation of EGFR has been experimentally validated (Tsai et al., 2019).

Prediction of transcription factors that regulate .

(A) Workflow for screening potential transcription factors (TFs) that regulate Egfr expression. (B) Predicted TF binding loci within the rat Egfr promoter. (C) TF-mediated activation or repression of Egfr reporter activity using dual-luciferase assays. The reporters and effectors were co-expressed in 293T cells, and both REN and LUC activity were measured. The relative LUC activities normalized to the REN activities are shown. Results are presented as the mean ± SD from three independent experiments, each with three replicates. *P < 0.05 (one-way analysis of variance followed by the Dunnett multiple comparison test). (D) Three potential conserved transcription factors of Egfr in humans, mice, and rats. The star (*) indicates support from ChIP-seq data in humans. (E) Relative mRNA expression of six candidate TFs at embryonic days 18 (E18d) and postnatal week 1 (P1w) and 8 (P8w) was compared with the expression at embryonic day 13 (E13d). Expression levels of target genes were normalized to Gapdh. Data are expressed as means ± SEM (n = 3).

Expression of six candidate TFs during spinal cord development

We investigated whether FOXA2, SOX17, BHLHA15, DDIT3, KLF1, and STAT3 were expressed during spinal cord development at key developmental stages using qRT-PCR and immunohistochemistry. According to the expression pattern of Egfr mRNA during spinal cord development (), we selected several time points (E13d, E18d, P1w, and P8w). The qRT-PCR demonstrated that Foxa2, Klf1, Sox17, and Stat3 were relatively highly expressed during embryonic development (), which then gradually decreased, while Ddit3 expression was relatively low at E13d and increased by E18d. Bhlha15 showed relatively high expression in adult rats (P8w). The differences in expression of some TFs and their effects on EGFR promoter-luciferase activity may be caused by the detection of pooled mixed mRNA expression levels from different regions of the spinal cord. Immunostaining indicated that FOXA2 was highly expressed at E13d, and expression decreased by E18d and P1w and was undetectable at P8w (Figure and ). Undetectable positive signals for the other five TFs are not shown. FOXA2 was co-localized within the cell nucleus, and FOXA2-positive cells were found in the floor plate suggesting the importance of FOXA2 in spinal cord development, which deserves further investigation.

Spatiotemporal dynamics of FOXA2 expression during spinal cord development.

(A) Immunostaining for FOXA2 during spinal cord development. Colors depict FOXA2-positive cells (red), neural stem cells (Nestin) or neurons (NeuN) (green), and Hoechst 33342+ nuclei (blue). Individual stains (b–d) and merged staining (e) of markers within the white boxed regions shown on the left (a) at different time points during spinal cord development. Left scale bar and right scale bars (b–e) are 100 μm. FP: Floor plate. (B) Histogram depicts the relative fluorescence intensity of FOXA2. E13d and E18d: embryo day 13 and 18, respectively; P1w and P8w: postnatal week 1 and 8, respectively. All values are presented as means ± SEM (n = 3). *P < 0.05 vs. E13d (Student’s t-test).

Discussion

EGFR belongs to the receptor tyrosine kinase superfamily and plays an essential role in embryonic development and adult tissue homeostasis, including differentiation, growth, cell maintenance, and repair of damaged tissues. Egfr expression has been explored at specific stages of spinal cord development (Romano and Bucci, 2020); however, evaluation of Egfr expression patterns spanning developmental stages from neural tube formation to the adult spine has not been reported. RNA-seq and qPCR results demonstrated a bell-shaped Egfr expression pattern, which peaked near birth. Immunofluorescence staining showed that EGFR protein levels accumulated in and were confined to the ventral horn of the adult spine. EGFR showed inconsistent expression patterns at the mRNA and protein level, but FOXA2 showed consistent expression patterns. Protein abundance is affected by a series of linked processes, including transcription, post-transcriptional modification, translation, localization, post-translational modification, and protein degradation (Vogel and Marcotte, 2012). FOXA2 is a developmental transcription factor that is specifically expressed in the spine floor plate and is tightly controlled because knockdown or misexpression causes an abnormal floor plate (Bayly et al., 2012). EGFR has a complex structure with an extracellular ligand-binding domain, a single transmembrane domain, and a cytoplasmic domain. It can homo- or hetero-dimerize or form dimers with other proteins, activate multiple signal pathways, and is widely expressed in multiple organs and cells. We hypothesized that inconsistent Egfr mRNA and protein levels may be caused by protein modifications or estimates from pooled mixed mRNA expression levels from very different regions of the whole spinal cord. MiRNA-mediated gene regulation has been implicated in development and various diseases. Studies of the regulatory relationships between miRNAs and Egfr have been mainly focused on humans and mice, and only one miRNA has been reported in the rat. We predicted several miRNAs that may regulate Egfr expression during spinal cord development in rats, and some of these miRNAs have been shown to inhibit cell proliferation in other species. EGFR is widely distributed in the CNS and plays an important role in the maintenance of the pool of neural stem cells, maturation and function of astrocytes, oligodendrogenesis, and neurite growth. The contradictory effects of EGFR on nerve regeneration largely rely on the specific target cells. The changes in intracellular transport and signal transduction of EGFR are associated with the occurrence and progression of neurodegenerative diseases and injuries (Le Pichon et al., 2013; Chen et al., 2019), suggesting that regulating Egfr expression or EGFR-mediated signal transduction through upstream regulators may contribute to facilitating nerve regeneration and combating neurodegenerative diseases.

Because of the importance of the ErbB gene family, several studies (Stein and Staros, 2006; Liu et al., 2013; Barberan et al., 2016) have discussed the origin and evolution of ErbB genes in vertebrate and invertebrate animals. A vast amount of published genomic data is available to provide details regarding ErbB evolution in animals. The decoded genomes of flatworms and earthworms (Eisenia andrei) showed expansion of Egfr and revealed its relationship with regeneration (Barberan et al., 2016; Shao et al., 2020). Liu et al. (2013) analyzed many vertebrate species but lacked species from Agnatha, which are related to all vertebrates and represent the ancestral vertebrate. In this study, we added three species from Agnatha (Eptatretus burgeri, Petromyzon marinus, and Lethenteron camtschaticum) and the whale shark (Rhincodon typus) (Read et al., 2017), which is the sister group of the remaining gnathostomes. Our analysis indicated that Egfr was duplicated in ancient vertebrates via small-scale duplication but then diverged into other ErbB members after the speciation of gnathostomes. This was generally consistent with the previous study (Liu et al., 2013) but provided a more accurate node (Gnathostomata) for its divergence.

Egfr is not only involved in the development of the CNS but also participates in the biological process of CNS injury and repair (Barberan et al., 2016). It has been shown that low expression of Egfr normally occurs in astrocytes but high expression is observed after injury (Scholze et al., 2014). EGFR can activate astrocytes, which contributes to the formation of glial scars (Erschbamer et al., 2007). Understanding the regulators of Egfr expression is important in CNS development, injury, and regeneration. Transactivation of EGFR has been reported to be mediated by an expanding repertoire of GPCRs and promotes cell survival (Grisanti et al., 2017). We identified several factors that may be involved in transactivation. Several TFs that may regulate Egfr expression were identified in this study. FOXA2 opens compacted chromatin for other transcription-related proteins by displacing nucleosomal core histones and participates in embryonic development and regulation of tissue-specific gene expression (Lee et al., 2019). Our experiments showed that FOXA2 was tightly related to spinal cord embryonic development.

This study has limitations. First, functional experiments, such as RNA interference (RNAi) and over-expression, were not conducted in this study to specifically determine the roles of EGFR and other factors in spinal cord injury. Second, the expression of transcription factors screened by dual-luciferase reporter assay was not consistent with EGFR expression during spinal cord development, which may have resulted from pooling mixed expression levels from different regions of spinal cord tissue. Further studies will be required to confirm and determine the regulatory activities of the identified TFs and miRNAs on Egfr expression in vivo during spinal cord development and injury. In future studies, we plan to determine the expression of these regulators after spinal cord injury and compare expression to that during spinal cord development. Development-specific regulators will be over-expressed in the spinal cord after injury and effects on spinal cord repair will be evaluated.

In conclusion, we elucidated the long-term expression patterns of EGFR during rat spinal cord development and the evolutional complexity and diversity of the ErbB gene family in Metazoa. Initially screened candidate regulators that mediated Egfr expression were identified and require further investigation, including overexpression and RNAi studies as mentioned above. Our study provides new insight into the evolution of the ErbB gene family and regulators in Egfr expression at the transcriptional or post-transcriptional levels based on the long-term profile of spinal cord development. Given the importance of Egfr in the development and regeneration of the nervous system, the regulatory factors identified in this study may provide new candidate targets for CNS regeneration and disease treatment in the future.

Additional files:

.

Pearson's correlation with Egfr expression based on spinal cord injury dataset

.

.

Table 1

Primers for quantitative reverse transcription-polymerase chain reaction used in this study

Gene	Sequence (5’–3’)	Product size (bp)
Gapdh	Forward: ACA GCA ACA GGG TGG TGG AC	20
	Reverse: TTT GAG GGT GCA GCG AAC TT	20
Egfr	Forward: CCA CCA AGA CAG GCG ACG	18
	Reverse: AGT AGC TTG GTT CTC GCA GTC	21
Foxa2	Forward: GGT GGG TAG CCA GAA AAA GGC	21
	Reverse: CAC GGC TCC CAG CAT ACT T	19
Sox17	Forward: TTC AGC CGT CCT ATT TCC CC	20
	Reverse: CCA CCA CCT CTC CTT TCA CC	20
Bhlha15	Forward: ACA ACC TCG TGG CTC TTT GT	20
	Reverse: AGG ATC ACC CTG GTT GGA ACA	21
Ddit3	Forward: TCA TAC ACC ACC ACA CCT GAA	21
	Reverse: AGG GAT GCA GGG TCA AGA GT	20
Klf1	Forward: CAC GCA CAC CGG AGA GAA G	19
	Reverse: GAG TGG GAC CAT CCT TTG GG	20
Stat3	Forward: CCG GAG CTG CAG TTT AGA CA	20
	Reverse: CAC TGT CTC TGG GGC TGA AG	20

Additional Table 2

Key genes involved in TMPS pathway and ligand-independent pathway were positively

GeneID	E9d	E11d	E14d	E18d	P1d	P3d	P1w	P2w	P3w	P4w	P8w	P12w	Pearson Correlation Coefficient with EGFR
Adamla	4.33	2.78	4.88	6.04	3.81	3.38	3.64	2.66	2.3	1.67	2.61	1.04	0.54
Adam4l1	2.33	2.26	3.14	4.61	4.1	4.42	3.58	3.24	3.35	3.23	3.42	2.74	0.47
Mmp16	9.36	11.77	13.45	18.54	16.72	14.01	13.42	9.78	7.46	5.11	6.11	5.99	0.78
Mmp17	1.82	2.87	3.66	4.89	6.02	8.43	6.89	5.09	3.7	4.19	3.15	3.02	0.54
Mmp23	0.55	1.47	0.57	0.26	0.35	0.32	0.44	0.35	0.3	0.16	0.13	0.12	0.41
Aktl	48.48	54.18	44.42	48.47	37.4	40.96	33.96	32.54	24.56	24.73	29.88	29.26	0.58
Sos2	3.41	4.6	4.54	5.29	4.96	4.48	3.83	4.38	4.02	6.33	3.81	1.39	0.33
Pik3cd	0.7	1.58	1.75	1.6	0.81	0.76	0.86	0.66	0.49	0.4	0.83	1.44	0.33
Rafl	37.19	39.36	43.79	37.06	35.45	28.3	24.82	21.81	16.68	16.12	16.56	12.13	0.48
Map2k2	33.52	36.83	32.8	33.97	33.75	34.66	28.06	22.92	19.25	19.79	29.77	28.88	0.5
Src	7.95	7.79	10.19	5.49	8.19	8.6	7.02	5.18	4.69	3.9	4.77	4.97	0.27
Fyn	23.67	22.22	49.59	52.92	65.71	69.96	77.95	76.94	59.5	39.67	40.52	26.04	0.36
Lck	0.59	0.07	1.4	1.21	1.3	1.31	1.31	0.39	0.23	0.08	0.23	0.46	0.43
EGFR	1.14	3.18	1.29	3.67	2.07	2.84	3	1.98	1.18	0.94	1.25	0.78	-

Additional Table 3

Predicted miRNA targeting Egfr

miRNA ID	Pearson’s correlation			Expression											Prediction
	cor	pval	E9d	E11d	E14d	E18d	P1d	P3d	P1w	P2w	P3w	P4w	P8w	P12w	PITA	miRanda	TargetScan	Number of support
mo-let-7a-5p	-0.62625	0.029354	6477.333	6403.333	16388.67	21683.33	17401.33	21453	37934.67	53476.33	86457.33	118870.3	106036.3	124852.7	1	0	1	2
mo-let-7b-5p	-0.627513	0.028932	223.6667	9.333333	554.3333	654.6667	620	868.6667	1761	2123.667	3526.667	4466.667	4260	5667.333	1	1	0	2
rno-let-7c-5p	-0.611147	0.03475	7232.667	2547.667	20106.67	25500	18315	23583	39117	49975.33	75005.33	97705	85402.33	108606.3	1	1	1	3
rno-let-7d-5p	-0.608652	0.035705	263	120.3333	586	878.6667	891.6667	1520	2732	4329	5791.333	6685.667	7998	9502	1	0	1	2
rno-let-7e-5p	-0.609466	0.035391	3816	11203	8543	9899	9281	9472.333	19049	27465	43790.33	63673.33	44901.67	56654.33	1	1	1	3
rno-let-7f-5p	-0.626573	0.029246	18170.33	12439.67	38217	68662.67	73109.67	78295.67	173574.3	251508	463744.3	561346.3	472311.3	526314.3	1	0	1	2
rno-miR-140-5p	-0.590189	0.043365	10.66667	8.666667	6.333333	5.666667	3.333333	6.666667	11.33333	11.66667	11.33333	20.33333	28.33333	35	1	1	0	2
rno-miR-150-5p	-0.582893	0.046689	33.66667	22.33333	78	139.6667	108.6667	243.3333	383.3333	707.6667	937	1075	1852	2634.333	1	1	0	2
rno-miR-1843a-5p	-0.603038	0.037923	51.66667	39.33333	143.3333	246.3333	272.6667	560.6667	667	1032.667	1771.333	2236	3131.667	4342	1	1	0	2
rno-miR-1843b-5p	-0.609339	0.03544	20	18.66667	58.33333	75.66667	99.66667	240.3333	299.6667	499.3333	793.3333	950.3333	1314	1811	1	1	0	2
rno-miR-21-3p	-0.62683	0.02916	12.33333	14.33333	4.333333	1.666667	4	7	10	24.33333	33.33333	30.33333	38	73.33333	1	1	0	2
rno-miR-26b-3p	-0.590143	0.043385	40	38.33333	16	8	6.666667	9	23.33333	32.66667	32.66667	38.66667	52	48.66667	1	1	0	2
rno-miR-29c-3p	-0.586028	0.04524	2.333333	6	5.333333	12	11.33333	32.33333	37.66667	91.33333	157	229.3333	362	623.3333	1	1	0	2
rno-miR-3068-3p	-0.631561	0.027609	149.6667	171.3333	171.3333	176.6667	203	288.3333	447	614	1067.333	1051	1090.333	1280.333	1	1	0	2
rno-miR-30a-3p	-0.609427	0.035407	2825	3624.667	2524.333	1981	1386.333	2170	3178.333	4582.333	4628	5750	6397	6921	1	1	0	2
rno-miR-30c-2-3p	-0.602544	0.038123	805.3333	884.3333	485.3333	532	364	486	642.6667	883.6667	1061	1433.667	1908.667	2270.667	1	1	0	2
rno-miR-30e-3p	-0.609498	0.035379	3039	3778	2565.333	2113	1621.333	2840.667	3867	5414	5898	6972	7936	8478	1	1	0	2
rno-miR-3102	-0.579426	0.04833	0.666667	1	0	1	1	0	0.666667	0.666667	3.333333	3	2.333333	4	1	1	0	2
rno-miR-330-5p	-0.604335	0.037402	76.66667	54.66667	62	61	74.33333	196	312.3333	696.6667	702.6667	833	1303	1552.667	1	1	1	3
rno-miR-339-5p	-0.665811	0.018102	16	4.666667	2	2.333333	2.333333	3.333333	10	15.66667	16.66667	21.66667	26.33333	47.33333	1	1	1	3
rno-miR-34a-5p	-0.594402	0.041523	3.666667	1	4	3.666667	11.66667	28.66667	52.33333	125.3333	239.6667	294	441	864	1	1	0	2
rno-miR-3552	-0.646476	0.02311	0	0	0	0.333333	1	1.333333	2.666667	9	11.33333	15	14.66667	22.66667	1	1	0	2
rno-miR-448-5p	-0.580258	0.047933	1	0	0	0.666667	3.666667	6	6.666667	15	34.33333	15	18	31.33333	1	1	0	2
rno-miR-6318	-0.666447	0.017952	5.333333	4.333333	3.333333	0.333333	0.333333	0.333333	2.333333	4.666667	5.666667	6.333333	2.333333	5.666667	1	1	0	2
rno-miR-98-5p	-0.631547	0.027614	535	221	1423.333	1830.333	1497.667	2475	4337.333	7388	12424	13643.33	13461.33	14301	1	0	1	2