RNA sequencing reveals small RNAs differentially expressed between incipient Japanese threespine sticklebacks

miRNA transcriptome

We conducted small RNA sequencing of four males and four females from both the Pacific Ocean and the Japan Sea forms using the Illumina Genome Analyzer IIx. After quality control of the sequence reads, data was collected from 26.2 ± 3.3 (mean ± SD) million reads per fish (Additional file 1: Table S1). Most of these reads (23.7 ± 2.9 million reads; 86.9–94.6% of the total reads) were located in the Ensembl stickleback miRNA database. More than 50% of these small RNAs were 22 nucleotides in length (Figure 1). In total, 1300 isoforms of miRNA were detected in the brain (924 in the Pacific Ocean males, 924 in the Pacific Ocean females, 916 in the Japan Sea males, and 884 in the Japan Sea females). miRNA expression profiles demonstrated that miRNAs homologous to mir21, mir100, let7, mir101, mir143, and mir9 were the most abundant in the stickleback brain, regardless of the species or sex (Figure 2). Most other miRNAs were expressed at relatively low levels (less than 3% of the annotated reads).

To elucidate the variation in the miRNA transcriptomes, we conducted principal component analysis (Figure 3; Additional file 2: Table S2). The miRNA transcriptomes of the Pacific Ocean and the Japan Sea forms make distinct clusters. Interestingly, the Japan Sea males and females make distinct clusters, whereas the miRNA transcriptomes of the Pacific Ocean males and females overlapped. These data suggest that the magnitude of sex differences of miRNA expression levels might differ between species.

Sex differences in the expression levels were identified for several miRNAs (Tables 2 and 3). Interestingly, all sex-biased miRNAs belonged to the let-7 family (Table 3) [69]. miRNAs of the let-7 family are highly conserved across taxa and are important during development [70, 71]. Two sexually dimorphic miRNAs (ENSGACT00000028241 and ENSGACT00000029064) were expressed from LG19; both of these were female-biased. In species with an XY-sex chromosome system, suppression of recombination can lead to degeneration of the Y chromosome [72, 73]. Unless dosage compensation mechanisms evolve, expression of genes located on the X-specific region becomes female-biased [74]. Therefore, we investigated the relationship between Y-degeneration and sex differences in miRNA expression levels. We found that expression of all miRNAs derived from the X-specific region (i.e. the counterpart region of the Y is likely degenerated) were female-biased, whereas expression of miRNA derived from the pseudoautosomal region was not necessarily female-biased (Figure 4), and log₂ (fold difference between male and female) significantly differed between the pseudo-autosomal and the X-specific regions (Mann–Whitney U test, U = 7, Z = 814.5, P < 0.001, N = 17). These results suggest that Y chromosome degeneration may have a substantial impact not only on mRNA expression [59], but also on miRNA expression.

Transcript ID	LG	Start position	Median reads per million (RPM) ± S.D.				Zebrafish homolog
			Pacific Ocean males	Pacific Ocean females	Japan Sea males	Japan Sea females
ENSGACT00000029075	IV	20044816	1606.0 (181.4)	2616.6 (434.6)	1421.9 (60.4)	2339.2 (396.8)	mirlet7i
ENSGACT00000028051	IV	19861309	20839.3 (1126.5)	22453.2 (1081.8)	18010.6 (611.4)	25558.9 (1427.5)	mirlet7a-3
ENSGACT00000028218	VII	18314432	19933.6 (1139.6)	21367.4 (1035.2)	17101.4 (549.5)	24133.6 (1426.5)	mirlet7a-4
ENSGACT00000028988	XII	12165589	20768.4 (1137.2)	22362.8 (1085.0)	17924.8 (600.1)	25445.9 (1434.9)	mirlet7a-6
ENSGACT00000028001	XVII	11962085	20520.0 (1120.3)	21925.0 (1071.1)	17736.6 (610.0)	24909.7 (1419.4)	mirlet7a-1
ENSGACT00000029000	XVII	3393574	20777.3 (1135.5)	22362.4 (1084.3)	17935.0 (602.4)	25456.6 (1431.4)	mirlet7a-6
ENSGACT00000029064	XIX	4582643	1608.3 (184.2)	2633.7 (436.6)	1411.8 (66.8)	2327.6 (395.7)	mirlet7i
ENSGACT00000028241	XIX	4785031	20810.0 (1125.6)	22433.4 (1083.4)	17991.0 (609.0)	25531.1 (1426.0)	mirlet7a-3

LG: linkage group.

Because miRNAs can target many mRNAs [75, 76], divergence in miRNA expression patterns may have substantial effects on the expression patterns of many mRNAs. Further experimental studies examining the roles of small RNAs in fish will be necessary to understand the functional effects of the miRNA transcriptome variation. This would be possible by developing either transgenic fish specifically overexpressing small RNAs or small RNA-deficient knockout fish [77–81].

Small RNAs homologous to repetitive sequences

Small RNAs with no matches in the stickleback non-coding RNA database were further analyzed. A histogram of the read length of these small RNAs revealed two peaks, with one peak at 22 nucleotides and the other at 27–29 nucleotides (Figure 5). The fraction of large-sized small RNAs may contain piRNAs, whereas the other fraction may correspond to novel miRNAs and/or endo-siRNAs. A histogram of reads per million (RPM) of the unidentified small RNAs revealed that most were expressed at low levels, with only a few expressed at high levels (Figure 6). Thirty-one novel small RNAs expressed at relatively high abundance (mean RPM > 50 in at least one of the four groups) could be classified into 17 isoforms on the basis of sequence identity (Additional file 3: Table S3).

The remaining six and three isoforms overlapped with repetitive sequences homologous to tRNA and rRNA, respectively. Previous studies also identified a number of tRNA-derived small RNAs in humans [83, 84], Giardia lamblia[85], and zebrafish [86]. These tRNA-derived small RNAs may contribute to gene regulation [84], although little is yet known about their functions. Interestingly, some transposons are derived from tRNAs [87, 88], so there appears to be an interesting link between tRNA and small RNAs.

Finally, we found that the iso-smRNA9, whose potential targets are predicted to be non-LTR transposons, was more highly expressed in the Japan Sea sticklebacks than in the Pacific Ocean sticklebacks (ANOVA, F_{1, 13} = 14.5, P = 0.002). Although expression levels of some other piRNAs may differ between different species and sexes, the differences were not significant after Bonferroni correction (Additional file 5: Table S5). None of the non-annotated small RNAs showed sex differences in the expression levels (Additional file 5: Table S5).

Thus, we identified small RNAs homologous to repetitive sequences such as RNA and DNA transposons. Hybrids between species often exhibit courtship dysfunction [89–94]. Abnormal transposon activity in hybrids may cause hybrid courtship dysfunction, but this has not been tested in any organism. Intra- and inter-population variation in the presence and absence of non-LTR retrotransposons has been found in sticklebacks [95]. In addition, our analysis involving whole genome sequence comparisons also revealed that DNA transposon insertion sites diverge between the Pacific Ocean and the Japan Sea sticklebacks (Kitano, unpublished data). Therefore, further studies examining variation in transposon activity between the different species and transposon activity in the hybrids will lead to a better understanding of speciation mechanisms.

Small RNA sequencing

Sympatric Pacific Ocean and Japan Sea sticklebacks were collected using minnow traps from the Bekanbeushi River System on Hokkaido Island, Japan in June 2007 [60, 62]. Fish were brought to a laboratory to examine the brain transcriptome of courting male and spawning female fish, and mating experiments were conducted in June and July 2007, as described previously [60, 63, 96]. Once the male fish had constructed a nest, a conspecific gravid female fish was placed in the same tank. Immediately after the female fish had inspected the nest, both the male and the female fish were removed from the tank prior to spawning. After immersing the fish in a lethal dose of tricaine methanesulfonate, the brains of each fish were dissected and stored separately at −70°C.

For RNA sequencing, we used four Pacific Ocean male, four Pacific Ocean female, four Japan Sea male, and four Japan Sea female fish (N = 16 in total). Total RNA was isolated using TRIzol Reagent (Life Technologies, Grand Island, NY, USA), and the quality of RNA was evaluated using the BioAnalyzer (Agilent, Santa Clara, CA, USA). RNA Integrity Number (RIN) ranged from 9.3 to 10 with the median of 10. Libraries were constructed using the TruSeq Small RNA Sample Preparation Kit (Illumina, San Diego, CA, USA). Small RNAs with 20–30 nucleotides were isolated according to the manufacturer’s instructions. Sequencing was performed using the Genome Analyzer IIx (Illumina, San Diego, CA, USA) at the University of Tokyo. We used 16 lanes of the Genome Analyzer IIx (one fish per one lane).

miRNA analysis

Sequence analyses were conducted using the CLC Genomics Workbench Software (CLC bio, Katrinebjerg, Denmark). First, we discarded reads with a low quality score (quality score on the Phred scale of less than 0.05), very short length (less than 14 bp), or three or more ambiguous nucleotides (Additional file 1: Table S1). Next, identical reads were clustered together to group different types of small RNAs. Next, the sequences were mapped against the Ensembl stickleback miRNA database (http://asia.ensembl.org/info/data/ftp/index.html) with two nucleotides mismatches allowed.

Principal component analysis (PCA) of reads per million (RPM) was conducted on a Pearson correlation matrix. To identify differentially expressed miRNAs between various species and sexes, statistical analyses were conducted on the square-root transformed RPM of each miRNA. Using the statistical package R [97], analysis of variance (ANOVA) was conducted to examine whether the species, sex, and their interactions significantly influenced RPM. Because the patterns of sex differences varied between species (see Figure 3), we included the interaction term for the analysis. The Bonferroni correction was used for multiple comparison correction. Shapiro-Wilk and Bartlett’s tests were used to test for normal distribution of the data and homogeneity of variances, respectively. None of the miRNAs violated the assumptions of the normal distribution or homogeneity of variances after Bonferroni correction. Only miRNAs with a mean RPM higher than 1000 in at least one of the four groups (Pacific Ocean male, Pacific Ocean female, Japan Sea male, or Japan Sea female) were used for PCA and ANOVA.

piRNA analysis

Small RNAs with no match to any sequences in the stickleback miRNA database were mapped against sequences in the Ensembl stickleback non-coding RNA database (http://asia.ensembl.org/info/data/ftp/index.html), which includes transfer RNA (tRNA), ribosomal RNA (rRNA), small cytoplasmic RNA, small nuclear RNA, small nucleolar RNA, microRNA precursors, long intergenic non-coding RNAs, and other miscellaneous RNA. The parameters used were the same as above. Small RNAs with no match to any sequences in the stickleback non-coding RNA database may contain piRNAs; thus, we further analyzed these small RNAs to identify stickleback piRNAs. We first examined non-annotated small RNAs with nucleotide lengths >25 because this fraction is more likely to contain piRNAs. Among these longer small RNAs, we examined small RNAs with a mean RPM higher than 50 in at least one of the four groups (Pacific Ocean male, Pacific Ocean female, Japan Sea male, or Japan Sea female). Although some small RNAs showed differences in length, identical sequences were observed (Additional file 3: Table S3); hence, identical reads of different sizes were considered to be the same isoform.

To examine whether these small RNAs contain homologous sequences against any repetitive sequences, the longest isoforms were blasted against the stickleback genome (BROADS 1.56) using the default parameters (match = 1; mismatch = −3; gap existence = 5; gap extension = 2) of the CLC Genomics Workbench software. Next, flanking sequences (3,000 bp of upstream and 3,000 bp of downstream) were downloaded using Perl script [98]. We then examined whether these regions contained any repetitive sequences using the CENSOR software [99] on the Genetic Information Research Institute website (http://www.girinst.org/). When hits of repetitive sequences were identified in the region, we investigated whether the small RNA sequences overlapped with repetitive sequences. We also blasted these sequences against the piRNABank database [100].

To identify differentially expressed piRNAs between different species and sexes, statistical analysis was conducted on the square-root transformed RPM. Square-root transformed RPM values were subjected to ANOVA, followed by Bonferroni correction. Shapiro-Wilk and Bartlett’s tests were used to test for normal distribution of the data and homogeneity of variances, respectively. Two small RNAs, iso_smRNA3 and iso_smRNA14, did not meet the assumptions of homogeneity of variance (Bartlett’s test, K-squared = 15.2, P = 0.0017) and normal distribution (Shapiro-Wilk test, W = 0.779, P = 0.0014), respectively. For these small RNAs, we also conducted a Mann–Whitney U-test to confirm that these two small RNAs did not exhibit any differences between sexes and species. Because the interaction between the sexes and different species showed no significant effect for any of these small RNAs, the interaction term was excluded.