Identification and differential expression of piRNAs in the gonads of Amur sturgeon (Acipenser schrenckii)

Ethics statement

The protocol was approved by the Committee on the Ethics of Animal Experiments of Guangdong Institute of Applied Biological Resources (GIABR2014008). Individual sturgeon were immersed in water with 10-4 (v/v) eugenol for approximately 1–3 min for euthanasia, according to the AVMA guidelines (Leary et al., 2013). All efforts were made to minimize suffering.

Sample and RNA preparation

In this study, we used 3-year-old Amur sturgeons (Acipenser schrenckii), whose sex could be identified accurately by laparoscopy and histology. The animals were obtained from the Engineering and Technology Center of Sturgeon Breeding and Cultivation of the Chinese Academy of Fishery Science (Beijing, China). The testes and ovaries of six 3-year-old Amur sturgeons (three males and three females) were collected. Total RNA was extracted from tissue samples separately with RNAiso reagent (TaKaRa, Shiga, Japan) according to the manufacturer’s instructions. The RNA concentrations were measured using a Qubit RNA Assay Kit in Qubit 2.0 Fluorometer (Life Technologies, Carslbad, CA, USA), and RNA purity was assessed using a Nano Photometer spectrophotometer (IMPLEN, Westlake Village, CA, USA). RNA integrity was inspected using an RNA Nano 6000 Assay Kit and Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA).

Small RNA library preparation and sequencing

Four RNA samples (two testes and two ovaries, 3 µg RNA of each) were used for the construction and sequencing of small RNA libraries. In brief, NEB 3′ SR Adaptor ligation, SR RT Primer hybridization and NEB 5′ SR Adaptor ligation were performed according to the NEBNext Multiplex Small RNA Sample Preparation Set (Illumina, San Diego, CA, USA) protocol. After first strand cDNA synthesis using M-MuLVA Reverse Transcriptase (RNase H⁻) by PCR, the 140 bp–160 bp products (with adaptors on both sides) were separated on an 8% PAGE gel and quantified using an Agilent Bioanalyzer 2100 system. Then, a cluster of index-coded samples was generated using a TruSeq SR Cluster Kit v3-cBot-HS (Illumina, San Diego, CA, USA) and sequenced on an Illumina HiSeq 2000 platform. Finally, 50 bp single-end reads were generated.

Small RNA annotation and piRNA identification

In this study, the testis and ovary transcriptomes of A. Schrenckii (Jin et al., 2015) were used as reference sequences for small RNA annotation. After removal of unclean reads (adapters, low quality reads, reads containing ‘n’, and redundant reads), clean unique reads were mapped onto the A. schrenckii transcriptome reference sequences using the program Bowtie (Langmead et al., 2009) with no mismatches allowed. Perfectly mapped reads were filtered by three successive steps to remove small RNA elements: (1) searched against Sanger miRBase (Release 19) to exclude conserved miRNAs; (2) screened against Rfam (http://rfam.sanger.ac.uk/) and RepeatMasker (http://www.repeatmasker.org/) with Bowtie to filter the sequences originating from rRNA, tRNA, snRNA, snoRNA and repeats; (3) analyzed by miREvo, mirdeep2 and MirCheck to remove the potential novel miRNA reads. The detailed process of small RNA filter was described in our previous study (Yuan et al., 2014).

The remaining reads with lengths of 26–32 nt were used for piRNA discovery. A k-mer scheme relied on the training sets of non-piRNA and the piRNA sequences of five model species (rat, mouse, human, fruit fly and nematode), was applied as previously described (Zhang, Wang & Kang, 2011). Then, it was compared with the existing ‘static’ scheme on the basis of the position-specific base usage. Putative novel piRNAs were scanned against piRNABank, (http://pirnabank.ibab.ac.in) using Bowtie with no mismatches allowed to identify the orthologs of known piRNAs. The relative frequencies of piRNA nucleotide utility reads were analyzed. Subsequently, the putative piRNAs were functionally annotated using the coding sequence (CDS), 5′ and 3′ untranslated region (UTR) database of the A. Schrenckii transcriptome with TransDecoder software (https://github.com/trinityrnaseq/trinityrnaseq/wiki/Coding-Region-Identification-in-Trinity-Assemblies).

piRNA microarray and data analysis

To validate the expression of putative piRNAs identified by Illumina sequencing, we selected: (1) 1,092 that exhibited significantly different expression in ovaries vs. testes with |logFC| ≥ 1.5, P ≤ 0.01, Padj ≤ 0.05, and read counts >10; and (2) 779 highly expressed putative piRNAs with read counts 50∼4,646 in ovaries or 50∼8,177 in testes. For each of the six RNA samples (three ovaries and three testes, extracted above) assessed, 4 µg of total RNA were used to hybridize with the microarray chip.

A piRNA microarray was manufactured by RIBOBIO (Guangdong, China), and each piRNA probe had three replicates. The chip was hybridized with single-color labeling (Cy5) RNAs according to the manufacturer’s protocol, with no modifications. Microarray results were extracted using a laser scanner (GenePix 4000B, Molecular Device) and digitized using Array-Pro image analysis software (Media Cybernetics, Rockville, MD, USA). Raw data were subtracted using the background matrix, and spots with CV [(standard deviation)/(signal intensity)] < 0.3 were normalized using a quantile normalization method to remove system related variations, including sample amount variations and signal gain differences of the scanners and to faithfully reveal the biological variations (Bullard et al., 2010). The medians of repeated data (normalized intensity) were used for statistical analyses with One-Way ANOVA method, and the expression of piRNAs was deemed significant with the criterion |logFC| ≥ 1 and P ≤ 0.05.

piRNA-generating gene prediction and annotation

Putative piRNAs were mapped onto the A. schrenckii transcriptome using the program Bowtie with no mismatches allowed to identify piRNA-generating genes. The piRNA-generating genes were mapped to Gene Ontology (GO) database (http://www.geneontology.org/) by Interproscan (Zdobnov & Apweiler, 2001) and the KEGG pathway database (http://www.genome.jp/kegg/) by BLASTX at E values 1e−10 (Kanehisa et al., 2006), respectively by two levels, (1) total of piRNA-generating genes and (2) sex-specific expression (only expressed in testes or ovaries and with a read count ≥ 10) or sex-biased expression of piRNA-generating genes in ovaries and testes (|logFC| ≥ 1 and P ≤ 0.05 obtained by microarray). Finally, the enriched functional groups or pathways were obtained with corrected P < 0.05.

Sequencing and statistics of small RNA reads

A total of 7.7 × 10⁶–8.6 × 10⁶ reads were sequenced from four small RNA libraries, with an error rate of 0.01%, Q30 > 97.3% and GC content of ∼48% (Table S1). Then, approximately 7.3 × 10⁶–8.1 × 10⁶ high-quality small RNA reads (>94% in each library) were obtained after removal of ambiguous reads (Table S1). The size distribution and frequency percentage of small RNA reads are shown in Fig. S1 and, of these, the potential piRNA reads (26–30 nt) were the major component (approximately 58%). About 3.1–3.4 × 10⁶ small RNA reads mapped to the A. schrenckii transcriptome reference sequences (Jin et al., 2015) were with a perfect match (Table S1). A total of 3.5–6.7 × 10⁴ reads were mapped to at least one putative Amur sturgeon miRNA precursor by searching against miRBase, corresponding to 1.0–2.2% perfectly matched small RNA reads (Table 1). Subsequent small RNA filtering indicated that other non-coding RNAs (rRNA, tRNA, snRNA and snoRNA) and repeat sequences were approximately 1.0–1.5% and 0.3–0.4%, respectively. Then, 5.3–6.6 × 10³ (0.2%) potential novel miRNA reads specific to A. schrenckii were detected with miREvo (Table 1).

piRNA discovery

According to a k-mer scheme analysis which relied on the training sets of non-piRNA and the piRNA sequences of five model species (Zhang, Wang & Kang, 2011), we identified total of 5.7–6.2 × 10⁵ piRNA reads in our sequencing results, that account for 16.6–19.8% of perfectly matched small RNA reads (Table 1). By following steps described in the methods section, 875,679 putative piRNAs from the testis and ovary libraries were predicted, including 93 piRNAs that were homologous to known piRNAs (four in fruit fly, 45 in zebra fish, 34 in human, four in mouse and six in rat, see Table S2). Sturgeon putative piRNAs have a strong U bias at the extreme 5′ position and are enriched for https://www.sciencedirect.com/topics/biochemistry-genetics-and-molecular-biology/adenine (A) at position 10 (Fig. 1A). Further analysis indicated that their rate of sex-specific expression (only expressed in ovary or testis) was up to 87% (767,805/875,679), and the ratio of ovary- vs. testis-specific putative piRNAs was nearly 1:1 (43.99%,385,222/875,679 vs. 43.69%, 382,583/875,679, see Fig. 1B). Moreover, we found a similar percentage of sense and antisense strand putative piRNAs (42.99%, 376,487/875,679 vs. 57.01%, 499,192/875,679, respectively). In total, approximately 42% (175,648 in ovaries and 195,552 in testes) of all putative piRNAs were derived from well-annotated gene regions, 58% were from unannotated genomic regions (Fig. 1C). Of these, only 31,949 (3.65%) and 73,802 (8.43%) matched the 5′ UTR and 3′ UTR, respectively. Abundant piRNAs in gene regions were found in CDS, which represented 265,499 (30.32%) putative piRNAs.

Sturgeon putative piRNAs are present in both male and female gonads

Small RNA reads of both male and female sturgeon gonads displayed a peak in reads of 26–30 nt in length (55.8% in testes and 56.5% in ovaries, shown in Fig. 2A). This is an expected length distribution for potential novel piRNA reads. Whereas Small RNA reads in five other somatic/non-germline tissues (pool of brain, heart, muscle, liver and spleen tissues of A. schrenckii (Yuan et al., 2014)) showed a peak in reads of 20–24 nt in length which is an expected length distribution for miRNA reads (∼52.9.4%). Moreover, both ovary and testis samples had a strong 5′ U bias in reads of 20–24 nt (miRNA reads) and 26–30 nt (piRNAs reads) in length, whereas only reads of 20–24 nt among somatic samples had a 5′ U bias (Figs. 2B–2D). The library size of testis and ovary is nearly equal (the ratio of total read counts ≈1:1, Table S1), thus we further compared piRNAs vs miRNAs in sturgeon testis and ovary samples, and obtained a similar ratio of 4.7:1∼5.5:1 (Table 1). This clearly indicates putative piRNAs are far more abundant class of small RNAs in both testis and ovaries compared to miRNAs. Moreover, both piRNA (1:1) and miRNA (1.2:1) showed an equal ratio in testes vs. ovaries (Table 1).

piRNA-generating genes

A total of 875,679 putative piRNAs were mapped onto the A. schrenckii transcriptome (including 122,381 unigenes in testis and 114,527 in ovaries, Jin et al., 2015). Total 49,390 piRNA-generating genes were obtained (23,079 in ovaries, 26,311 in testes, and 9634 in both). The mean length of piRNA-generated genes was 967 nt (median 506), with a range of 150–16,256 nt. There was a mean of 18 unique piRNAs reads per gene (scale of 1 up to 4,181 piRNAs per gene).

To illustrate the functions of the sturgeon putative piRNAs, gene function of the piRNA-generated genes was annotated by searching against GO database with Interproscan and KEGG pathway database by BLASTX at E values 1e−10. The GO analysis showed 264 enriched GO terms for biological process, molecular function and cellular component, in them, 13 GO terms related to reproduction (such as sexual reproduction, mating/reproductive behavior, post-mating regulation of female receptivity) were obtained (Table S3). Moreover, 29 enriched KEGG pathways, including six involved in aquaculture, such as pathogen infection, multiple diseases, and antigen processing, were identified (Table S3). We further investigated the function of piRNA-generated genes which with sex-specific piRNAs located (1,078 only detected in testis and 1,237 only detected in ovary, with read counts ≥ 10), and more GO terms related to reproductive processes were identified (Fig. 3 and Table S4).

Validation of putative piRNA expression by microarray

We used an independent microarray platform to validate the expression of 1,871 putative piRNAs (1,092 differently expressed and 779 highly expressed piRNAs, denoted ASY-piRNA-X, where Xs are numerals) identified by Illumina sequencing. A total of 311 showed significantly biased expression between sexes (deemed sex-biased expression) and were clustered into ten clades including five clades up-regulated in ovaries and five clades up-regulated in testes (Fig. S2 and Table S5). Of these, 124 were specifically up-regulated in ovaries with 1 < |logFC| < 8.5 and P ≤ 0.05, and 187 were specifically up-regulated in testes with 1 < |logFC| < 4.1 and P ≤ 0.05 (Table S5). In sturgeon testes, the putative piRNA with the highest expression level was ASY-piR-1255 with median of normalized intensity 4.04 in testes vs. 0.72 in ovaries (Table S5). In contrast, the most expressed in the ovaries was ASY-piR-122 with median of normalized intensity 13.33 in ovaries vs. 4.91 in testes (Table S5). Moreover, 1,335 were co-expressed in both the testes and ovaries (median >1 in both and P > 0.05). Further analysis indicated no significant difference in size distribution among sex-biased piRNAs (average 27.6 nt in testes and 27.8 nt in ovaries) and co-expressed piRNAs in both gonads (average 27.5 nt). The Illumina sequencing and microarray expression of putative piRNAs was especially correlated for 311 sex-biased piRNAs (r = 0.489, P = 0.000). GO and KEGG annotation indicated that the 311 sex-biased expressed piRNAs mainly participate in the metabolic processes of sturgeon gonads.