http://orcid.org/0000-0003-3054-0123
Figures
Abstract
The European beaver (Castor fiber L.) is an important free-living rodent that inhabits Eurasian temperate forests. Beavers are often referred to as ecosystem engineers because they create or change existing habitats, enhance biodiversity and prepare the environment for diverse plant and animal species. Beavers are protected in most European Union countries, but their genomic background remains unknown. In this study, gene expression patterns in beaver testes and the variations in genetic expression in breeding and non-breeding seasons were determined by high-throughput transcriptome sequencing. Paired-end sequencing in the Illumina HiSeq 2000 sequencer produced a total of 373.06 million of high-quality reads. De novo assembly of contigs yielded 130,741 unigenes with an average length of 1,369.3 nt, N50 value of 1,734, and average GC content of 46.51%. A comprehensive analysis of the testicular transcriptome revealed more than 26,000 highly expressed unigenes which exhibited the highest homology with Rattus norvegicus and Ictidomys tridecemlineatus genomes. More than 8,000 highly expressed genes were found to be involved in fundamental biological processes, cellular components or molecular pathways. The study also revealed 42 genes whose regulation differed between breeding and non-breeding seasons. During the non-breeding period, the expression of 37 genes was up-regulated, and the expression of 5 genes was down-regulated relative to the breeding season. The identified genes encode molecules which are involved in signaling transduction, DNA repair, stress responses, inflammatory processes, metabolism and steroidogenesis. Our results pave the way for further research into season-dependent variations in beaver testes.
Citation: Bogacka I, Paukszto Ł, Jastrzębski JP, Czerwińska J, Chojnowska K, Kamińska B, et al. (2017) Seasonal differences in the testicular transcriptome profile of free-living European beavers (Castor fiber L.) determined by the RNA-Seq method. PLoS ONE 12(7): e0180323. https://doi.org/10.1371/journal.pone.0180323
Editor: Maoteng Li, Huazhong University of Science and Technology, CHINA
Received: February 10, 2017; Accepted: June 14, 2017; Published: July 5, 2017
Copyright: © 2017 Bogacka et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: All relevant data are within the paper and its Supporting Information files. The raw sequencing data are available in the National Center for Biotechnology Information Short Read Archive (http://www.ncbi.nlm.nih.gov/Traces/sra/sra.cgi) under accession PRJNA355849.
Funding: This study was supported by the Polish National Science Centre (project No. 2012/07/B/NZ9/01335).
Competing interests: The authors have declared that no competing interests exist.
Introduction
The European beaver (Castor fiber L.) is the largest rodent in Eurasia and the world’s second largest rodent after the capybara. In the past, the species was widely distributed across Eurasian forests, between western Europe and eastern Siberia, but it was nearly driven to extinction after centuries of hunting for its fur and castoreum, the exudate from the beaver scent gland. By the end of the 19th century, hunting pressure combined with the decline in beaver habitats reduced the Eurasian beaver population to 1,200 individuals in several remote regions [1]. Fortunately, the introduction of a hunting ban and other conservation measures enabled the beaver population to survive and spread to other regions through natural and artificial recolonization.
The beaver population in Europe has increased significantly in the last two decades, and the conservation status of the species on the Red List of Threatened Species of the International Union for Conservation of Nature (IUCN) has been changed from endangered to least concern. The beaver is often referred to as an ecosystem engineer in Eurasian temperate forests [2, 3, 4] due to its ability to create or modify habitats. Beavers enhance biodiversity and prepare the ecosystem for the emergence of other species of plants, water-dwelling invertebrates, fish, amphibians, reptiles, birds and mammals [3–6]. Beaver dams retain sediments and organic material [3], they influence the nutrient cycle and speed up the decomposition of organic matter. The dams also induce changes in the chemical and physical parameters of water [3, 7].
Beavers are monogamous animals, and both males and females co-parent their young. Beavers usually live in small family groups of 2–14 individuals which consist of an adult pair, offspring younger than one year, yearlings and, sporadically, subadults [8]. They have a seasonal pattern of reproduction. Similarly to hamsters and horses, beavers are long-day breeders, and their cyclicity is dependent on photoperiod. The reproductive activity of beavers peaks in winter, and the animals mate under the ice between January and March. In summer, beavers take care of their offspring and store food reserves. Beavers do not hibernate, and they accumulate subcutaneous fat depots in autumn [1, 9]. Both males and females exhibit reproductive seasonality, but males appear to have a longer breeding season because the completion of spermatogenesis usually requires more time than oogenesis [10] and earlier recrudescence of male breeding activity facilitates sexual competition for mates [11]. Testis size, the ability to produce gametes, hormonal activity and gene expression change throughout the year in seasonally active males, especially in wild animals [12–15].
In view of the above findings, we hypothesized that changes in gene expression between breeding and non-breeding seasons can also be detected in beaver testes. In the presented study, beaver testes harvested in two periods of reproductive activity were subjected to comprehensive transcriptomic analyses with the use of the RNA-Seq approach. A differential expression analysis of testes harvested in breeding (April) and non-breeding (November) seasons was carried out to determine the impact of reproductive activity on transcriptome profiles. Differentially expressed genes and their products were linked with biological functions.
Materials and methods
Animals
The study was performed on European beavers captured in various habitats in the Region of Warmia and Mazury in north-eastern Poland, as described in our previous studies [16, 17]. The animals were captured upon the prior consent of the Regional Directorate for Environmental Protection in Olsztyn, Poland (decision No. RDOS-28-OOP-6631-0007-638/09/10/pj). All procedures were conducted in accordance with the ethical standards of the Institutional Animal Ethics Committees (local approvals: SGGW/11/2010 and UWM/87/2012/DTN). Beavers were captured by the same group of qualified staff over a period of two years (2011–2012), in the same period of the year, during two different stages of reproductive activity: in April–‘breeding period’ (n = 4) and November–‘non-breeding’, sexual silence (n = 4). Every experimental group consisted of different individuals. The animals were placed in cages and transported to the laboratory of the Research Station for Ecological Agriculture and Preservation of Native Breeds of the Polish Academy of Sciences in Popielno. Each beaver was weighed in the cage, and the weight of an empty cage was subtracted. Beavers were anesthetized by injection of xylazine (3 mg/kg BW; Sedazin®, Biovet Puławy, Poland) and ketamine (15 mg/kg BW; Bioketan, Vetoquinol Biowet, Poland) and sacrificed by exsanguination. The testes were sampled, frozen immediately in liquid nitrogen and stored at -80°C until total RNA isolation.
RNA isolation and transcriptome sequencing
Total RNA for transcriptome sequencing was isolated from four testicular compartments (n = 2 for each experimental group) using the Qiagen RNeasy Kit and the Qiagen RNase-Free DNase Set for cleanup (Qiagen, USA) according to the manufacturer’s recommendations. The concentration and purity of total RNA isolated from testicular tissues were analyzed with the Tecan Infinite M200 plate reader (Tecan Group Ltd, Switzerland), and RIN values were evaluated in the Agilent Bioanalyzer 2100 (Agilent Technology, USA). Total RNA was isolated using an OpenExome kit (Poland) to prepare sequencing libraries. The cDNA library was developed based on total RNA with the TruSeq Stranded mRNA LT Sample Prep Kit v3, using the appropriate index adapters for each sample. The concentration of mRNA ranged from 305 to 348 ng/μl. Libraries were pooled and sequenced to produce 2 x 100 bp paired-end (PE) reads in the Illumina HiSeq 2000 high-throughput sequencer. All samples were sequenced in a single lane with human reference mRNA as positive control. Libraries were quantified with the use of the KAPA library quantification kit.
De novo transcriptome assembly
The quality of raw paired-end reads was evaluated with the use of the FastQC software v.0.10.0 (www.bioinformatics.babraham.ac.uk). Adaptor sequences were removed from raw sequencing reads with the AlienTrimmer tool in default mode [18]. Low-quality reads (PHRED <20) were removed from the dataset. Trimmed sequences from 4 testis samples were de novo assembled in the Trinity software v. r20140717 [19] with a 20-core processor and a 90 GB RAM server of the Regional IT Center (Olsztyn, Poland). The reconstructed contigs were used to develop a reference transcriptome of Castor fiber with a minimum transcript length of 500 nt. De novo assembly was performed with the Trinitystats.pl script, and the assembly statistics included average contig length, GC content, total assembled bases and N50 parameters. Duplicate transcripts were removed from the reference transcriptome in the CD-HIT-EST v.4.6 software (http://weizhong-lab.ucsd.edu/cd-hit/) to minimize data redundancy (by clustering sequences at 95% identity) [20]. The Benchmarking Universal Single-Copy Orthologs (BUSCO) v.1.1 tool [21] was used to check the completion of transcriptome assembly based on the percentage of transcripts identified as putative core eukaryotic genes (CEGs).
Comprehensive transcriptome analysis
The assembled transcripts with FPKM > 2 in at least two probes were used in a comprehensive analysis because this approach proved to be highly effective in other studies [22, 23]. Selected contigs were blasted (blastx, cut-off E value of 1e-5) against five ENSEMBL databases of protein sequences from Mus musculus, Rattus norvegicus, Homo sapiens, Ictidomys tridecemlineatus and Dipodomys ordii. The highly expressed transcripts (FPKM >2) were aligned with the NCBI database of rodent proteins with the use of the CloudBlast (blastx-fast) tool implemented in BLAST2GO to obtain genes specific for beaver testes [24]. Transcripts were annotated to three main gene ontology categories: biological process (BP), cellular component (CC) and molecular function (MF). The search for conserved protein motifs and the prediction of protein families were performed in the InterProScan v.5 program [25] for each highly expressed transcript.
The obtained beaver contigs were compared with the rat genome to determine the abundance and distribution of transcripts. Contigs with FPKM >2 and longer than 500 bp were used as search queries in blastx against the NCBI nr protein database to identify testicular unigenes. Unigenes were annotated with the BLAST tool against the nr database with cut-off E-value of 10−5. The obtained unigenes were mapped onto the Rattus norvegicus genome, a model organism that is most closely related to the beaver, with the use of the GMAP software [26]. Each unigene was compared with the rat genome. Based on this alignment, the reciprocal best-hits method with a 100 kb window was used to determine the distribution of beaver transcripts across the rat genome. The RBH graph was generated in the Circos program [27].
Analysis of gene expression
Trimmed paired-end reads of each sample were realigned to a non-redundant beaver transcriptome in the Bowtie software [28]. Mapped contigs were filtered to remove transfrags with low expression levels. Non-normalized read counts and TMM normalized FPKM were calculated in the eXpress software [29] to estimate transcript expression levels.
The expression of contigs with absolute log fold-change (logFC) of 2 or higher and the corrected p-value (FDR) of less than 0.05 was significantly different. Two statistical methods were used to confirm differentially expressed transcripts in the testes (autumn vs. spring): DESeq [30] and edgeR [31]. Differentially expressed contigs were aligned to the NCBI nr protein database with blastx (cut-off E-value of 10−5) to identify the coding DNA sequence (CDS). Differentially expressed genes (DEGs) were visualized in a heatmap plot with gplots in R package (http://www.r-project.org/). All DEGs were compared against the Homo sapiens database in the KEGG Orthology Based Annotation System (KOBAS) for metabolic pathway and Gene Ontology enrichment analyses with p-value < 0.05 and q-value < 0.9.
RNA isolation, cDNA synthesis and qRT-PCR for validation of transcriptome sequencing data
Total RNA was isolated from beaver testes collected during breeding (n = 4) and non-breeding seasons (n = 4) using the A&A mini column kit (A&A Biotechnology, Poland) with a DNase treatment step. The concentration and the amount of the isolated RNA were determined spectrophotometrically (Infinite M200 PRO, Tecan, Switzerland), and RNA integrity was verified on 1.5% agarose gel. One μg of RNA was reverse transcribed into cDNA with a total volume of 20 μL with the use of the QuantiTect® Reverse Transcription Kit (Qiagen, USA). The reaction included two steps: preliminary genomic DNA elimination (2 min at 42°C) and main reverse transcription (30 min at 42°C). The reaction was terminated by incubating the samples for 3 min at 95°C in a thermal cycler (Eppendorf, Germany).
To validate transcriptome sequencing results, the expression of differentially regulated genes, such as PRRKG1, DSG2, SMARCA2, BMX and AGT, was determined by qRT-PCR with the use of the 7300 PCR System and the Power SYBR Green Master Mix (Applied Biosystems, USA). The constitutively expressed ACTB and GAPDH genes were used as reference (housekeeping) genes according to a previously described procedure [32, 33]. Specific primer sequences (S1 Table) for amplifying target and reference genes were designed in the Primer Express Software 3 (Applied Biosystems). PCR reaction mixtures with a final volume of 25 μL consisted of cDNA (20 ng for target genes and 5 ng for reference genes), 400 nM of the primers, 12.5 μL of the Power SYBR Green PCR Master Mix (Applied Biosystems, USA), and RNase-free water. The qRT-PCR reaction conditions were as follows: enzyme activation and initial denaturation at 95°C for 10 min, followed by 40 cycles of denaturation at 95°C for 15 s, annealing for 1 min at 60°C for all tested genes. In no template controls (NTC), cDNA was replaced with RNase-free water. The efficiency of qRT-PCR for each tested cDNA was determined at 100%. All samples were run in duplicates. The specificity of amplification was tested at the end of the reaction by analyzing the melting-curve. The relative expression of the tested genes was calculated with the use of the comparative cycle threshold method (ΔΔCt) described by Schmittgen & Livak [34], where ΔΔCt is obtained by subtracting ΔCt of the geometric mean, calculated from the Ct value of both reference genes (GAPDH and ACTB), from the corresponding ΔCt of each experimental sample of the target gene. The group with the lowest expression was selected as the calibrator. The results of qRT-PCR were processed statistically in the Statistica software (Statoft Inc. Tulsa, USA) with the use of the Student’s t-test and were expressed as means ± SEM. The results were regarded as statistically significant at p < 0.05.
Results
De novo transcriptome assembly
A total of 380.22 (2 x 190.11) million of paired-end reads were sequenced in the Illumina HiSeq 2000 sequencer for the de novo assembly of the transcriptome of Castor fiber testicular tissue (Table 1).
After processing of raw reads (adaptor sequences were cut and low-quality reads with Phred value < 20 were removed), 373.06 million (98.12%) clean reads were obtained and assembled de novo using the Trinity software. Assembly results are summarized in Table 2. The assembly was performed on 260,311 contigs. To minimize redundancy, contigs with sequence similarity higher than 95% were clustered in 231,045 non-redundant contigs grouped into 130,741 unigenes. Each unigene was described by a group of related transcripts included in the same de Bruijn graph, and each cluster was regarded as a single gene. Non-redundant contigs ranged from 500 nt to 20,720 nt, with N50 value (the length of the contig describing half of the sum of the lengths of all contigs) of 1,734 nt. The longest transcript was identified as fibrous sheath interacting protein 2 with 99% coverage and 69% sequence identity with Marmota marmota (using blastx, NR Rodentia). More than 100 k contigs were longer than 1000 nt. A total of 39,657 unigenes with ORFs longer than 300 nt were detected.
Comprehensive analysis of beaver testicular transcriptome
In the set of all non-redundant contigs, 88% (2663/3023) complete and 4.4% (136/3023) fragmented core eukaryotic (vertebrate) genes were found. More than 92% of the predicted housekeeping genes confirmed the completeness of transcriptome assembly.
The whole set of clean reads was aligned back (with MAPQ >30) to the non-redundant Castor fiber transcriptome, and 88.7% (Cf_a1) to 91.5% (Cf_s2) of reads were successfully mapped (Table 2). These data were used for transcriptome profiling (comprehensive analysis) and differential expression analyses. Gene counts and TMM normalized FPKM values were calculated with the eXpress tool, and transcript expression levels were estimated. To remove weakly expressed transcripts from assembly datasets, only contigs with FPKM > 2 in at least 50% of the samples were used in comprehensive analysis. In the group of 26,228 unigenes that fulfilled the above requirement, the longest isoforms were extracted and blasted (cut-off E-value of 10e-5) with five protein Ensembl databases. The blastx search revealed 14,734, 14,696, 14,043, 13,072 and 14,065 highly expressed homologs for Rattus norvegicus, Ictidomys tridecemlineatus, Mus musculus, Dipodomys ordii and Homo sapiens, respectively. In this group, 258 corresponding homologs were identified exclusively in rat, 33 in squirrel, 16 in mouse, and 1 in kangaroo rat. The remaining 172 beaver transcripts were not found in Rodentia, and they had only human homologs (Fig 1).
Based on sequence homology, Castor fiber transcripts (FPKM >2) were matched to 5 species (4 Rodentia species and Homo sapiens). The numbers in the overlapping areas denote BLASTX matches with two or more species. The numbers in the non-overlapping areas denote uniquely blasted homologs.
Highly expressed transcripts were compared with the rat genome (the best described Rodentia genome with the closest phylogenetic relationship with Castor fiber). A total of 15,774 unigenes were annotated using the NCBI nr database. The results of contig mapping to the rat genome are visualized in Fig 2. Approximately 60.14% of highly expressed transcripts longer than 500 bp with FPKM>2 were mapped to orthologs localized on the rat genome. A comprehensive analysis revealed that transcriptome assembly had a broad representation and could be used as the reference transcriptome for Castor fiber testicular tissue.
The outer track shows the density of mapped unigenes of C. fiber to the R. norvegicus genome, in both forward (orange histogram) and reverse (green histogram) strands. The middle track represents exome density throughout R. norvegicus chromosomes in both forward (blue histogram) and reverse (red histogram) strands. The dark blue circles in the inner track depict the percentage of all unigenes aligned to each chromosome (the size of circle corresponds to the percentage of aligned unigenes), whereas the light blue circles depict the density of identified unigenes throughout chromosomes (the size of circle corresponds to the amount of aligned unigenes per chromosome length). The chromosomes length is expressed in million base pairs (Mb).
Functional classification
The study revealed 8,188 highly expressed genes mapped to GO terms, where 3918 were assigned to the Biological Process sub-ontology (Fig 3). The annotated pathways included metabolic processes represented by 212 unigenes (5.41%), transcription, DNA-templated–represented by 200 unigenes (5.10%), regulation of transcription, DNA-templated–represented by 160 unigenes (4.08%), and positive and negative regulation of transcription from RNA polymerase–represented by 190 and 149 unigenes, respectively (4.85% and 3.80%). The most interesting top-hit biological process was spermatogenesis, represented by 53 corresponding unigenes (S2 Table). Other biological processes related to the male gonadal tissue included 23 (spermatoid development), 17 (steroid hormone mediated signaling pathway), 15 (sperm motility), 13 (male gonad development), 5 (Sertoli cell development), 3 (Sertoli cell proliferation), 3 (sperm axoneme assembly) and 2 (spermatid differentiation) highly expressed unigenes (S2 Table).
Annotations (received by Blast2GO tool from ‘no.’ NCBI database limited to [rodentia, taxa:9989] using blastx-fast algorithm) are divided into three GO categories: Biological Process (the upper part, in blue), Cellular Component (the middle part, in red), Molecular Function (the lower part, in green). Values on the vertical axis indicate the number of sequences assigned to GO terms. Top 25 GO terms are presented on the vertical axis.
In the cellular component category, the highest number of transcripts matched the term integral component of membrane (793), cytoplasm (713), nucleus (641) and extracellular exosome (579). The molecular function category consisted of 1885 terms which were assigned to 3819 unigenes. In this category, the majority of unigenes (457) were associated with ATP binding (Fig 3).
A total of 423 unigenes were classified into 116 KEGG pathways (S3 Table). The top three KEGG pathways involved purine metabolism (312, 73.75%), thiamine metabolism (245, 57.91%) and aminobenzoate degradation (73, 17.26%) (S1 Fig). A metabolic pathway analysis revealed 12 highly expressed transcripts that were involved in steroid hormone biosynthesis, 3 transcripts that were involved in steroid biosynthesis and 2 transcripts that were involved in steroid degradation (S4 Table).
The majority of functional annotation results (3,699 out of 8,188) matched Ictidomys tridecemlineatus. InterPro domain annotations showed functional information for 13,069 unigenes. Most of the unigenes annotated in the InterProScan were found in protein domain and functional databases: 10,578 PFAM, 3249 PANTHER and 5279 SMART. The identified proteins (encoded by unigenes) also encompassed 115 oxidoreductases, 403 transferases, 326 hydrolases, 30 lyases, 35 isomerases and 54 ligases.
Seasonal differences in beaver testicular transcriptome
A comparison of the testes harvested during breeding (April) and non-breeding (November) seasons revealed 152 transcripts (including 42 protein coding genes) that were significantly differentially expressed (log2FC >1; FDR 0.05) (Fig 4). The remaining non-annotated differentially expressed contigs were putative noncoding RNAs or novel genes. Of these, 37 genes were up-regulated and 5 genes were down-regulated in beaver testes from the non-breeding season relative to the samples obtained in the breeding season (Table 3).
MA plot (A) and volcano plot (B) describe log2 fold change between spring and autumn according to DESeq. Each point represents one contig. The red points represent the contigs of p adjusted < 0.05 and log2 fold change > 1. (C)–represents the Venn diagram of the quantity of DEGs calculated by both methods: edgeR and DESeq (log2 fold change > 1, p adjusted < 0.05). Values in the "up" line denote the number of up-regulated genes, values in the "down" line denote the number of down-regulated genes, and values in the "all" line denote the total number of down- and up-regulated genes. (D)–represents the matrix of Pearson correlation coefficients calculated for 152 contigs designated as significantly expressed using the DESeq statistical method. Blue color represents high correlation and red–low samples correlation. Samples are clustered by both rows and columns. Cyan and red bars indicate the samples grouped by seasons. E–represents the heat map of differences in the expression of 42 genes between spring and autumn. The rows represent selected genes; in each column, individual samples are grouped into two seasons (upper dendrogram): spring (red bar) and autumn (cyan bar). The colors and intensity of fields in the heat map indicate transformed fold change (log2 fold change values from -3,767 to 3,767): red—underexpression, green—overexpression. The gene order was generated automatically in the dendrogram (left side of the graph).
The best-saturated mammalian functional annotations (human reference) were used to associate DEGs with 2041 GO terms and 28 DEGs with 41 KEGG (http://www.genome.jp/kegg-bin/) and 208 Reactome (http://www.reactome.org/cgi-bin) metabolic pathway processes (S5 Table and S6 Table).
Comparison of RNA-Seq and qRT-PCR values
An analysis of mRNA abundance in genes evaluated by qRT-PCR revealed a correlation with RNA-Seq results (Fig 5). The mRNA levels PRKG1, DSG2, SMARCA2 and BMX were higher and AGT expression was lower in November than in April.
The expression of PRKG1, DSG2, SMARCA2 and BMX genes was up-regulated, whereas AGT gene expression was down-regulated in beaver testes during the non-breeding season relative to the breeding season. Different letters indicate statistically significant differences (p < 0.05) between the tested groups of beavers.
Discussion
To the best of our knowledge, this is the first study providing a comprehensive analysis of transcriptome profiles of European beaver testes with the use of the RNA-Seq technique. It is also the first study to describe variations in transcriptome profiles between breeding and non-breading seasons. The RNA-Seq technique is a powerful tool for generating large-scale transcriptome data that were useful for presenting differentially expressed genes in various species, tissues or cell types. According to some studies, Illumina sequencing data are replicable with relatively little technical variation [35, 36].
In the current study, we obtained a total of 373.06 million high-quality reads (98.12% Q2 bases). De novo assembly of contigs yielded 130,741 unigenes which can be useful in future studies of functional genomics in the European beaver. The average length of unigenes with N50 value of 1,734 was 1,369.3, and it was higher than that reported by other researchers [37–40]. The above indicates that our transcriptome sequencing results were successfully assembled and that they are of high quality.
A comprehensive analysis of the testicular transcriptome revealed more than 26,000 highly expressed unigenes which exhibited the highest homology with Rattus norvegicus and Ictidomys tridecemlineatus genomes and lower homology with other rodent species such as Dipodomys ordii or Mus musculus. More than 8,000 highly expressed beaver genes were found to be involved in fundamental biological processes, cellular components or molecular pathways. The largest group of genes encoding biological processes comprises genes responsible for the regulation of metabolic processes, transcription and regulation of transcription. A set of genes encoding male reproductive functions was also identified. This group comprises genes controlling the development and proliferation of Sertoli cells, spermatoid differentiation, sperm development, sperm motility and the steroid hormone mediated signaling pathway. Some testis-specific markers involved in beaver reproductive functions could be omitted if they did not meet stringent search criteria, i.e. if they did not show at least a two-fold change in the expression level.
The present study revealed 42 genes that were differentially regulated in the testes during the breeding season than in the non-breeding season. Surprisingly, during the non-breeding season (November), a much larger group of genes (37) was significantly more highly expressed, and a smaller group of genes (5) was less expressed than during the breeding season (April). The up-regulation of a larger number of genes in mouse testes after androgen withdrawal was determined in a microarray analysis by Petrusz et al. [41] and Sadate-Ngatchou et al. [42]. In the present study, the observed variations in the transcript profiles of beaver testes could be linked with plasma androgens concentrations. However, in our previous study of beavers, plasma testosterone levels did not differ significantly between the analyzed reproductive seasons [16].
In the group of differentially regulated high-impact genes, 18 genes encoded signaling molecules, 7 genes encoded transcription factors and DNA repair molecules, 9 genes encoded cell surface proteins, 10 genes encoded stress response or inflammatory process, 12 genes encoded ion channels and 4 genes encoded extracellular matrix components. Only the most interesting genes are briefly discussed below.
Several up-regulated transcripts from the non-breading season could be involved in the regulation of spermatogenesis in beaver testes, including VLDLR which belongs to the low-density lipoprotein receptor (LDLR) family. It has been reported that VLDLR overexpression in germ cells disrupted spermatogenesis in transgenic mice [43], which suggests that the VLDLR transgene inhibits the development of sperm cells. Similar conclusions could be drawn from our results which showed markedly higher transcript levels in beaver testes during the non-breeding season when spermatogenesis is limited.
SMAD proteins are intracellular mediators of the transforming growth factor-β (TGFβ). These transcription factors are implicated in bone morphogenetic protein (BMP) and activin signaling, they are highly involved in the regulation of male fertility, and they influence germ and somatic cells during fetal and postnatal life [44, 45]. In 15-day-old mice, the SMAD5 transcript was identified by in situ hybridization in spermatogonia, whereas a weak signal was detected in spermatocytes. In adult testes, the signal was strongest in spermatogonia and spermatocytes and weaker in round spermatids and in Sertoli cells [46]. Interestingly, our findings revealed up-regulation of the SMAD5 transcript in beaver testes during the non-breeding season, which could indicate that the BMPs/SMADs pathway could play a role in seasonally dependent reproductive activity of the beaver.
The present findings indicate seasonally dependent changes in the expression of genes encoding transporter proteins. The expression of SLC7A1 and KCNMB2 (BK channels) in beaver testes was significantly higher in November. The SLC7A1 gene mediates the transport of cationic amino acids across cell membrane. It has been reported that rat Sertoli cells express SLC7A1 which relies on the NA+-independent transport system to deliver arginine or other cationic amino acids to germ cells [47]. The KCNMB2 gene is a calcium-gated potassium channel that has been described in various mammalian endocrine cells, including human and hamster testicular Leydig cells [48]. It has been suggested that these channels contribute to Leydig cell steroidogenesis. Interestingly, the presence of iberiotoxin, a specific channel blocker, did not induce changes in testosterone production by hamster Leydig cells in vitro under basal conditions, but a significant increase in testosterone levels was reported when hCG was added to culture media [48]. According to Gong et al. [49], these ion channels could also play an important role in spermatogenesis control. Sperm cells could possess a Ca2+-activated K+ channel which has been implicated in sperm activation and gamete interaction [50]. The above findings indicate that changes in the membrane potential of germ cells could be an important element of the signaling mechanism.
The large heparan sulfate proteoglycan (HSPG2, perlecan) is an extracellular matrix component that is normally expressed in the basement membrane (BM) underlying epithelial and endothelial cells. It was detected in the seminiferous tubule BM and in the interstitium, where the protein was localized around blood vessels and Leydig cells [51]. HSPG2 could also play an important role in self-renewal and differentiation of spermatogonial stem cells in the testes [51, 52].
Receptor tyrosine kinase ErbB4 has been detected on circulating human monocytes and neuronal macrophages, which points to its involvement in inflammatory processes. ErbB4 was also highly expressed in mouse Sertoli cells and in germ and Leydig cells [53]. Interestingly, targeted inactivation of ErbB4 in Sertoli cells induced a reduction in testis size, decreased testicular androgen production and delayed spermatogenesis [53].
The down-regulated transcripts identified in our study encode proteins that could be involved in stress responses and inflammatory processes. IRAK1BP1 is an inhibitory component of TLR, IL-1 and TNFR-related pathways [54, 55]. Overexpression of IRKA1BP1 in macrophages increased the expression of anti-inflammatory IL-10, but decreased the synthesis of proinflammatory IL-6. IRKA1BP1 also contributes to LPS-induced tolerance by influencing NF-ĸB [55].
In this study, the identified down-regulated genes also included angiotensinogen gene (AGT) which is a part of the rennin-angiotensin system. In rat testes, only renin mRNA was detected, whereas both renin and angiotensinogen mRNA were found in mouse testes [56]. The AGT transcript was identified as one of the up-regulated genes (2.64-fold increase, microarray analysis) in the testes of transgenic mice carrying the androgen-binding protein (ABP) gene [41]. The importance of the angiotensin-converting enzyme (ACE, a component of the rennin-angiotensin system) for sperm functions has been demonstrated in male mice by targeting genes that lack somatic ACE but retain testicular ACE [57]. Interestingly, we found that the LNPEP (IRAP) gene, which is involved in the rennin-angiotensin system, was up-regulated in beaver testes outside the breeding season. The LNPEP gene encodes zinc-dependent aminopeptidase, cleaves vasopressin, oxytocin and bradykinin, and catalyzes the final step in the conversion of angiotensinogen to angiotensin IV. The results of our study indicate that angiotensin could be produced locally in beaver testes, and they provide important insights into the biology of the renin-angiotensin system in this species.
In summary, this is the first study to generate a representative testicular transcriptome of the European beaver with the use of the powerful RNA-Seq technique. It is also the first study to describe seasonally dependent variations in the transcriptome profiles of beaver testes. We identified a set of 42 genes encoding molecules that are involved in signal transduction, DNA repair, stress responses, inflammatory processes, metabolism and steroidogenesis. Our findings pave the way for further research into the processes that occur in beaver testes during various periods of reproductive activity.
Supporting information
S1 Fig. KEGG analysis for highly expressed unigenes.
https://doi.org/10.1371/journal.pone.0180323.s001
(TIF)
S1 Table. Characteristics of primers used for real time PCR to validate RNA-Seq transcriptome sequencing.
https://doi.org/10.1371/journal.pone.0180323.s002
(DOCX)
S2 Table. Highly expressed genes (FPKM >2) involved in the biological processes.
Gene Ontology terms related to spermatogenesis, steroid hormone mediated signaling pathway, sperm motility, male gonad development, Sertoli cell development, spermatid development, spermatid differentiation, sperm axoneme assembly and Sertoli cell proliferation using the Blastx-fast algorithm in the Blast2GO tool. The columns contain the following information: “SeqName”–the ID of contigs, “Description”–full protein names, “Length”–contig length values, “sim mean”–average similarity between the aligned sequences.
https://doi.org/10.1371/journal.pone.0180323.s003
(DOCX)
S3 Table. Highly expressed unigenes (FPKM >2) assigned to 116 KEGG pathways.
https://doi.org/10.1371/journal.pone.0180323.s004
(DOCX)
S4 Table. Highly expressed genes (FPKM >2) associated with the steroid hormone biosynthesis KEGG pathway (Blastx-fast algorithm in the Blast2GO tool).
The columns contain the following information: “SeqName”–the ID of contigs, “Description”–full protein names, “Length”–contig length values, “sim mean”–average similarity between the aligned sequences.
https://doi.org/10.1371/journal.pone.0180323.s005
(DOCX)
S5 Table. Functional annotation of differentially expressed genes based on Blast2GO results.
https://doi.org/10.1371/journal.pone.0180323.s006
(DOCX)
S6 Table. Functional annotation of differentially expressed genes based on GO, KEGG and Reactome human databases.
The columns contain the following information: “Term”–names of terms in each database, “Database”–database name, “ID”–the id of terms in each database, “Input number”–number of assigned unigenes, “Background number”–number of all genes assigned to a given term in each database, “p-Value”,”Corrected p-Value”–statistical significance of enrichment, “Input”–short names of assigned DEGs.
https://doi.org/10.1371/journal.pone.0180323.s007
(XLS)
Acknowledgments
The authors would like to thank Jan Goździewski of the Polish Hunting Association in Suwałki for capturing and delivering the animals.
Author Contributions
- Conceptualization: IB TK BK.
- Data curation: JPJ ŁP.
- Formal analysis: JPJ ŁP JC KC.
- Funding acquisition: TK IB BK.
- Investigation: JC KC AK NS.
- Project administration: IB TK BK.
- Resources: ZG KC JC TK.
- Software: JPJ ŁP.
- Supervision: IB TK.
- Validation: JPJ ŁP JC KC.
- Writing – original draft: IB TK JPJ ŁP.
- Writing – review & editing: IB JPJ ŁP JC KC BK AK NS ZG TK.
References
- 1. Nolet BA, Rosell F. Comeback of the beaver Castor fiber: an overview of old and new conservation problems. Biol Conserv. 1998; 83: 165–173.
- 2. Jones CG, Lawton JH, Shachak M. Organisms as ecosystem engineers. Oikos. 1994; 69: 373–386.
- 3. Rosell F, Bozsér O, Collen P, Parker H. Ecological impact of beavers Castor fiber and Castor canadensis and their ability to modify ecosystems. Mamm Rev. 2005; 35(3–4): 248–276.
- 4. Ciechanowski M, Kubic W, Rynkiewicz A, Zwolicki A. Reintroduction of beavers Castor fiber may improve habitat quality for vespertilionid bats foraging in small river valleys. Eur J Wildl Res. 2011; 57: 737–747.
- 5. Naiman RJ, Johnston CA, Kelly JC. Alteration of North American streams by beaver. BioScience. 1988; 38: 753–762.
- 6. Collen P, Gibson RJ. The general ecology of beavers (Castor spp.), as related to their influence on stream ecosystems and riparian habitats, and the subsequent effects on fish. Rev. Fish Biol. Fish. 2001; 10: 439–461.
- 7. Obidzinski A, Orczewska A, Cieloszczyk P. The impact of beavers' (Castor fiber L.) lodges on vascular plant species diversity in forest landscape. Pol J Ecol. 2011; 59(1): 63–73.
- 8. Wilsson L. Observations and experiments on the ethology of the European beaver (Castor fiber L.). Viltrevy. 1971; 8: 115–260.
- 9. Zurowski W. Building activity of beavers. Acta Theriol. 1992; 37: 403–411.
- 10. Schlatt S, De Geyter M, Kliesc S, Nieschlag E, Bergmann M. Spontaneous recrudescence of spermatogenesis in the photoinhibited male Djungarian hamster, Phodopussungorus. Biol Reprod. 1995; 53: 1169–1177. pmid:8527522
- 11. Prendergast BJ. Internalization of seasonal time. Horm Behav. 2005; 48: 503–511. pmid:16026787
- 12. Johnson L, Thompson DL. Age-related and seasonal variation in the Sertoli cell population, daily sperm production and serum concentrations of follicle-stimulating hormone, luteinizing hormone and testosterone in stallions. Biol Reprod. 1983; 29(3): 777–789. pmid:6414546
- 13.
Glover TD, D’Occhio MJ, Millar RP. Male life cycle and seasonality. Marshall’s Physiology of Reproduction Churchill Livingstone, Edinburgh, 1990.
- 14. Zamiri MJ, Khalili B, Jafaroghli M, Farshad A. Seasonal variation in seminal parameters testicular size, and plasma testosterone concentration in Iranian Moghani rams. Anim Reprod Sci. 2010; 94: 132–136.
- 15. Li Q, Zhang F, Zhang S, Sheng X Han Y, Yuan Z, et al. Seasonal expression of androgen receptor, aromatase, and estrogen receptor alpha and beta in the testis of the wild ground squirrel (Citellusdauricus Brandt). Eur J Histochem. 2015; 59: 9–16. pmid:25820559
- 16. Chojnowska K, Czerwinska J, Kaminski T, Kaminska B, Panasiewicz G, Kurzynska A, et. al. Sex- and seasonally related changes in plasma gonadotropins and sex steroids concentration in the European beaver (Castor fiber). Eur J Wildl Res. 2015; 61: 807–811.
- 17. Czerwinska J, Chojnowska K, Kaminski T, Bogacka I, Panasiewicz G, Smolinska N, et al. Plasma-Glucocorticoids and ACTH Levels During Different Periods of Activity in the European Beaver (Castor fiber L.). Folia Biol (Krakow). 2015; 63(4): 229–234.
- 18. Criscuolo A, Brisse S. Alien Trimmer: a tool to quickly and accurately trim off multiple short contaminant sequences from high-throughput sequencing reads. Genomics. 2013; 102: 500–506. pmid:23912058
- 19. Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, et al. De novo transcript sequence reconstruction from RNA-Seq: reference generation and analysis with Trinity. Nat Protoc. 2013; 8: 1494–1512. pmid:23845962
- 20. Fu L, Niu B, Zhu Z, Wu S, Li W. CD-HIT: accelerated for clustering the next-generation sequencing data. Bioinformatics. 2012; 28: 3150–3152. pmid:23060610
- 21. Simão FA, Waterhouse RM, Ioannidis P, Kriventseva EV, Zdobnov EM. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics. 2015; 31(19): 3210–3212. pmid:26059717
- 22. Moghadam HK, Harrison PW, Zachar G, Szekely T, Mank JE. The plover neurotranscriptome assembly: transcriptomic analysis in an ecological model species without a reference genome. Mol Ecol Resour. 2013; 13: 696–705. pmid:23551815
- 23. Chen Y-C, Harrison PW, Kotrschal A, Kolm N, Mank JE, Panula P. Expression change in Angiopoietin-1 underlies change in relative brain size in fish. Proc Biol Sci. 2015; 7; 282(1810). pii: 20150872. pmid:26108626
- 24. Conesa A, Götz S, García-Gómez JM, Terol J, Talon M, Robles M. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005; 21: 3674–3676. pmid:16081474
- 25. Jones P, Binns D, Chang HY, Fraser M, Li W, McAnulla C et al. InterProScan 5: genome-scale protein function classification. Bioinformatics. 2014; 30: 1236–1240. pmid:24451626
- 26. Wu TD, Watanabe CK. GMAP: a genomic mapping and alignment program for mRNA and EST sequences. Bioinformatics. 2005; 21: 1859–1875. pmid:15728110
- 27. Krzywinski M, Schein J, Birol I, Connors J, Gascoyne R, Horsman D. Circos: an information aesthetic for comparative genomics. Genome Res. 2009; 19: 1639–1645. pmid:19541911
- 28. Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009; 10,R25–R25.10. pmid:19261174
- 29. Roberts A, Pachter L. Streaming fragment assignment for real-time analysis of sequencing experiments. Nat Methods. 2012; 10: 71–73. pmid:23160280
- 30. Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010; 11, R106. pmid:20979621
- 31. Robinson MD, McCarthy DJ, Smyth GK. Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010; 26: 139–140. pmid:19910308
- 32. Chojnowska K, Czerwinska J, Kaminski T, Kaminska B, Kurzynska A, Bogacka I. Leptin plasma concentrations, leptin gene expression and protein localisation in the hypothalamic-pituitary-gonadal/adrenal axes of the European beaver (Castor fiber). Theriogenology. 2017; 87: 266–275. pmid:27780608
- 33. Czerwinska J, Chojnowska K, Kaminski T, Bogacka I, Smolinska N, Kaminska B. Orexin receptor expression in the hypothalamic-pituitary-adrenal and hypothalamic-pituitary-gonadal axes of free-living European beavers (Castor fiber L.) in different periods of the reproductive cycle. Gen Comp Endocrinol. 2017: 240: 103–113. pmid:27664717
- 34. Schmittgen TD, Livak KJ. Analyzing real-time PCR data by the comparative C(T) method. Nat Protoc. 2008; 3: 1101–1108. pmid:18546601
- 35. Marioni JC, Mason CE, Mane SM, Stephens M, Gilad Y. RNA-seq: an assessment of technical reproducibility and comparison with gene expression arrays. Genome Res. 2008; 18(9): 1509–1517. pmid:18550803
- 36. Gunawan A, Sahadevan S, Neuhoff C, Große-Brinkhaus C, Gad A, Frieden L, et al. RNA deep sequencing reveals novel candidate genes and polymorphisms in boar testis and liver tissues with divergent androstenone levels. PLoS One. 2013; 8(5), e63259. pmid:23696805
- 37. Nie Q, Fang M, Jia X, Zhang W, Zhou X, Xiaomei He, et al. Analysis of muscle and ovary transcriptome of Sus scrofa: assembly, annotation and marker discovery. DNA Res. 2011; 18(5): 343–351. pmid:21729922
- 38. Liu H, Wang T, Wang J, Quan F, Zhang Y. Characterization of Liaoning cashmere goat transcriptome: sequencing, de novo assembly, functional annotation and comparative analysis. PLoS One. 2013; 8(10), e77062. pmid:24130835
- 39. Salem M, Paneru B, Al-Tobasei R, Abdouni F, Thorgaard GH, Rexroad CE, et al. Transcriptome assembly, gene annotation and tissue gene expression atlas of the rainbow trout. PLoS One. 2015; 10(3), e0121778. pmid:25793877
- 40. Deng T, Pang C, Lu X, Zhu P, Duan A, Tan Z, et al. De Novo Transcriptome Assembly of the Chinese Swamp Buffalo by RNA Sequencing and SSR Marker Discovery. PLoS One. 2016; 11(1), e0147132. pmid:26766209
- 41. Petrusz P, Jeyaraj DA, Grossman G. Microarray analysis of androgen-regulated gene expression in testis: the use of the androgen-binding protein (ABP)-transgenic mouse as a model. Reprod Biol Endocrinol. 2005; 3,70. pmid:16336681
- 42. Sadate-Ngatchou PI, Pouchnik DJ, Griswold MD. Identification of testosterone-regulated genes in testes of hypogonadal mice using oligonucleotide microarray. Mol Endocrinol. 2004; 18(2): 422–433. pmid:14605096
- 43. Tacken PJ, Van der Zee A, Beumer TL, Florijn RJ, Gijpels MJJ, Havekes LM, et al. Effective generation of very low density lipoprotein receptor transgenic mice by overlapping genomic DNA fragments: high testis expression and disturbed spermatogenesis. Transgenic Res. 2001; 10(3): 211–221. pmid:11437278
- 44. Mendis SH, Meachem SJ, Sarraj MA, Loveland KL. Activin A balances Sertoli and germ cell proliferation in the fetal mouse testis. Biol Reprod. 2011; 84(2): 379–391. pmid:20926807
- 45. Barakat B, Itman C, Mendis SH, Loveland KL. Activins and inhibins in mammalian testis development: new models, new insights. Mol Cell Endocrinol. 2012; 359(1–2): 66–77. pmid:22406273
- 46. Itman C, Loveland KL. SMAD expression in the testis: an insight into BMP regulation of spermatogenesis. Dev Dyn. 2008; 237(1): 97–111. pmid:18069690
- 47. Cérec V, Piquet-Pellorce C, Aly HA, Touzalin AM, Jégou B, Bauché F. Multiple pathways for cationic amino acid transport in rat seminiferous tubule cells. Biol Reprod. 2007; 76(2): 241–249. pmid:17065601
- 48. Matzkin ME, Lauf S, Spinnler K, Rossi SP, Köhn FM, Kunz L, et al. The Ca2+-activated, large conductance K+-channel (BKCa) is a player in the LH/hCG signaling cascade in testicular Leydig cells. Mol Cell Endocrinol. 2013; 367(1–2): 41–49. pmid:23267835
- 49. Gong XD, Li JC, Leung GP, Cheung KH, Wong PY. A BK(Ca) to K(v) switch during spermatogenesis in the rat seminiferous tubules. Biol Reprod. 2002; 67(1): 46–54. pmid:12079998
- 50. Wu WL, So SC, Sun YP, Zhou TS, Yu Y, Chung YW, et al. Functional expression of a Ca2+-activated K+ channel in Xenopus oocytes injected with RNAs from the rat testis. Biochim Biophys Acta. 1998; 1373(2): 360–365. pmid:9733997
- 51. Loveland K, Schlatt S, Sasaki T, Chu ML, Timpl R, Dziadek M. Developmental changes in the basement membrane of the normal and hypothyroid postnatal rat testis: segmental localization of fibulin-2 and fibronectin. Biol Reprod. 1998; 58(5): 1123–1130. pmid:9603244
- 52. Glattauer V, Irving-Rodgers HF, Rodgers RJ, Stockwell S, Brownlee AG, Werkmeister JA et al. Examination of basement membrane components associated with the bovine seminiferous tubule basal lamina. Reprod Fertil Dev. 2007; 19(3): 473–481. pmid:17394796
- 53. Naillat F, Veikkolainen V, Miinalainen I, Sipilä P, Poutanen M, Elenius K. ErbB4, a receptor tyrosine kinase, coordinates organization of the seminiferous tubules in the developing testis. Mol Endocrinol. 2014; 28(9): 1534–1546. pmid:25058600
- 54. Conner JR, Smirnova II, Poltorak A. Forward genetic analysis of Toll-like receptor responses in wild-derived mice reveals a novel antiinflammatory role for IRAK1BP1. J Exp Med. 2008; 205(2): 305–314. pmid:18268037
- 55. Conner JR, Smirnova II, Moseman AP, Poltorak A. IRAK1BP1 inhibits inflammation by promoting nuclear translocation of NF-kappaB p50. Proc Natl Acad Sci USA. 2010; 107(25): 11477–11482. pmid:20534545
- 56. Dzau VJ, Ellison KE, Brody T, Ingelfinger J, Pratt RE. A comparative study of the distributions of renin and angiotensinogen messenger ribonucleic acids in rat and mouse tissues. Endocrinology. 1987; 120(6): 2334–2338. pmid:3552634
- 57. Hagaman JR, Moyer JS, Bachman ES, Sibony M, Magyar PL, Welch JE, et al. Angiotensin-converting enzyme and male fertility. Proc Natl Acad Sci USA. 1988; 95(5): 2552–2557.