Characterisation and expression of microRNAs in developing wings of the neotropical butterfly Heliconius melpomene

Background Heliconius butterflies are an excellent system for studies of adaptive convergent and divergent phenotypic traits. Wing colour patterns are used as signals to both predators and potential mates and are inherited in a Mendelian manner. The underlying genetic mechanisms of pattern formation have been studied for many years and shed light on broad issues, such as the repeatability of evolution. In Heliconius melpomene, the yellow hindwing bar is controlled by the HmYb locus. MicroRNAs (miRNAs) are important post-transcriptional regulators of gene expression that have key roles in many biological processes, including development. miRNAs could act as regulators of genes involved in wing development, patterning and pigmentation. For this reason we characterised miRNAs in developing butterfly wings and examined differences in their expression between colour pattern races. Results We sequenced small RNA libraries from two colour pattern races and detected 142 Heliconius miRNAs with homology to others found in miRBase. Several highly abundant miRNAs were differentially represented in the libraries between colour pattern races. These candidates were tested further using Northern blots, showing that differences in expression were primarily due to developmental stage rather than colour pattern. Assembly of sequenced reads to the HmYb region identified hme-miR-193 and hme-miR-2788; located 2380 bp apart in an intergenic region. These two miRNAs are expressed in wings and show an upregulation between 24 and 72 hours post-pupation, indicating a potential role in butterfly wing development. A search for miRNAs in all available H. melpomene BAC sequences (~ 2.5 Mb) did not reveal any other miRNAs and no novel miRNAs were predicted. Conclusions Here we describe the first butterfly miRNAs and characterise their expression in developing wings. Some show differences in expression across developing pupal stages and may have important functions in butterfly wing development. Two miRNAs were located in the HmYb region and were expressed in developing pupal wings. Future work will examine the expression of these miRNAs in different colour pattern races and identify miRNA targets among wing patterning genes.


Background
Neotropical butterflies of the genus Heliconius provide striking examples of both divergence and convergence in their wing colour patterns. Distributed throughout the tropical forests of central and southern America, they signal their distastefulness to predators through brightly coloured wings. Many species take part in Müllerian mimicry 'rings', where multiple species converge in wing pattern and thereby benefit through protection from predators [1]. Wing patterns are also used in courtship and mate recognition, meaning they are both adaptive and contribute to genetic isolation and speciation [2,3]. Genetic crosses have shown that most phenotypic wing pattern and colour variation is controlled by a few Mendelian loci [4,5]. For example, in H. melpomene, genes in linkage group 15 control the yellow and white pattern elements (HmYb/Sb/N) and genes in linkage group 18 control red pattern elements (HmB/D) [6][7][8]. The HmYb locus, which controls the presence or absence of a hindwing yellow bar, is orthologous to Cr in the mimetic species H. erato. The loci are found in the same genomic location in these two species and interestingly, also in that of the P supergene locus of H. numata (which controls a polymorphic whole-wing patterning system) [9]. This suggests that in Heliconius the same genetic loci are involved in the repeated evolution of adaptive traits.
In other butterfly species, such as Bicyclus anynana, conserved developmental pathways appear to have been co-opted to a role in development of wing pattern elements like eyespots [10,11]. Key transcription factors are involved, such as Notch, Hedgehog and Engrailed [12][13][14], which have possibly evolved their new role through cis-regulatory changes. The developmental basis of wing colour patterning in Heliconius has yet to be elucidated. Positional cloning and sequencing of the HmYb and HmB/D loci and their orthologous loci in H. erato have revealed genes that have not been implicated previously in butterfly wing patterning [15][16][17]. Work is ongoing to further identify the switch genes within the HmYb and HmB/D regions using population genetics and gene expression approaches. Genetic changes at these switch genes among different colour pattern races are likely to involve cis-regulatory or coding sequence changes, changes to post-transcriptional control or a combination of these. MicroRNAs (miRNAs) are important post-transcriptional regulators of gene expression that have been particularly implicated in the fine-tuning of transient and complicated developmental processes. Therefore, they could have a role in the regulation and development of wing patterning that occurs during larval/pupal transitional stages in butterflies. miRNAs are 19-25 nucleotides long, endogenously expressed non-coding RNAs (for a review see [18]). In animals, mature miRNAs are derived from transcribed hairpin structures (pre-miRNAs) of around 70-90 nucleotides in length that are processed by specialised proteins (Dicer proteins). One strand of the resulting miRNA duplex is incorporated into the RNA-induced silencing complex (RISC). Post-transcriptional silencing is then mediated by binding of this complex to the 3' UTR of the target messengerRNA (mRNA), which causes degradation or translational repression of the gene. Because the 5' seed sequence (at nucleotide positions 2-8) that determines the miRNA and mRNA pairing is just seven nucleotides long in animals [19], one miRNA can have many potential targets and each mRNA can be targeted by more than one miRNA [20].
Some miRNAs are remarkably conserved across distant orders, suggesting conserved evolutionary functions. For example, over half of Caenorhabditis elegans miRNAs share sequence homology with those found in human and Drosophila [21]. One of the first miRNAs to be discovered, let-7, appears to act as an evolutionary conserved developmental timer, with temporal expression patterns being co-ordinated with progression to an adult fate [22]. Loss of let-7 function in Drosophila leads to widespread defects during metamorphosis, including small wings [23]. Several other studies have shown a key role for miRNAs in metamorphosis and many miRNAs differ in their expression patterns across life stages [24][25][26][27]. In the hemimetabolan insect Blattella germanica, prevention of miRNA processing by silencing of Dicer-1 inhibits metamorphosis, with individuals retaining nymphoid features [28]. Newly emerged miRNA genes have been detected in Drosophila and these genes seem to be evolving adaptively and sometimes rapidly [29,30]. Hence, in addition to performing conserved roles, miRNAs could also be involved in the fine-tuning of gene expression patterns underlying the evolution of novel phenotypic traits.
In this study we generated and sequenced small RNA libraries from mixed larval and pupal wings of two colour pattern races of H. melpomene; H. m. rosina, which has the yellow hindwing bar encoded by the HmYb locus and H. m. melpomene, which does not. Our aims were to characterise the first miRNAs in Heliconius butterflies, to examine differences in expression between two colour pattern races and to identify miRNAs encoded within the HmYb region and elsewhere in the genome.

miRNA characterisation in Heliconius
A total of 6,895,260 processed sequences were obtained from the small RNA libraries developed for both Heliconius races; with more reads being obtained for H. m. melpomene (3,967,516) than for H. m. rosina (2,927,744). These data have been submitted to the NCBI Gene Expression Omnibus (GEO) [31,32] under accession number GSE23292. Sequences were between 16 and 27 nucleotides in length, with 23 nucleotides being the most abundant length for both races ( Figure  1). A comparison of these sequences to miRNAs housed in miRBase (release 15) revealed 142 different sequences with homology to previously identified miRNAs. 49.06% and 45.91% of small RNA sequences were identified as miRNAs for H. m. melpomene and H. m. rosina respectively. The remaining sequences are likely to include as yet unidentified Heliconius miRNAs and other types of small regulatory RNAs. A full list of detected miRNAs, along with normalised abundances for each colour pattern race is given in Additional file 1. The sequences obtained for each race along with read counts are given in Additional file 2 (H. m. melpomene) and Additional file 3 (H. m. rosina).

Expression of miRNAs
By far the most abundant miRNA, miR-31, accounts for > 50% of all sequenced miRNAs in both colour pattern races. Altogether, 79 miRNAs were present at more than five reads per million sequenced for at least one race and 55 had more than five reads per million for both races. Considering only these 55 miRNAs, 40 had differences in abundance greater than 20% between the two races. Nine of these (chosen for showing large differences in abundance and/or high abundance) were subject to Northern blot analysis across time-staged pupal hindwing development ( Figure 2). For each race, hindwing tissue collected from two biological replicates at 24, 48 and 72 hours (± 30 minutes) post-pupation was compared to forewing and thorax tissue collected at the same time (shown in Additional file 4) and to a U6  non-coding small nuclear RNA loading control. In all cases the nine miRNAs were expressed in the three tissue types, i.e. none were either wing, or hindwing specific in their expression. miR-31, miR-10 and miR-308 showed a ubiquitous pattern of expression across stages and races. miR-276 showed an upregulation in both races at 48 and 72 hours compared to 24 and this difference seems to be more pronounced in forewings than hindwings (Additional file 4). For all other miRNAs, expression patterns varied across developmental stage and were mirrored in hindwings and forewings of both races, but not thorax. miR-184 and miR-317 were upregulated at 24 and 48 hours compared to 72; miR-263 was up-regulated at 48 and 72 hours compared to 24; miR-275 was up-regulated at 48 hours compared to 24 and 72 and miR-277 was up-regulated at 72 hours compared to 24 and 48.

Identification of miRNAs in HmYb
When short RNA sequence reads were compared to all available H. melpomene BAC sequences (currently this consists of~2.5 Mb and includes the HmYb/Sb and HmB/D regions [15,17]) two miRNAs were identified. A comparison of these sequences to miRNAs in miRBase identified them as hme-miR-193 and hme-miR-2788. miR-193 is a conserved miRNA found in mammals, birds, fish and insects. miR-2788 is a newly described miRNA, identified at this time only in B. mori [33]. Their predicted stem-loop structures are shown in Figure 3 and an alignment of sequenced reads to the BAC sequence, along with each read count are given in Additional file 5. Interestingly, they are located just 2380 bp apart in the HmYb region, in an intergenic location between genes HM00025; a putative member of the fizzy family, and HM00026; a putative homologue of poly(A)-specific ribonuclease (PARN). The butterfly hme-miR-2788 has three mismatches with its closest sequence, found in B. mori. Sequence abundances per million reads calculated allowing for these three mismatches between our sequences and those of mature miRNAs present in miRBase (and not two as shown in Additional file 1) were 8 Figure 4). There was no evidence that they were differentially expressed between colour pattern races, however both show an upregulation from 24 to 72 hours post-pupation. A further search for predicted Heliconius miRNAs within the BAC sequences did not identify any new miRNAs, although this analysis did again identify the two miRNAs described above as well as their star sequences (the 'passenger' strand of the pre-miRNA duplex that is not incorporated into the RISC).

Characterisation of butterfly miRNAs
Here we identify the first miRNAs in butterflies and characterise their expression in developing larval and pupal wings. The size distribution of sequenced reads showed peaks at 23 nucleotides for both Heliconius races. Average miRNA length is 22 nucleotides in animals [18]. Using similar deep sequencing techniques, other studies on insects show small RNA sequence size distributions with a peak at 22 nucleotides [34][35][36] and miRNA sequence sizes appear to peak at 20-22 nucleotides in B. mori [24]. In our data the peak at 23 hme-miR-193 hme-miR-2788 nucleotides is due to miR-31, which is most highly abundant as a 23mer in both races. If this miRNA is excluded, both races show peaks at 22 nucleotides, with another peak at 27 nucleotides that possibly represents Piwi-interacting RNAs; piRNAs ( Figure 1). Sequences were identified as miRNAs only if they shared homology with other miRNAs deposited in miR-Base, or if their precursor sequences were identified in Heliconius genomic DNA. This approach identified 142 butterfly miRNAs. Recently deep sequencing, or a combination of deep sequencing and bioinformatics has identified 149 miRNAs in the pea aphid, Acyrthosiphon pisum [34], and 65 and 77 miRNAs in Aedes and Culex mosquitoes respectively [35]. Two different studies of B. mori yielded 101 conserved and 14 novel miRNAs [24] and 257 miRNAs, of which 202 are described as unique to this species [33]. At this time, 171 Drosophila melanogaster miRNAs have been deposited in miRBase. All of our identified miRNAs show homology to other miR-NAs in miRBase, as a search of genomic DNA sequences failed to reveal any predicted Heliconius specific miRNA precursor sequences. However, only 2.5Mb of genomic BAC clone sequence is available at this time. The completion of the Heliconius genome-sequencing project in the near future will facilitate the identification of further miRNAs. These could include butterfly or Heliconius specific miRNAs, or miRNAs that were not detected in this study as they are not expressed in wings.

Expression of butterfly miRNAs in developing wings
The most highly abundant miRNAs, miR-31, miR-184 and miR-263 account for around 60% of detected miR-NAs in both races. All three of these miRNAs are also highly conserved and may have important roles in butterfly wing development. Remarkably, miR-31 represented 51% and 56% of all miRNAs sequenced in H. m. melpomene and H. m. rosina respectively and was expressed in all tissues and developmental stages examined. This miRNA has been best studied in humans, where its dysregulation is implicated in breast, lung, colorectal and head and neck cancers [37][38][39][40]. In Drosophila it is expressed in a pair rule pattern of 14 stripes, in the anterior endoderm and the hindgut [  depleted embryos complete development, but have severe segmentation defects including abnormal cuticle patterns and a complete loss of alternating segments [42]. In B. mori it shows highest expression during newly hatched larva and moulting larval stages [26] leading to the suggestion that it controls epithelial metabolism during moulting. The very high abundance of this miRNA in Heliconius during 5 th larval instar and pupal developmental stages suggests a role in butterfly metamorphosis.
In the Drosophila embryo, miR-184 is expressed in the central nervous system, brain and imaginal discs [43]. Expression continues into the adult stages and it is also expressed throughout the life cycle of B. mori [25] and Aedes and Culex mosquitoes [44]. It is the most highly expressed miRNA in the mosquito A. aegypti [35]. Hence this miRNA may have a general function in insects. However, it shows a decreasing pattern of expression from 24 to 72 hours post-pupation in both races of Heliconius (this pattern is not observed for thorax) indicating it may possibly have a more specialised role in butterfly wing development.
miR-263 is preferentially and abundantly expressed in B. mori pupae compared to adults [24]. In Drosophila it appears to have a conserved affiliation with sensory organ differentiation and is expressed in peripheral sense organs [45]. Drosophila deletion mutants show a sporadic loss of sensory bristles, revealing the role of this miRNA in ensuring the developmental 'robustness' of patterning processes (for example by protecting interommatidial bristles during the programmed cell death that refines the pattern of the pupal retina). Intriguingly, this miRNA increases in expression dramatically in Heliconius between 24 and 48 hours post-pupation in both forewings and hindwings. Insect sensory bristles and lepidopteron wing scales are widely believed to be homologous structures (with scales being evolutionarily derived from bristles) and share both a specification mechanism and a Notch-mediated lateral inhibition mechanism of patterning [46]. Early scale-forming cell determination occurs 16 hours after pupation, hence the timing of expression of miR-263 in Heliconius could be suggestive of a shared role in the developmental pathway of insect sensory organ bristle and lepidopteron wing scale cell patterning.
Most of the miRNAs examined in this study showed dynamic changes in expression pattern throughout pupal wing development, indicating that miRNAs in general have an important role in this process. Changes were mirrored in forewings and hindwings of both colour pattern races (with the exception of miR-276) indicating a role for these miRNAs in general wing development rather than in specific colour patterning. This was in contrast to the sequencing results, where differences in normalised read counts of sequences between races indicated that some miRNAs might be differentially regulated in wings with different colour patterns. Great care was taken to standardise the pool of RNA used for sequencing each race, such that each pool contained the same amount of RNA from each developmental stage represented. However, 5 th instar larval wings in particular are difficult to stage with a high degree of accuracy and development is both dynamic and very rapid. Therefore it seems likely that most differences in miRNA abundances observed from the sequence data are due to differences in developmental staging. Future studies of the role of miRNAs in wing colour patterning and pigmentation could investigate miRNAs that are less abundantly expressed, using more sensitive methods for analysing changes in expression (such a quantitative PCR) and also examine specific areas of developing wings (such as that destined to become a yellow bar) to look for regional-specific up-or down-regulation of miRNAs. An understanding of the dynamics of miRNA expression in adult wings would also require further work. It seems likely that different miRNAs would be involved than those expressed during wing development, however once scale cells are developed and pigmented the wing tissue may not have high levels of transcriptional activity (scale cells are not replaced if lost or removed and wings generally fade in colour over time).

miRNAs in HmYb
Two miRNAs, hme-miR-193 and hme-miR-2788, and their star sequences were discovered in the HmYb region. These are located together in an intergenic region 276bp downstream of the stop codon of PARN (HM00026) and may be co-transcribed. Indeed, their sequences are both represented at a 1.8 fold increase in abundance in H. m. melpomene. The other flanking gene (Fizzy; HM00025) shows differences in splicing that associate with wing colour pattern and is a possible candidate gene for the HmYb switch locus. In contrast to the nine other miRNAs examined by Northern blot (described above) these miRNAs did not appear to be expressed in thorax tissue. This implies a wing-specific role. The Northern blots show no evidence for differential expression of the two miRNAs between colour pattern races. Furthermore, both are expressed in forewings and hindwings indicating that their expression is not specific to the yellow hindwing bar pattern encoded by HmYb. However, they are strongly up-regulated between 24 and 72 hours in wings, when patterning and pigmentation processes are specified. A role for these two miR-NAs in wing colour patterning will be the subject of future research.

Conclusions
Here we present the first butterfly miRNAs and investigate their role in wing development. Several conserved miRNAs are highly expressed in Heliconius pupal wings and one, miR-263, could be involved in wing scale cell patterning in a developmental mechanism that is shared with insect sensory organ bristle patterning. Other miR-NAs show dynamic changes in expression in developing pupal wings, indicating an important role for miRNAs in this process. This study did not detect changes in expression of miRNAs associated with wing colour pattern, however two miRNAs were detected in the HmYb region that contains colour pattern switch loci. These miRNAs are located close together near a gene that associates with colour pattern and are up-regulated in wings between 24 and 72 hours post-pupation. Future work will investigate possible roles for these and other miRNAs in wing colour patterning.

Sample collection
Wild H. m. rosina were collected in Gamboa, Panama (09°07'N/79°42'W) and H. m. melpomene were collected in Darien, Panama (09°10'N/78°43'W). Individuals for the study were reared in insectaries at the Smithsonian Tropical Research Institute in Gamboa, Panama. Females were allowed to freely oviposit on Passiflora menispermifolia. Eggs were collected daily and larvae were raised in individual pots supplied with fresh young leaves of P. biflora. For library development, whole hindwings and forewings were dissected from 5 th instar larvae and visually staged for development according to standardised laboratory criteria and those outlined in [47] and [48]. For pupae, wings were collected from two developmental stages, early pupae (where no pigments are visible on wings) and mid-melanin (where ommochrome and melanin pigments are visible). Tissues were stored in RNAlater ® (Ambion). For Northern blot analyses forewing, hindwing and thorax tissue samples were collected separately from a total of nine individuals for each of the two races (three individual replicates for each of three time points post-pupation). Dissections were performed at 24, 48 and 72 hours post-pupation (± 30 minutes) using a Leica Stereozoom 4 microscope under RNase free conditions. All tissues were placed in RNAlater ® (Ambion) and stored at -80°C degrees until extraction.

RNA extraction, miRNA sequencing and Northern blot
Total RNA (including small RNAs) was extracted using the mirVana™ kit (Ambion) according to the manufacturer's recommendations. RNA was quantified and checked for purity and integrity using the Agilent 2100 Bioanalyzer. For each race, 100 μg of total RNA from 11 individuals was pooled such that each pool contained 50% larval and 50% pupal RNA in equal amounts for each developmental stage (the detailed composition for each pool was 4.1% larval stage < 1; 2% larval stage 1-1.75; 2.9% larval stage 2-2.5; 22% larval stage 2.75-3; 19% larval stage > 3; 25% early pupae; 25% mid-melanin pupae). Both larval and pupal wing stages were included as the most dynamic changes in wing development begin occurring during the final (5th) instar and continue throughout pupation. Small RNA fractions of between 19-24 nucleotides were then isolated and used for short RNA library generation as described in [49]. The libraries were sequenced using an Illumina Genome Analyzer II at BaseClear (Leiden, The Netherlands). 2.5 μg of total RNA was used for Northern blot analysis as described previously [50]. Briefly, RNA from two biological replicates for each race/developmental stage was separated on a 15% denaturing polyacrylamide gel and blotted to Hybond NX membranes (Amersham). Expression of small RNAs was assessed by hybridisation to a [P 32 ]-labelled (Perkin Elmer, UK) nucleic acid oligonucleotide probe and compared to a U6 loading control. Probe sequences are given in Additional file 6.

Computational analysis
Raw Illumina reads were processed by first converting FASTQ to FASTA format, and then removing any adaptor sequences with exact matches to the first eight bases of the 3' adaptor. Any sequences without adaptor matches were excluded from further analyses. Normalised miRNA expression levels were compared between the two Heliconius races using miRProf [51] allowing up to two mismatches to mature miRNA sequences present in miRBase [52]. miRNA predictions were performed by running miRCat [51]  The predicted miRNAs miR-193 and miR-2788 (and their corresponding sequences and star sequences) were submitted to miRBase [52] and named hme-miR-193 and hme-miR-2788 respectively.

Additional material
Additional file 1: Heliconius miRNAs identified by deep sequencing of two colour pattern races. miRNAs identified in H. m. melpomene and H. m. rosina. Two or more miRNAs are listed when a sequence matches to two different, but related miRNAs in miRBase (the sequences and read counts of all miRNA variants detected are given in Additional files 2 and 3). Normalised abundance (number of reads per million; calculated allowing for two mismatches to mature miRNA sequences given in miRBase) is given for each race along with fold change (melpomene/ rosina).
Additional file 2: miRNA sequences and read counts detected in H. m. melpomene. FASTA file of all identified miRNA sequence variants detected in H. m. melpomene. x refers to the number of reads obtained for each sequence. Additional file 6: Probe sequences used for Northern blots. Probe sequences for ten miRNAs analysed by Northern blot.