Abstract Drosophila Blanks is a testes-specific RNA-binding protein required for post-meiotic spermiogenesis. However, Blanks's role in regulating RNA populations in the testes remains unknown. We performed small RNA and mRNA high-throughput sequencing in blanks mutant testes and controls. We identified two miRNAs, one siRNA, and hundreds of mRNAs that are significantly upregulated or downregulated in blanks mutant testes. Pathway analysis revealed that differentially expressed mRNAs are involved in catabolic and metabolic processes, anion and cation transport, mating, and reproductive behavior. Our results reveal that Blanks plays important roles in defining testicular small RNA and mRNA profiles. Keywords: Genetics, Molecular biology 1. Introduction RNA binding proteins (RBPs) are involved in every aspect of RNA biology. While most ubiquitously expressed RBPs are involved in essential cellular processes, tissue-specific RBPs play crucial roles in regulating developmental processes by shaping tissue-specific RNA profiles. Of the tissue-specific RBPs characterized thus far, the majority are expressed in the testes [27][1]. These testes-specific RBPs regulate the testicular RNA landscape, which includes multiple classes of small RNAs (microRNAs [miRNAs], short-interfering RNAs [siRNAs], and Piwi-interacting RNAs [piRNAs]) and mRNAs (Aravin et al., 2003; Aravin et al., 2006). Testes-specific RBPs are required for spermatogenic processes, including stem cell maintenance, cell division, and differentiation [[28]2, [29]3, [30]4, [31]5]. Blanks is an RNA-binding protein with two dsRNA-binding domains (dsRBDs) expressed in Drosophila testes ([32]Fig. 1) [[33]6,[34]7]. Characterization of blanks mutant flies demonstrated that Blanks plays important roles in post-meiotic spermiogenesis and male fertility [[35]6, [36]7]. Blanks^KG00804 homozygous mutant males, but not blanks^KG00804/TM3 heterozygous mutant control males, lack mature sperm in their seminal vesicles due to a defect in late stage, post-meiotic spermiogenesis, resulting in male-sterility [[37]6, [38]7]. A Blanks transgene rescued male fertility in blanks^KG00804 homozygous mutants [39][7]. However, Blanks transgenes containing missense point mutations within the first dsRBD failed to rescue male fertility in blanks^KG00804 homozygous mutants. These results suggest that Blanks RNA-binding activity is required for Blanks function in spermiogenesis. Fig. 1. [40]Fig. 1 [41]Open in a new tab Domain structure of Drosophila Blanks. Despite its biological importance in spermatogenesis, the role of Blanks in regulating RNAs in the testes remains unclear. A previous study examined levels of two miRNAs, bantam and miR-7, in blanks mutant flies and found no change in their levels compared with control flies [42][7]. However, whether Blanks regulates levels of other small RNAs, including other miRNAs, siRNAs, and piRNAs, is not known. Another study examined levels of mRNAs using a microarray-based approach, which encompassed only a fraction of the transcriptome and was limited by its reliance on probe-hybridization [43][6]. Therefore, it remains unknown whether Blanks plays roles in shaping testicular small RNA and mRNA profiles. To address this question, we used high-throughput sequencing to examine the small RNA profile and mRNA profile in blanks mutant testes in an unbiased, genome-wide approach. We sequenced RNA isolated from blanks homozygous (blanks^KG00804) and heterozygous (blanks^KG00804/TM3) mutant testes. We identified small RNAs (two miRNAs and an endogenous siRNA) that are downregulated in blanks homozygous mutant testes. We also found hundreds of mRNAs whose expression levels are upregulated or downregulated in blanks homozygous mutant testes. GO term/KEGG pathway analysis using differentially expressed genes identified biological processes, cellular compartments, molecular functions, and pathways that are likely dysregulated in blanks homozygous mutants. 2. Materials and methods 2.1. RNA preparation Testes were dissected from blanks^KG00804 homozygous and blanks^KG00804/TM3 heterozygous adult male flies. Total RNA from testes was purified using miRVana (ThermoFisher). Tens of testes were used for each biological replicate preparation, resulting in tens of microgram of purified total testes RNAs. The blanks^KG00804 mutant allele has a SUPor-P element inserted 29 nt downstream of the blanks transcriptional start site [[44]6, [45]7]. Blanks protein level in blanks^KG00804 homozygous mutant testes is reduced >400-fold compared with blanks^KG00804/TM3 heterozygous control [[46]6, [47]7]. 2.2. High-throughput sequencing of small RNAs (small RNA-seq) Three biological replicates each of small RNA libraries were prepared using size-selected 18–30 nt long RNAs by gel purification, sequenced at the Beijing Genomics Institute using HiSeq4000 platform (Illumina), and analyzed as previously described [[48]8, [49]9, [50]10, [51]11, [52]12, [53]13, [54]14, [55]15, [56]16]. Approximately 27–57 million reads were obtained for each library. Approximately 80–87% of reads were mapped to the dm6 Drosophila genome. Approximately 27–55% of non-rRNA-mapping reads were mapped to miRNA hairpins. The abundance of miRNAs, siRNAs, and piRNAs normalized by the sequencing depth (non-rRNA-mapping reads) in each library was calculated. Then their mean abundance among three biological replicates was calculated. To eliminate miRNAs with very low expression levels, which are unlikely to have a physiological role, only miRNAs (n = 94) whose mean abundance was more than 100 reads per million total reads in either blanks heterozygous or homozygous mutants as well as three endo-siRNAs (esi-1.1, esi-1.2, and esi-2.1) were analyzed in [57]Fig. 2. Fig. 2. [58]Fig. 2 [59]Open in a new tab Small RNA-seq of blanks mutant testes. (A) Scatter plot of normalized miRNA and endo-siRNA reads (reads per million non-rRNA-mapping reads) in blanks homozygous (blanks^KG00804) and heterozygous (blanks^KG00804/TM3) mutant testes determined by high-throughput sequencing. Each point represents a unique miRNA or endo-siRNA. Green points represent miRNAs (miR-14-3p and miR-317-5p) whose levels are significantly changed in blanks homozygous mutants, while black points represent miRNAs whose levels are not significantly changed. An orange point represents endo-siRNA (esi-2.1) whose level is significantly changed while gray points represent endo-siRNAs (esi-1.1 and esi-1.2) whose level are not significantly changed. Mean ± SD for three biological replicates. (B) Relative abundance of significantly dysregulated miRNAs (miR-14-3p and miR-317-5p) and an endo-siRNA (esi-2.1). Mean ± SD for three biological replicates. Adjusted P-value <0.05 is indicated by * (Student's t-test). 2.3. High-throughput sequencing of mRNAs (mRNA-seq) Poly-A+ mRNA-seq libraries were constructed as previously described [[60]8, [61]16, [62]17]. Paired-end 100 nt sequencing (2 × 100 bp) was performed by the Beijing Genomics Institute using HiSeq4000 platform (Illumina). Approximately 93–95% of paired reads were mapped to the dm3 Drosophila genome using TopHat on the Galaxy platform [[63]18, [64]19]. The differential expression of a total 15863 annotated nuclear-encoded genes was analyzed using Cufflinks and Cuffdiff on the Galaxy platform [[65]18, [66]19]. GO term and KEGG pathway enrichment analyses were performed with WebGestalt ([67]http://www.webgestalt.org) using 388 significantly differentially expressed genes. 2.4. qRT-PCR qRT-PCR was performed as previously described [[68]12, [69]13]. RNAs were treated with Turbo DNase (Thermo Fisher Scientific) and were reverse-transcribed into cDNA using an oligo-dT primer and AMV Reverse Transcriptase (NEB). qPCR was performed using SsoAdvanced Universal SYBR Green Supermix on CFX96 (Biorad) and results were analyzed using CFX Manager (Biorad). The primers used were as follows. gapdh, 5′-TGATGAAATTAAGGCCAAGGTTCAGGA-3′ and 5′- TCGTTGTCGTACCAAGAGATCAGCTTC-3′. rp49, 5′-CTGCCCACCGGATTCAAG-3′ and 5′-CGATCTCGCCGCAGTAAAC-3′. actin5c, 5′-AAGTTGCTGCTCTGGTTGTCG-3′ and 5′-GCCACACGCAGCTCATTGTAG-3′. CG11598, 5′-CCAAACACCTAAACACAAGCGT-3′ and 5′-GAATGGTCTGAACTTGCCGC-3′. CG43061, 5′-TGGGGAAAGTTGTAAGTCGCA-3′ and 5′-TAGTGGCCCAAGACAAACGG-3′. MtnD, 5′-GGAACAAACTGCCAGTGCTC-3′ and 5′-CCTTGGGTCCGTTCTAGCAG-3′. CG13905, 5′-CAGGGATTCCAAGCGGAGTT-3′ and 5′-GCCAGAGTTCCGACAGTTGA-3′. CG13215, 5′-TGCAACCAATCGCAGATCCT-3′ and 5′-ATTTGCCAAAGCTGTCGGTG-3′. 2.5. Deposited sequenced libraries The SRA accession number for the small RNA-seq and mRNA-seq libraries reported in this manuscript is SRP137802. 2.6. Statistical test Unpaired two-tailed Student's t-test was performed to determine statistical significance of changes in small RNA and mRNA abundance. The Bonferroni method and the Benjamini-Hochberg method were used to adjust P-values for multiple comparisons in small RNA-seq and mRNA-seq analyses, respectively. 3. Results 3.1. Small RNA-seq We performed high-throughput sequencing of small RNAs prepared from blanks homozygous (blanks^KG00804) and heterozygous (blanks^KG00804/TM3) mutant testes using three biological replicates each. Two miRNAs, miR-14-3p and miR-317-5p, were significantly downregulated (became 0.07-fold and 0.23-fold respectively) in blanks homozygous mutants compared with the heterozygous controls ([70]Fig. 2 and Supplementary Table S1). In addition, esi-2.1, the most abundant endo-siRNA in flies, was significantly downregulated (became 0.30-fold) in blanks homozygous mutants. The levels of other miRNAs and endo-siRNAs esi-1.1 and esi-1.2 were unchanged, consistent with a previously published study that showed no change in bantam and miR-7 levels [71][7]. No noticeable changes in piRNA levels, small RNAs expressed exclusively in gonadal tissues [72][20], were observed ([73]Fig. 3). Fig. 3. [74]Fig. 3 [75]Open in a new tab No change in piRNA abundance in blanks mutant testes. (A) Scatter plot of normalized number of reads of transposon-sense mapping piRNAs grouped based on transposon families. (B) Scatter plot of normalized number of reads of transposon-antisense mapping piRNAs grouped based on transposon families. 3.2. mRNA-seq We performed high-throughput sequencing of poly-A+ mRNAs prepared from blanks homozygous (blanks^KG00804) and heterozygous (blanks^KG00804/TM3) mutant testes (Supplementary Table S2). We identified 388 genes whose mRNA levels were significantly changed in blanks homozygous mutants compared with heterozygous mutant controls; 362 of them were downregulated and 26 upregulated ([76]Fig. 4A). The 24 genes with the highest fold changes in abundance are listed in [77]Fig. 4B (downregulated) and 4C (upregulated). The top 4 downregulated genes were Hsp70Ba (Heat-shock-protein-70Ba. Heat shock protein 70 chaperone), CG11598 (contains a partial AB-hydrolase lipase domain), CG42540 (band 7/mec-2 family protein), and CG43061 (contains a Conotoxin I-superfamily domain). The top 4 upregulated genes were MtnD (Metallothionein D. Metal ion binding protein), CG13905 (Insect allergen-related protein), CG13215, and CG9466 (Lysosomal α-mannosidase V). Fig. 4. [78]Fig. 4 [79]Open in a new tab mRNA-seq of blanks mutant testes. (A) Volcano plot showing the log[2](fold-change) and –log[10](p-value) of each mRNA level in blanks homozygous (blanks^KG00804) mutants relative to heterozygous controls (blanks^KG00804/TM3), determined by high-throughput sequencing of mRNAs. mRNA data with p-value <0.001 (= q-value <0.05) are depicted in red. (B) Top 24 downregulated genes. Genes whose differential expression was validated by qRT-PCR in (D) were indicated by *. (C) Top 24 upregulated genes. Genes whose differential expression was validated by qRT-PCR in (D) were indicated by *. (D) Relative abundance of mRNAs normalized by the gapdh mRNA level determined by qRT-PCR. Mean ± SD for three biological replicates. P-value <0.05 is indicated by * (Student's t-test). To validate differentially expressed mRNAs identified by mRNA-seq, we measured levels of some representative mRNAs using qRT-PCR. We chose 2 downregulated mRNAs and 3 upregulated mRNAs out of each top 4. Consistent with our mRNA-seq results, our qRT-PCR analysis showed that CG11598 and CG43061 mRNAs were significantly downregulated and MtnD, CG13905, and CG13215 mRNAs were significantly upregulated in the blanks homozygous mutants compared with the heterozygous controls ([80]Fig. 4D). As a negative control, we measured rp49 mRNA level and it did not exhibit significant change between the blanks homozygous mutant and control ([81]Fig. 4D), consistent with our mRNA-seq analysis (Supplementary Table S2). Thus, our qRT-PCR analysis validated the mRNA-seq results. 3.3. GO term/KEGG pathway enrichment analysis Using the 388 differentially expressed genes in the blanks homozygous mutat testes, we performed GO term/KEGG pathway enrichment analyses to determine pathways altered in blanks homozygous mutant testes ([82]Table 1). GO term biological processes of catabolic and metabolic processes, anion and cation transport, and mating and reproductive behavior were enriched. GO term cellular components of plasma membrane and lipid particle were enriched. GO term molecular functions of transmembrane transporter activity, oxidoreductase activity, carbohydrate binding, lipase activity, iron ion binding, serine hydrolase activity, and carboxylic ester hydrolase activity, were enriched. KEGG pathways of starch and sucrose metabolism, metabolic pathways, and pentose and glucuronate interconversions were enriched. Table 1. GO term enrichment analysis of dysregulated mRNAs. p-value GOterm Biological Process Multi-multicellular organism process 0.0000004 Single-organism catabolic process 0.0000376 Anion transport 0.0000549 Lipid metabolic process 0.0000743 Mating 0.0004445 Reproductive behavior 0.0005331 Organic acid metabolic process 0.0005615 Cation transport 0.0010928 GOterm Cellular Component Plasma membrane region 0.0002045 Lipid particle 0.0003110 GOterm Molecular Function Active transmembrane transporter activity 0.0000091 Anion transmembrane transporter activity 0.0000211 Oxidoreductase activity, acting on CH-OH group of donors 0.0000588 Carbohydrate binding 0.0000763 Lipase activity 0.0001259 Cation transmembrane transporter activity 0.0004441 Cofactor binding 0.0004733 Iron ion binding 0.0013241 Serine hydrolase activity 0.0016117 Carboxylic ester hydrolase activity 0.0016939 Oxidoreductase activity, acting on single donors with incorporation of molecular oxygen 0.0018319 Oxidoreductase activity, acting on the aldehyde or oxo group of donors 0.0018756 Oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen 0.0032525 KEGG pathway Starch and sucrose metabolism 0.0001429 Metabolic pathways 0.0002715 Pentose and glucuronate interconversions 0.0005717 [83]Open in a new tab Significantly enriched GO terms and KEGG pathways associated with differentially expressed genes identified in [84]Fig. 4A. 4. Discussion Our high-throughput sequencing studies identified specific small RNAs and mRNAs which are dysregulated in blanks mutant testes, thereby demonstrating that Blanks plays important roles in regulating small RNA and mRNA profiles in testes. Our small RNA-seq study identified two miRNAs (miR-14-3p and miR-317-5p) and one endo-siRNA (esi-2.1) downregulated in blanks mutant testes ([85]Fig. 2), suggesting that Blanks is directly or indirectly required for efficient production and/or stability of these small RNAs. Previous studies showed that a Blanks protein complex can bind siRNAs in vitro [86][6]. Therefore, Blanks might bind and stabilize miR-14-3p, miR-317-5p, and esi-2.1. None of the predicted target mRNAs of miR-14-3p, miR-317-5p, and esi-2.1 were upregulated in blanks mutant (Supplemental Table S2). As miRNAs can cause both translational suppression and mRNA destabilization, which mode is predominant or occurs upstream depends on biological situation [[87]21, [88]22, [89]23, [90]24]. It will be interesting to examine whether the protein abundance of predicted targets of miR-14-3p, miR-317-5p, and esi-2.1 are upregulated in blanks mutant testes. Our mRNA-seq study identified 26 genes whose mRNAs are upregulated and 362 genes whose mRNAs are downregulated in blanks mutant testes ([91]Fig. 3). We validated the changes in mRNAs using qRT-PCR. Blanks may play roles in regulating these mRNA levels transcriptionally and/or post-transcriptionally. Blanks interacts with heterochromatin factors and insulators [[92]6, [93]25, [94]26] and thereby may regulate transcription of these mRNAs. Alternatively, at the post-transcriptional level, Blanks might bind these mRNAs and change their levels. These dysregulated genes are significantly associated with GO terms of catabolic and metabolic processes, anion and cation transport, and mating and reproductive behavior, suggesting that these biological processes are dysregulated in blanks mutant testes. These pathways may underlie previously reported defects in post-meiotic spermiogenesis and male-sterility of blanks mutant flies [[95]6, [96]7]. Future studies, including identification of RNAs bound by Blanks, are required to understand the molecular mechanisms by which Blanks regulate testicular small RNA and mRNA profiles. Declarations Author contribution statement Susan E. Liao: Conceived and designed the experiments; Performed the experiments; Analyzed and interpreted the data; Wrote the paper. Yiwei Ai: Performed the experiments. Ryuya Fukunaga: Conceived and designed the experiments; Performed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper. Funding statement This work was supported by the grants from American Heart Association [15SDG23220028] and the National Institutes of Health [R01GM116841] to RF. Competing interest statement The authors declare no conflicts of interest. Additional information No additional information is available for this paper. Acknowledgements