Background
==========

MicroRNAs (miRNAs) are a class of short (18-25 nucleotides) non-coding RNAs that regulate gene expression posttranscriptionally. Their regulatory potential relies heavily on the recognition of binding sites that are located mainly in the 3\'-untranslated regions (3\'-UTRs) of target messenger RNAs (mRNAs) \[[@B1]\]. Currently, numerous miRNAs with diverse sequences are being characterized in a wide range of species \[[@B2]\], suggesting that this small RNA molecule has a major effect on phylogeny. The importance of miRNAs is also suggested from recent research demonstrating that miRNA-guided gene regulation is involved in diverse biological functions, such as cell differentiation, development, carcinogenesis, and tumour suppression \[[@B3]-[@B6]\]. For example, phylogenetically conserved miRNAs (e.g., *let-7*, *miR-1*, *miR-124*, and *miR-125*) are involved in cell differentiation and development \[[@B7]-[@B10]\]. In this case, *let-7*regulates the expression of *RAS*proteins known as critical oncogene products \[[@B11]\]. Moreover, *miR-34*, another evolutionarily conserved miRNA, is a direct downregulator of p53 and is involved in a genetic pathway that promotes cell-cycle progression \[[@B12]\].

In recent years, more than 700 miRNAs have been identified in humans \[[@B13]\], and this number is increasing. In a recent report by Friedman et al., the expression of a large number of target genes is predicted to be regulated by miRNAs \[[@B14]\]; however, relatively few of these have been verified experimentally. To overcome this problem, a series of computational methods has been developed to predict a large number of miRNA targets; e.g., TargetScan \[[@B14]\], RNAhybrid \[[@B15]\], MicroTar \[[@B16]\], PITA \[[@B17]\], miRanda \[[@B18]\], and PicTar \[[@B19]\]. Nevertheless, these computational approaches often provide numerous target candidates with a large number of false positives because of the weak complementarity between miRNAs and 3\'-UTRs \[[@B20]\]. Recently, a phylogenetic profiling approach has been applied to overcome this limitation. For example, studies of the evolution of orthologous target sites have provided insights into the prediction of efficient miRNA targets \[[@B21]\]. As for miRNAs, many miRNA families are found among various bilaterian animals, suggesting that several miRNAs and their target genes may have co-evolved; however, these features have yet to be systematically characterized.

In this study, we hypothesized that the core regulatory relationship between miRNAs and their target genes were conserved throughout the evolution of bilaterian animals. In addition, by predicting these relationships, we sought to elucidate the core function of miRNAs in the primitive gene-regulatory network of the common bilaterian ancestor. Accordingly, we focused on five miRNAs (*let-7*, *miR-1*, *miR-124*, *miR-125/lin-4*, and *miR-34*) that are conserved among bilaterian species (*Homo sapiens*(*H. sapiens*), *Mus musculus*(*M. musculus*), *Gallus gallus*(*G. gallus*), *Drosophila melanogaster*(*D. melanogaster*), and *Caenorhabditis elegans*(*C. elegans*)) and designed a procedure to extract conserved miRNA/target-gene pairs. We extracted evolutionarily conserved miRNA/target-gene pairs based on hybridization patterns and orthologous information. In addition, we experimentally verified several candidate pairs to support our methodology. Our results suggest a functional role of three major miRNAs (*let-7*, *miR-1*, and *miR-124*) that regulated genes related to development, muscle formation, and cell adhesion. These results suggest a new role for the core function of miRNAs in the primitive gene-regulatory network of the common bilaterian ancestor.

Results and Discussion
======================

Extraction of evolutionarily conserved miRNAs among five bilaterian animals
---------------------------------------------------------------------------

To extract conserved miRNA/target-gene pairs, we chose five model species (*H. sapiens*, *M. musculus*, *G. gallus*, *D. melanogaster*, and *C. elegans*) among bilaterian animals, for which there exists a vast array of data on both miRNAs and mRNAs \[[@B13]\]. Previously, several important features were described to classify miRNAs into families. It is well known that the seed sequence (the 5\' side of the miRNA sequence) is important for interaction with the target mRNAs \[[@B22]\]. Many miRNA target prediction software programs were developed using the features of seed sequences \[[@B14]-[@B19]\]. Moreover, several features have been proposed to identify conserved miRNA families, such as conservation of the mature miRNA sequence (features of the earliest miRNA classification in miRBase) \[[@B23]\] and information on the phylogenetic relationship among miRNAs \[[@B24]\]. By focusing on these features, we proposed the following criteria for extracting well-conserved miRNA families among five species: (1) complete seed sequence matching, (2) mature miRNA sequence identity exceeding 75%, and (3) high conservation among miRNA families, considering the phylogenetic relationship among miRNAs (category I) \[[@B24]\]. Consequently, from 2,404 mature miRNA sequences, we extracted five miRNA families (*let-7, miR-1*, *miR-124*, *miR-125/lin-4*, and *miR-34*) conserved evolutionarily among the five bilaterian animals (Table [1](#T1){ref-type="table"}). The sequence identity among most of the conserved miRNA families was over 80%. In particular, the sequences of *let-7*and *miR-1*family members showed very high mature miRNA sequence identity (exceeding 90%) among all five bilaterian species, which suggests that a strong selective pressure exists for the nucleotide sequence. Huang et al. described the extraction of 15 conserved miRNA families among six bilaterian animals (*H. sapiens*, *M. musculus*, *G. gallus*, *D. rerio*, *D. melanogaster*, and *C. elegans*) based on their original classification method \[[@B24]\], which included our five miRNA families. In this manuscript, we devised more stringent criteria based on nucleotide conservation to extract highly conserved miRNA families. *Nematostella vectensis*and *Amphimedon queenslandica*, which diverged before the emergence of bilaterian animals, reportedly express various types of miRNAs \[[@B25]\]; however, none of these miRNAs are sequentially similar to the five evolutionarily conserved miRNA families found in the current study (data not shown), which suggests that these evolutionarily conserved miRNAs appeared after the divergence of bilaterian animals or were lost in *N. vectensis*and *A. queenslandica*.

###### 

List of miRNAs conserved among various bilaterian animals

  Species                                                                                 
  --------------------- ------------ --------------- ------------ ----------- ----------- ------
  ***let-7***           *let-7a*     *let-7a*        *let-7a*     *let-7*     *let-7*     95.5
  ***miR-1***           *miR-1*      *miR-1*         *miR-1a*     *miR-1*     *miR-1*     90.9
  ***miR-124***         *miR-124*    *miR-124*       *miR-124a*   *miR-124*   *miR-124*   82.6
  ***miR-125/lin-4***   *miR-125b*   *miR-125b-5p*   *miR-125b*   *miR-125*   *lin-4*     81.8
  ***miR-34***          *miR-34a*    *miR-34a*       *miR-34a*    *miR-34*    *miR-34*    79.2

Filtering and enrichment of the evolutionarily conserved miRNA/target-gene pairs
--------------------------------------------------------------------------------

To extract genes targeted through evolution by the five conserved miRNAs, we designed a procedure that comprised three screening steps (Figure [1](#F1){ref-type="fig"} and Table [2](#T2){ref-type="table"}). In step 1, we extracted potential target genes based on optimal free-energy information with a requirement of complete seed-sequence-matching using the RNAhybrid software \[[@B15]\], which predicts potential binding sites of short RNAs among target sequences. For this purpose, the optimal free energy was determined beforehand to efficiently cover experimentally validated miRNA/target-gene pairs. We calculated the free energy of 139 experimentally validated miRNA/target-gene pairs with a complete match of seed sequence, as assessed using RNAhybrid; subsequently, we defined the optimal free energy as \< -17 kcal/mole of potential miRNA/target-gene pairs (data not shown). The evaluation of target-gene extraction during each step was carried out using an index termed \"Enrichment\" (see Methods section). As a result, in step 1, the number of potential target-gene candidates decreased from 357,430 to 153,387, with an enrichment index of 1.8. This prediction set contained 112 of the 145 experimentally validated target-gene pairs.

###### 

Summary of target gene extraction after each screening step

                                             ***H. sapiens***   ***M. musculus***   ***G. gallus***   ***D. melanogaster***   ***C. elegans***   Ratio of predicted miRNA/target-gene pairs (%)   Ratio of experimentally verified miRNA/target-gene pairs (%)   Enrichment
  ------------------------------------------ ------------------ ------------------- ----------------- ----------------------- ------------------ ------------------------------------------------ -------------------------------------------------------------- ------------
  **All possible miRNA/gene pairs**          97820              98790               53995             52850                   53975              357430/357430 (100)                              145/145 (100)                                                  1.0
  **STEP 1 (Free energy and seed region)**   56666              53420               20410             14098                   8793               153387/357430 (42.9)                             112/145 (77.2)                                                 1.8
  **STEP 2 (Optimal parameter)**             12405              11256               2524              1118                    898                28201/357430 (7.9)                               76/145 (52.4)                                                  6.6
  **STEP 3**                                                    31                                    31/10356\* (0.3)        4/52 (7.7)         25.7                                                                                                            

The total number of possible target genes for the five miRNAs (*let-7*, *miR-1, miR-124, miR-125/lin-4*, and *miR-34*) is represented for the five model species at each screening step (see Figure 1). Asterisks indicate the number of orthologous genes conserved among at least four species. Enrichment was calculated as the \"ratio of experimentally verified miRNA/target-gene pairs\" divided by the \"ratio of predicted miRNA/target-gene pairs\".

![**Computational extraction of conserved miRNA/target gene pairs among bilaterian animals**. Evolutionarily conserved miRNAs were extracted from the five model species (*H. sapiens*, *M. musculus*, *G. gallus*, *D. melanogaster*, and *C. elegans*). For each miRNA, potential target genes were predicted using the following criteria: optimal free-energy threshold and complete matching of nucleotide sequences between the seed sequence of miRNA/mRNA duplexes (step 1), binding pattern of the 3\'-UTR of miRNA/mRNA duplexes (step 2), and orthologous gene information (step 3).](1471-2164-11-101-1){#F1}

In step 2, we considered the 3\'-UTR binding pattern of the miRNA/target mRNA. We defined four binding parameters (i.e., number of mismatches of mRNA within the whole miRNA sequence, number of mismatches of miRNA within the whole miRNA sequence, number of G-U wobble pairs within the whole miRNA sequence, and number of G-U wobble pairs within the seed sequence) of the hybridization pattern (Additional file [1](#S1){ref-type="supplementary-material"}) for optimisation of the thresholds for each of the features used to predict reliable miRNA/mRNA pairs. The ranges of the four binding parameters were determined by calculating the coverage of 112 miRNA-mRNA pairs verified experimentally and of 153,387 miRNA-mRNA pairs predicted for each binding feature (Figure [2A-D](#F2){ref-type="fig"}). Five hundred parameter combinations were plotted on a 2D graph, using \"Enrichment\" on the X-axis and \"Ratio of experimentally verified miRNA/mRNA\" on the Y-axis (Figure [2E](#F2){ref-type="fig"}). From these parameter combinations, we defined the optimal combination of binding parameters for efficient screening based on a maximum EC value of 184.7 (number of mismatches of mRNA within the whole miRNA sequence: 12; number of mismatches of miRNA within the whole miRNA sequence: 10; number of G-U wobble pairs within the whole miRNA sequence: 4; and number of G-U wobble pairs within the seed sequence: 0). Accordingly, the number of potential target-gene candidates was reduced from 153,387 to 28,201 and the experimentally validated target genes decreased from 112 to 76 after introduction of the criterion of optimal hybridization pattern (Enrichment index, 6.6) (Table [2](#T2){ref-type="table"}).

![**Parameters used for the prediction of miRNA/mRNA pairs and their coverage**. To optimize the binding parameters of miRNA/mRNA duplexes, we determined the coverage of four binding parameters (mismatch of mRNA within the whole miRNA sequence (A), mismatch of miRNA within the whole miRNA sequence (B), G-U wobble pairs within the whole miRNA sequence (C), and G-U wobble pairs within the seed sequence (D)). Calculation of the coverage was performed using 112 experimentally verified miRNA/mRNA pairs (orange line) and 153,387 predicted miRNA/mRNA pairs (blue line). The yellow bar indicates differences in coverage between experimentally verified and computationally predicted miRNA/mRNA pairs. Four or five points chosen from the highest yellow bar were used as the range of each of the four binding parameters (blue squares) used in this study. We determined parameter space using the binding patterns of the miRNA/mRNA pairs based on four features (E). Five hundred parameter combinations were plotted on a 2D graph using \"ratio of experimentally verified miRNA/mRNA\" on the Y-axis and \"Enrichment\" on the X-axis. Black circles (57.3% of the coverage and 3.2-fold of the Enrichment) indicate the point that corresponded to optimized parameters for the prediction of final conserved miRNA/target pairs: 12 for the number of mismatches in the mRNA, 10 for the number of mismatches in the miRNA, 4 for the number of G-U wobble pairs within the whole miRNA sequence, and 0 for the number of G-U wobble pairs within the seed sequence (see Methods section).](1471-2164-11-101-2){#F2}

Finally, in step 3, we incorporated orthologous gene information and extracted genes that were evolutionarily conserved among more than four diverse bilaterian animals, including *H. sapiens*. As a result, the number of predicted miRNA/target-gene pairs was minimized substantially, from 10,356 to 31, using a significantly high Enrichment index of 25.7 (Table [2](#T2){ref-type="table"}). The number of predicted miRNA/target-gene pairs was especially high for the three miRNAs *let-7*(eight targets), *miR-1*(seven targets), and *miR-124*(eleven targets) compared with *miR-125/lin-4*(three targets) and *miR-34*(two targets) (Additional file [2](#S2){ref-type="supplementary-material"}). This suggests that *let-7*, *miR-1*, and *miR-124*may have played a major role in primordial miRNA gene regulation in the common bilaterian ancestor. To verify the significance of conserved miRNA-target gene pairs, we performed same sequence analysis (from step 1 to step 3) against total 25 species-specific miRNAs (5 miRNAs each from 5 species) as a control experiment supposing that these miRNAs are also conserved in other bilaterians. For example, target prediction of hsa-miR-2277, a species-specific miRNA in human was performed in all 5 species (step 1 and step 2) and conserved targets were extracted (step 3). As a result, 11 out of 25 non-conserved miRNAs did not show any conserved miRNA-target gene pair. Furthermore, average number of the miRNA-target gene pairs of the negative control was 2.4, which is statistically lower than that of conserved miRNA-target gene pairs 6.2 based on the Welch\'s t-test (P \< 0.05). These results support that number of genes achieved from the prediction of conserved miRNAs target genes in this study is indeed significant. In summary, we developed a new filtering method for extracting evolutionarily conserved miRNA/target-gene pairs, which was used to extract 31 reliable miRNA/target-gene pairs among the five families of miRNAs.

We discovered that only one orthologous target gene, calponin-3 (*CNN3*), was conserved completely among the five bilaterian animals (Figure [3A](#F3){ref-type="fig"}). As for target genes conserved in four species, we found, for example, the La-related protein 4 (*LARP4*), ETS domain-containing protein Elk-3 (*ELK3*), argonaute-4 (*EIF2C4*), transgelin-2 (*TAGLN2*), and V-type proton ATPase subunit B brain isoform (*ATP6V1B2*) genes (Figure [3B](#F3){ref-type="fig"} and Additional file [3](#S3){ref-type="supplementary-material"}). Of note, the same approximate position of the predicted target site was observed in the orthologous 3\'-UTR of *CNN3*(120 nucleotides (nt)), *LARP4*(3,100 nt), *EIF2C4*(220 nt), and *TAGLN2*(50 nt) among vertebrates (Figure [3](#F3){ref-type="fig"} and Additional file [3](#S3){ref-type="supplementary-material"}). According to Bartel et al., the distribution of miRNA target sites within the 3\'-UTR is biased near the mRNA stop codon or poly-A tail compared with the middle portion of 3\'-UTR \[[@B26]\]. Our results show that target site distribution varied according to the type of miRNA target gene. The target sites on 3\'-UTR of *CNN3*and *TAGLN2*were biased near the stop codon, from *H. sapiens*to *C. elegans*(Figure [3A](#F3){ref-type="fig"} and Additional file [3](#S3){ref-type="supplementary-material"}). Regarding the other candidates, we observed all types of target site distribution on 3\'-UTR. A future statistical analysis of miRNA target-site distribution among conserved miRNA/target-gene pairs is required to substantiate this view. With the exception of 3\'-UTR of the *LARP4*gene, most of the binding patterns of evolutionarily conserved target sites were sequentially different, without taking the seed region into consideration. The target-site binding patterns within 3\'-UTR of the *LARP4*gene were identical between *H. sapiens*and *M. musculus*(Figure [3B](#F3){ref-type="fig"}), although the similarity of the two 3\'-UTR sequences was \~70% (data not shown). A recent study reported on cooperative regulation by multiple miRNAs \[[@B19]\]. Likewise, the band 4.1-like protein 4B (*EPB41L4B*) gene was an orthologous target of two different types of miRNAs: *miR-1*and *miR-124*(Table [3](#T3){ref-type="table"}). Our analysis suggests that multiple miRNA regulation may have already existed in the era of ancestral bilaterian species.

###### 

Evolutionarily conserved genes regulated by miRNAs

  Target gene                                                                                                                                                                  
  -------------------------------------------- ---------- -------------------------------------------------------------------------- ------- ----------- ----------- --------- -------
  **Candidates targeted by multiple miRNAs**                                                                                                                                   
                                                                                                                                                                               
  ENSG00000095203                              EPB41L4B   Band 4.1-like protein 4B                                                           **+\*\***   **+**                 
                                                                                                                                                                               
  **Candidates targeted by single miRNAs**                                                                                                                                     
                                                                                                                                                                               
  ENSG00000187772                              LIN28B     Lin-28 homolog B                                                           **+**                                     
  ENSG00000198799                              LRIG2      Leucine-rich repeats and immunoglobulin-like domains protein 2 precursor   **+**                                     
  ENSG00000086544                              ITPKC      Inositol-trisphosphate 3-kinase C                                          **+**                                     
  ENSG00000196233                              LCOR       Ligand-dependent corepressor                                               **+**                                     
  ENSG00000139263                              LRIG3      Leucine-rich repeats and immunoglobulin-like domains protein 3 precursor   **+**                                     
  ENSG00000134698                              EIF2C4     Eukaryotic translation initiation factor 2C 4                              **+**                                     
  ENSG00000170456                              DENND5B    MGC24039 protein                                                           **+**                                     
  ENSG00000136231                              IGF2BP3    Insulin-like growth factor 2 mRNA-binding protein 3                        **+**                                     
  ENSG00000158710                              TAGLN2     Transgelin-2                                                                       **+\***                           
  ENSG00000143549                              TPM3       Tropomyosin alpha-3 chain                                                          **+\***                           
  ENSG00000117519                              CNN3       Calponin-3                                                                         **+**                             
  ENSG00000071073                              MGAT4A     Alpha-1,3-mannosyl-glycoprotein 4-beta-N-acetylglucosaminyltransferase A           **+**                             
  ENSG00000161813                              LARP4      La-related protein 4                                                               **+**                             
  ENSG00000147416                              ATP6V1B2   Vacuolar ATP synthase subunit B, brain isoform                                     **+**                             
  ENSG00000135862                              LAMC1      Laminin subunit gamma-1 precursor                                                              **+\***               
  ENSG00000125695                              STRADA     STE20-related adapter protein                                                                  **+**                 
  ENSG00000131459                              GFPT2      Glucosamine\--fructose-6-phosphate aminotransferase                                            **+**                 
  ENSG00000111145                              ELK3       ETS domain-containing protein Elk-3                                                            **+**                 
  ENSG00000151726                              ACSL1      Long-chain-fatty-acid\--CoA ligase 1                                                           **+**                 
  ENSG00000164144                              ARFIP1     Arfaptin-1                                                                                     **+\*\***             
  ENSG00000080819                              CPOX       Coproporphyrinogen III oxidase, mitochondrial precursor                                        **+**                 
  ENSG00000138279                              ANXA7      Annexin A7                                                                                     **+**                 
  ENSG00000150093                              ITGB1      Integrin beta-1 precursor                                                                      **+\*\***             
  ENSG00000093167                              LRRFIP2    Leucine-rich repeat flightless-interacting protein 2                                           **+**                 
  ENSG00000116141                              MARK1      Serine/threonine-protein kinase MARK1                                                                      **+**     
  ENSG00000166797                              FAM96A     Protein FAM96A                                                                                             **+**     
  ENSG00000131914                              LIN28      Lin-28 homolog A                                                                                           **+\***   
  ENSG00000113282                              CLINT1     EPN4_HUMAN Isoform 2 of Q14677 - Homo sapiens                                                                        **+**
  ENSG00000137872                              SEMA6D     Semaphorin-6D precursor                                                                                              **+**

List of conserved targets regulated by conserved miRNAs. Transcript of Band 4.1-like protein 4B is only regulated by two miRNAs (*miR-1*and *miR-124*). \"\*\" and \"\*\*\" indicates experimentally verified miRNA-mRNA and possible miRNA-mRNA candidates revealed by microarray data, respectively.

![**Two examples of miRNA target sites in orthologous gene transcripts**. Potential target sites of *miR-1*(blue arrows) in the 3\'-UTR sequences of the orthologous *CNN3*(A) and *LARP4*(B) transcripts are shown. Predicted duplexes formed by the 3\'-UTR sequences (top) and miRNAs (bottom) are shown in dotted boxes for each potential target site. The green bar on the *H. sapiens*3\'-UTR sequence indicates a DNA region used for the construction of the reporter plasmid pLuc-CNN3 (Figure 4A). See Additional file [3](#S3){ref-type="supplementary-material"} for other candidates. (\*) The length of the CNN3 3\'-UTR is currently registered as a little shorter than that indicated (527 nt in size) and contains the *miR-1*binding site (Ensembl release 53).](1471-2164-11-101-3){#F3}

![**Example of the 3\'-UTR reporter plasmid and experimental validation**. The 3\'-UTR sequences of *DLL1, CNN3, LARP4, ELK3, EIF2C4, TAGLN2*, and *ATP6V1B2*were subcloned into the *Xho*I/*Not*I site of the psiCHECK™-2 vector. *CNN3*was chosen as representative of the eight candidates listed above (see Methods section). (B) HeLa cells were cotransfected with each combination of 100 ng of reporter plasmid and the indicated amounts of each miRNA (*DLL1*, 5 pmol of *miR-34*(a); *CNN3*, 60 pmol of *miR-1*(b); *LARP4*, 20 pmol of *miR-1*(c); *ELK3*, 60 pmol of *miR-124*(d); *EIF2C4*, 60 pmol of *let-7*(e), *TAGLN2*, 5 pmol of *miR-1*(f), *ATP6V1V2*, 60 pmol of *miR-1*(g), and *LARP4*, 60 pmol of *let-7*(h)). Colours depict each miRNA: *miR-34*(red), *miR-1*(blue), *miR-124*(orange), *let-7*(green), and negative control (black). The relative expression of the luciferase gene was measured 24 h after transfection. The normalized luciferase activity of the control vector was set as 1.0. The data represent the average of three experiments and SDs. \* *P*\< 0.01.](1471-2164-11-101-4){#F4}

Experimental validation of miRNA target genes
---------------------------------------------

To validate the evolutionarily conserved miRNA/target-gene candidates, we performed transfection and luciferase assays on 6 of the 31 identified evolutionarily conserved miRNA/target-gene candidates. Four of the six candidates (*CNN3, LARP4, TAGLN2*, and *ATP6V1B2*) were predicted to be regulated by *miR-1*, while one candidate (*ELK3*) was predicted to be targeted by *miR-124*, and the final candidate (*EIF2C4*) was predicted to be the target of *let-7*. We used the well-established downregulation of the delta-like protein 1 precursor (*DLL1*) gene by *miR-34*as a positive control \[[@B22]\], and the *LARP4*/*let-7*pair, which was extracted up to step 2 in our analysis, was chosen as a non-evolutionarily conserved pair. We subcloned the 3\'-UTR sequence downstream from the *Renilla*luciferase gene (Figure [4A](#F4){ref-type="fig"}) and co-transfected 100 ng of the 3\'-UTR reporter construct into HeLa cells using 5, 20, and 60 pmol of miRNA Mimics (hsa-let-7a, hsa-miR-1, hsa-miR-124, hsa-miR-34a, and miRIDIAN microRNA Hairpin Inhibitor Negative Control \#1) (see Methods section). We observed the downregulation of six out of six candidates and of the positive control compared with the negative controls (Figure [4B, a-g](#F4){ref-type="fig"}). Typical results of the reporter gene assay are shown in Figure [4B](#F4){ref-type="fig"} for the indicated amounts of miRNAs (5, 20, and 60 pmol). The downregulation of these candidates was significant (*P*\< 0.01), although some of these pairs represented an inhibition of only 30% under the current conditions. Among these candidates, *TAGLN2*was previously suggested to be downregulated by *miR-1*, as assessed by microarray analysis \[[@B27]\]. This feature was recently confirmed using the \"pulsed stable isotope labelling with amino acids in cell culture\" (pSILAC) method and a reporter gene assay \[[@B28]\]. Regarding the *let-7*/*LARP4*combination, the expression of *LARP4*was not downregulated after *let-7*transfection, which was supported statistically (Figure [4B, h](#F4){ref-type="fig"}). These experimental results suggest that our new method has the potential for efficiently extracting reliable miRNA/target-gene pairs and may be effective in the elucidation of the primordial regulatory relationships between miRNAs and their target genes during the early stage of bilaterian evolution.

Possible regulation of evolutionarily conserved miRNA targets in bilaterian animals
-----------------------------------------------------------------------------------

To provide further insight into the primary functions of evolutionarily conserved miRNAs (Table [3](#T3){ref-type="table"}), we next focused on the functions of the target genes and found that evolutionarily conserved miRNA/target genes could be largely classified into four functional categories: development, differentiation, muscle movement, and gene regulation. First, we describe the function of evolutionarily conserved genes involved in development and differentiation. The laminin subunit gamma-1 precursor (*LAMC1*) gene, which was possibly regulated by *miR-124*, is one of the major components of the basement membrane. According to Smyth et al. \[[@B29]\], null mutation of *LAMC1*causes embryonic lethality because of the absence of the basement membrane and failure to differentiate the endoderm. Among other candidates regulated by *miR-124*, the expression of the leucine-rich repeat flightless-interacting protein 2 (*LRRFIP2*) gene induces an extra axis in *Xenopus laevis*embryos \[[@B30]\]. Moreover, the semaphorin-6D precursor (*SEMA6D*) gene, plays an important role in cardiac morphogenesis during chick embryonic dev which was predicted as a candidate of *miR-34*targeting, elopment \[[@B31]\]. The *lin-28*gene, which regulates developmental timing in *C. elegans*, is reportedly controlled by *lin-4*, as assessed using *in vivo*experiments \[[@B32]\]. Subsequently, the *lin-28*gene was also found to be regulated by *miR-125*, which is an orthologous miRNA of *lin-4*, in *H. sapiens*and *M. musculus*\[[@B33]\]. Our prediction confirmed the regulation of orthologous *lin-28*genes by *lin-4/miR-125*miRNA in *H. sapiens*, *M. musculus*, and *C. elegans*and further suggested that a similar regulatory relationship was conserved in *G. gallus*. Interestingly, our prediction showed that the orthologous *let-7*/*lin-28*pair was also evolutionarily conserved among bilaterians. This feature had been validated experimentally in *H. sapiens*\[[@B34]\]; however, in *C. elegans*, *lin-28*is expressed in an early stage of development, while *let-7*is expressed in a later stage of development, suggesting that regulation of *lin-28*by *let-7*is subtle in *C. elegans*\[[@B7]\]. It would be of interest to analyse the interaction of the *let-7*/*lin-28*pair in other species, such as *M. musculus*and *G. gallus*, to understand whether stage-specific expression of *let-7*was present in the common ancestor of bilaterian animals or if it is a trait acquired later during evolution, as many of the evolutionarily conserved miRNA/target-gene pairs were related to an essential function involved in differentiation and development.

Next, we focused on the tissue-specific miRNA/target genes. Among the 30 evolutionarily conserved target candidates, approximately one-third were expressed in a tissue-specific manner in humans, according to the BioGPS portal <http://biogps.gnf.org>. *miR-1*is highly expressed in muscle tissues \[[@B8]\]. Here, three candidate genes regulated by *miR-1*(i.e., *TAGLN2*, *CNN3*, and *TPM3*) are also expressed in muscle tissues, according to BioGPS. *TAGLN2*is a homolog of *TAGLN*, which encodes an actin-binding protein and is a diagnostic marker of breast and colon carcinoma in humans \[[@B35]\]. The *CNN3*gene also encodes an actin-binding protein that represses bone morphogenetic protein (BMP) signalling in chondrocytes, which is important for bone formation \[[@B36]\]. Moreover, the *TPM3*gene encodes yet another actin-binding protein that modulates muscle contraction. The other miRNA, *miR-124*, is expressed in the nervous system \[[@B9]\]. Similarly, one of the evolutionarily conserved target candidates, the glucosamine-fructose-6-phosphate aminotransferase 2 (*GFPT2*) gene, is expressed in the central nervous system \[[@B37]\]. In addition, the annexin A7 (*ANXA7*) gene, which is another candidate target of *miR-124*, is involved in the development of the murine brain \[[@B38]\]. The striking overlap between the tissue specificity of evolutionarily conserved miRNA and that of their target genes suggests that one of the main functions of primordial miRNAs may have been the regulation of genes implicated in the temporary control of the development of muscle and of the nervous system, in a tissue-specific manner.

Finally, we found two interesting candidate genes, *EIF2C4*and *LARP4*, which encode translation-related proteins. It is well accepted that miRNAs are regulators of gene expression, mostly at the translational level \[[@B1]\]. *EIF2C4*is also known as Argonaute 4 (*AGO4*). Although the function of AGO4 is unknown, other Argonaute protein family members are involved in the RNA-induced silencing complex (RISC), which is essential for the miRNA or siRNA pathways. A previous microarray analysis performed in HepG2 cells revealed that the *EIF2C4*gene was affected by *let-7*\[[@B6]\]. In the present study, we demonstrated for the first time the direct downregulation of *EIF2C4*by *let-7*, as assessed using a reporter gene assay in HeLa cells (Figure [4B, e](#F4){ref-type="fig"}); therefore, we speculate that negative-feedback regulation of *EIF2C4*by *let-7*exists in the miRNA pathway. It has been reported that regulation of the *AGO1*mRNA, which is a major component of the RISC complex, in the miRNA pathway by *miR-168*controls plant development in *Arabidopsis thaliana*\[[@B39]\]. Another candidate, *LARP4*, encodes a member of La-motif protein family that controls translational efficiency \[[@B40]\]. We also demonstrated the downregulation of the *LARP4*gene via *miR-1*using a reporter gene assay (Figure [4B, c](#F4){ref-type="fig"}), which further supports our contention that some of the evolutionarily conserved miRNAs may play an important role in the regulation of translation by controlling the expression levels of translation factors and by negatively regulating their own miRNA pathway.

Conclusions
===========

We developed a procedure to extract potential evolutionarily conserved miRNA/target-gene pairs based on orthologous gene information from five bilaterian animals, and efficiently extracted 31 evolutionarily conserved miRNA/target-gene pairs from 357,430 pairs. We experimentally validated the downregulation of six candidate pairs (out of six tested pairs) in HeLa cells, which suggests that our method using orthologous information was efficient in extracting evolutionarily conserved miRNA target-gene candidates. Our findings reveal that miRNA target sites were conserved among various species, and demonstrate that especially *let-7*-, *miR-1*-, and *miR-124*-mediated gene regulation may have played an important role throughout evolution, in processes such as development, differentiation, and muscle movement. Moreover, our results indicate that miRNA-mediated translational regulation as well as tissue-specific expression of miRNA/target-gene pairs may have already existed in the common bilaterian ancestor. In conclusion, our study will provide new insights into the early stages of miRNA function.

Methods
=======

miRNA and 3\'-UTR sequence data
-------------------------------

We downloaded 2,404 mature miRNA sequences (885 for *H. sapiens*, 689 for *M. musculus*, 520 for *G. gallus*, 153 for *D. melanogaster*, and 157 for *C. elegans*) from the miRBase, version 13.0 <http://microrna.sanger.ac.uk/sequences/>\[[@B13]\]. We downloaded sequences corresponding to 3\'-UTR using the Ensembl transcript ID annotation in FASTA format (40,498 transcripts for *H. sapiens*, 3,332 for *M. musculus*, 13,089 for *G. gallus*, 16,822 for *D. melanogaster*, and 13,560 for *C. elegans*) from the Ensmart database Ensembl release 53 <http://www.ensembl.org/index.html>\[[@B41]\]. Orthologous gene information was also downloaded from the Ensmart database Ensembl release 53. We obtained 145 experimentally verified miRNA/target-gene pairs from TarBase Version 5.0.1 <http://diana.cslab.ece.ntua.gr/tarbase/>\[[@B42]\].

Identification of miRNAs conserved among bilaterian animals
-----------------------------------------------------------

The 2,404 miRNA sequences were aligned using ClustalX \[[@B43]\] with the following alignment parameters: gap opening, 22.50; gap extension, 0.83; and bootstrap value, 100. We checked the conservation of 2,404 miRNA sequences to extract evolutionarily conserved miRNAs. We defined the conservation threshold as an \"overall sequence identity \> 75% with complete matching of the seed sequence (1-7, 2-8, or 3-9 nucleotides from the miRNA\'s 5\' end)\"; furthermore, we introduced information on the phylogenetic relationship among miRNAs to extract reliably conserved miRNAs and used highly conserved miRNA families (category I) \[[@B24]\].

Extraction of evolutionarily conserved miRNA/target-gene pairs among bilaterian animals
---------------------------------------------------------------------------------------

To extract evolutionarily conserved miRNA/target-gene pairs among bilaterian animals, we devised a three-step filtering approach (Figure [1](#F1){ref-type="fig"}). In step 1, we predicted genes targeted by each of the five miRNAs (*let-7*, *miR-1*, *miR-124*, *miR-125/lin-4*, and *miR-34*) using RNAhybrid, which is fast and flexible software for miRNA target prediction, with the free-energy option and the seed-sequence option \[[@B15]\]. The RNA duplex free-energy filter was defined as the appropriate value that led to the efficient extraction of experimentally verified miRNA/target-gene pairs. We also considered a complete match across the seed sequence (1-7, 2-8, or 3-9 nucleotides from the miRNA\'s 5\' end), which was used as a filter by adding the seed option of RNAhybrid.

In step 2, we used four binding parameters of the hybridization pattern of the miRNA/mRNA duplexes. According to a recent study, a binding rule is likely to exist for the recognition of target mRNAs by miRNAs \[[@B34]\]. Moreover, G-U wobble pairs within miRNA/mRNA duplexes play a key role in the interaction with target mRNAs \[[@B44]\]. Subsequently, potential candidates were extracted using four binding parameters (number of mismatches of mRNA within the whole miRNA sequence, number of mismatches of miRNA within the whole miRNA sequence, number of G-U wobble pairs within the whole miRNA sequence, and number of G-U wobble pairs within the seed sequence) of the miRNA/mRNA duplexes (Additional file [1](#S1){ref-type="supplementary-material"}). We used the hybridization pattern of experimentally verified and predicted miRNA/mRNA pairs to calculate coverage, by changing these binding parameters one by one (28-0 for the number of mismatches of mRNA within the whole miRNA sequence; 14-0 for the number of mismatches of miRNA within the whole miRNA sequence; 10-0 for the number of G-U wobble pairs within the whole miRNA sequence; and 6-0 for the number of G-U wobble pairs within the seed sequence) (Figure [2A-D](#F2){ref-type="fig"}). The range of each of the four binding parameters was determined based on the coverage of experimentally verified and predicted miRNA/mRNA pairs. Parameter combinations were then plotted on a 2D graph by calculating the \"ratio of experimentally verified miRNA/mRNA\" and \"Enrichment\" using the points in the four parameter ranges (Figure [2E](#F2){ref-type="fig"}). The criterion \"Enrichment\" was defined and calculated as the value of \"ratio of experimentally verified miRNA/target-gene pairs\" divided by the \"ratio of predicted miRNA/target-gene pairs\". We obtained the most effective combination of four binding parameters for extracting miRNA/target-gene pairs based on the EC value. Parameter conbination with highest EC value was selected.

In step 3, orthologous gene information was used to extract orthologous genes targeted by the same type of miRNA. For the retrieval of evolutionarily conserved miRNA/target-gene pairs from various bilaterian animals, we set the orthologous gene information criteria as orthologous genes conserved in at least four species, each containing the miRNA target site of interest.

Expression vectors
------------------

To construct target-site reporter plasmids, each DNA fragment (3\'-UTR sequence of the *DLL1*gene (668 nt; accession no. AF003522), *ELK3*gene (519 nt; accession no. BC017371), *EIF2C4*gene (2148 nt; accession no. AB046787), *TAGLN2*gene (1391 nt; accession no. D21261), *LARP4*gene (1678 nt; accession no. AY004310), *CNN3*gene (1391 nt; accession no. BC025372), and *ATP6V1B2*gene (1208 nt; accession no. L35249)) was amplified from HeLa genomic DNA via polymerase chain reaction using site-specific primers and was inserted into the *Xho*I/*Not*I sites of the psiCHECK-2 plasmid vector (which encodes both firefly and *Renilla*luciferases; Promega, Madison, WI, USA) (Figure [4A](#F4){ref-type="fig"}). The oligonucleotide was designed to introduce *Xho*I and *Not*I sites at the 5\' and 3\' termini, respectively. The resulting plasmids were termed pLuc-DLL1, pLuc-CNN3, pLuc-LARP4, pLuc-ELK3, pLuc-EIF2C4, pLuc-TAGLN2, and pLuc-ATP6V1B2, respectively. miRIDIAN™ miRNA Mimic for hsa-let-7a, hsa-miR-1, hsa-miR-124, hsa-miR-34a, and negative control (miRIDIAN microRNA Hairpin Inhibitor Negative Control \#1) were purchased from Dharmacon. miRNA Mimic molecules are chemically modified double-stranded RNA oligonucleotides.

The sequences of the oligodeoxyribonucleotides used for PCR were as follows:

\(1\) DLL1_S: 5\'-TTACTCGAGAATGGAAGTGAGATGGCAAGAC-3\'

\(2\) DLLl_A: 5\'-TTAGCGGCCGCTTGTCATTCATAAAATTTATTT-3\'

\(3\) CNN3_S: 5\'-TTACTCGAGTCCACACAGAAGGAGCTCAG-3\'

\(4\) CNN3_A: 5\'-TTAGCGGCCGCCAGGAAGAGCAAATGCATCA-3\'

\(5\) LARP4_S: 5\'-TTACTCGAGGTTCCCATTTGATGGCATGT-3\'

\(6\) LARP4_A: 5\'-TTAGCGGCCGCCATAGCACCTTGGCGATGTT-3\'

\(7\) ELK3_S: 5\'-TTACTCGAGCGTCTGGCCACAATTAAGGA-3\'

\(8\) ELK3_A: 5\'-TTAGCGGCCGCTGCTTTCATATTGCCCACTG-3\'

\(9\) EIF2C4_S: 5\'-TTACTCGAGTTCTACCAGCAGCTCGGAAT-3\'

\(10\) EIF2C4_A: 5\'-TTAGCGGCCGCTTGGATTCCAGCAAGTCCTC-3\'

\(11\) TAGLN2_S: 5\'-TTACTCGAGCCCTCCCACGAATGGTTAAT-3\'

\(12\) TAGLN2_A: 5\'-TTAGCGGCCGCATGGAAAATGAGAAGCCACG-3\'

\(13\) ATP6V1B2_S: 5\'-TTACTCGAGTCCGCGCTCTTGTGAAATAC-3\'

\(14\) ATP6V1B2_A: 5\'-TTAGCGGCCGCATAATCATGCTGACTCCCCC-3\'

Transfection and luciferase reporter assay
------------------------------------------

Transient transfection and luciferase assays were performed as described previously, with slight modifications \[[@B22]\]. Briefly, HeLa cells were grown in 10% FBS in DMEM and seeded in 24-well plates 24 h before transfection. Cells were transfected with the indicated amounts of reporter and miRNA Mimic (100 ng of target reporter and 5, 20, and 60 pmol of miRNA Mimic) in the presence of Lipofectamine 2000 (Invitrogen, Carlsbad, CA, USA). Firefly and *Renilla*luciferase activities were measured consecutively using the Dual-luciferase assay system (Promega) 24 h after transfection, according to the manufacturer\'s instructions.

List of abbreviations used
==========================

miRNA: microRNA; 3\'-UTR: 3\'-untranslated region; EC value: multiplied value of Enrichment and Coverage; nt: nucleotides; SD: standard deviation.

Authors\' contributions
=======================

KT participated in all aspects of the study. KF contributed to the overall conception of this work and performed the functional profiling of miRNA-targeted genes. YW supported the design of the target-gene extraction procedure and validated new miRNA target candidates using transfection and luciferase reporter assay experiments. AS constructed the experimental system. NS participated in the validation of new miRNA target candidates using transfection and luciferase reporter assay experiments. MT supervised the project. AK participated in the experimental design and drafted the manuscript. All authors read and approved the final manuscript.

Supplementary Material
======================

###### Additional file 1

**Basic concept of miRNA/mRNA duplex formation**. An example of the binding pattern of miRNA (bottom)/mRNA (top) duplexes is shown using *cel-let-7*and die-1 3\'-UTR sequences from *C. elegans*. The green and blue squares depict mismatched nucleotide sequences of the mRNA and miRNA, respectively. The red square depicts G-U wobble pairs within the whole miRNA sequence and the black arrow pinpoints a G-U wobble pair within the seed sequence.

###### 

Click here for file

###### Additional file 2

Summary of the number of target genes in each extraction step

###### 

Click here for file

###### Additional file 3

**Additional examples of miRNA target sites in orthologous gene transcripts used for experimental verification**. Potential target sites of *miR-124*(orange arrows) in the 3\'-UTR sequences of orthologous *ELK3*transcripts (A). Potential target sites of *let-7*(green arrows) in the 3\'-UTR sequences of orthologous *EIF2C4*transcripts (B). Potential target sites of *miR-1*(blue arrows) in the 3\'-UTR sequences of orthologous *TAGLN2*transcripts (C) and *ATP6V1B2*transcripts (D). 3\'-UTR sequences and miRNAs are shown in dotted boxes for each potential target site; the colours of dotted boxes and arrows correspond to those of each miRNA.

###### 

Click here for file

Acknowledgements
================

We are grateful to Mizuki Komasa, Motomu Matsui, and Hikaru Taniguchi (Keio University, Japan) for useful discussions. This research was supported in part by a Grant-in-Aid from the 21st Century Centre of Excellence (COE) Program entitled \"Understanding and Control of Life\'s Function via Systems Biology\" (Keio University); research funds from the Yamagata Prefectural Government and Tsuruoka City in Japan; and a grant from Keio University.
