Skip to main content
Download PDF
  • Research article
  • Open access
  • Published:

Dynamic expression and functional analysis of circRNA in granulosa cells during follicular development in chicken

BMC Genomics volume 20, Article number: 96 (2019) Cite this article

Abstract

Background

Circular RNA (circRNA) is a type of noncoding RNA involved in a variety of biological processes, especially in post-transcriptional regulation. The granulosa cells of follicles play a determining role in ovarian development. However, the function of circRNA in chicken follicles is unclear. To better understand the molecular mechanism underlying follicular development and granulosa cell function, we performed a strategy of second-generation sequencing and linear RNA depletion for granulosa cells from small yellow follicles (SYF, 5–8 mm), the smallest hierarchal follicles (F6, 9–12 mm), and the largest hierarchal follicles (F1, ~ 40 mm).

Results

We predicted a total of 11,642 circRNAs that distributed on almost all chromosomes. The majority of the splice lengths of circRNAs were 200–500 nt and mainly produced from intron and CDS regions. During follicle growth, differentially expressed (DE) circRNAs showed dynamic changes which were tissue- and stage-specific. The host genes of DE circRNAs were functionally enriched in GTPase activity and several pathways involved in reproduction. Moreover, bioinformatic prediction analysis for circRalGPS2 demonstrated that circRNAs from the same genes may share common miRNA to act as a sponge. The predicted target genes were enriched in various biological processes including cognition, cell communication, and regulation of signaling, and several pathways related to reproduction such as tight junction, oocyte meiosis, progesterone-mediated oocyte maturation, and GnRH signaling.

Conclusions

This study provides a starting point for further experimental investigations into chicken circRNAs and casts a light on the understanding of follicle development.

Background

Follicle development is a key determinant of reproduction performance and is affected by physiology and the environment. Despite a different follicle situation between birds and mammals, the characteristics of maturing hen ovarian follicles are equivalent to those occurring in mammals, thus hen ovarian follicles can be used as a model for the study of follicle development [1]. Ovulation is a process that occurs via the initial recruitment of primordial follicles, which is followed by primary follicles, small white follicles, larger white follicles, small yellow follicles (SYFs), and hierarchal follicles that progress toward maturity [2]. One of the striking characteristics of the bird follicle is the presence of 5–7 hierarchal follicles before ovulation. Hierarchal follicles develop from a pool of SYFs each day during the peak laying period [3].

The basic structure of the follicular wall consists of the theca and granulosa. The ovaries can increase in follicle number when treated by exogenous follicle-stimulating hormone (FSH) [4,5,6]; considering that FSH receptor (FSHR) mRNA in the theca layer changes little during follicle development [6], it is concluded that the granulosa cell (GC) contributes as a decisive factor in follicle development. The development of GCs involves multiple complex biological processes, affected by many molecular components including those at the genome and RNA transcript level [7, 8]. The expression of AMH mRNA increases with follicle development, and markedly decreases during follicle selection [2] and may affect SYF number [9]. It has been shown that epidermal growth factor (EGF) suppresses GC differentiation and decreases the abundance of FSHR mRNA [10]. FOXL2 plays different roles in prehierarchical GCs and pre-ovulatory GCs [11]. In brief, the growth, differentiation, and function of GCs are dependent on FSH and related functional genes. Elucidating the features of the molecular mechanism during GC development could provide a reference to understanding molecular development and improving assisted reproduction in humans.

There are two types of cellular RNA during transcription processes: coding and noncoding. In recent years, more and more noncoding RNA investigations have attracted researchers focused in this area. Notably, circular RNA (circRNA) is becoming a hot topic in life science since Salzman identified ~ 80 circRNAs in 2012 [12]. The circRNA exhibit characteristics without a 5′-cap structure or 3′-polyadenylated tail, mainly derived from exons, introns, intergenic, and antisense transcripts and are widespread in mammals and plants [13, 14]. The circRNA are recognized as having the following biological functions: an miRNA sponge, transcription template, regulating gene transcription, and RNA-binding proteins [15, 16]. In particular, the expression patterns of circRNAs are often tissue- and stage-specific [17], and highly conserved across the human brain and other species [18]. With the fast development of bioinformatics techniques and deep sequencing, a number of circRNAs has been detected in animal follicle development. Tao et al. reports that the host genes of the circRNAs in goat pre-ovulatory ovarian follicles are involved in the ovarian steroidogenesis pathway and the p53 signaling pathway [19], chi-circ_0008219 modulates follicle growth by sponging three miRNAs to achieve. circRNA_103827 and circRNA_104816 detected in human GCs could be a potential indicator of a compromised follicular microenvironment [20]. However, data associated with the function of circRNAs in chicken is still scarce, especially with respect to follicle development. This led us to systematically study the expression of circRNAs in chicken follicle GCs.

In the present study, we performed RNA-seq for circRNAs during chicken follicle development to explore the structure and expression profile of GCs.

Methods

Granulosa cells obtained for circular RNA-seq

Sixteen generations of Jinghai Yellow Chicken were reared on Jiangsu Jinghai Poultry Industry Group Co., Ltd. (Nantong City, Jiangsu Province, China). Hens were caged individually while being transferred to the laying house, with a 16 L:8D light regime and ad libitum access to water and a restricted diet. Body weight was measured, and based on the pedigree record, three half-sib hens with an average body weight at 27 weeks of age were humanely sacrificed using 60–70% carbon dioxide. Only hens with a soft eggshell egg in the oviduct were considered in our study. Follicles including SYF (5–8 mm in diameter), hierarchal follicles from smallest to largest including those from F6 to F1 (smallest follicles were 9–12 mm in diameter, largest follicles were ~ 40 mm in diameter), large white follicles (LWF), stroma ovary (O) and postovulatory (Po) follicle, uterus (Ut) and oviduct (Ov) were collected immediately. Detailed methods for collecting follicle GCs and theca cells can be found in a previous paper [21]. We removed the yolk from follicle walls carefully and separated GCs and theca cells, and rinsed them with PBS. The GCs were frozen in liquid nitrogen and used for RNA-seq analysis and experimental validation.

RNA-seq preparation

Total RNA was isolated by TRIzol reagent (Invitrogen, USA). A series of methods and experimental instruments including 1% agarose gels, Qubit® RNA Assay Kit in Qubit® 2.0 Fluorometer (Life Technologies, USA), the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, USA) and NanoPhotometer® spectrophotometer (IMPLEN, USA) were used to assess the contamination and degradation, concentration, integrity and purity of total RNA.

A total of 5 μg RNA from each sample was used to construct the circRNA library. First, ribosomal RNA was removed using an Epicentre Ribozero™ rRNA Removal Kit (Epicentre, USA). Second, the linear RNA was digested with 3 U of RNase R (Epicentre) per μg of RNA. The sequencing libraries were generated by NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina® (NEB, USA) following manufacturer’s recommendations. Briefly, first strand cDNA was synthesized using a random hexamer primer and M-MuLV Reverse Transcriptase (RNase H). Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. After adenylation of the 3′ ends of DNA fragments, NEBNext Adaptor with hairpin loop structures were ligated to prepare for hybridization. After PCR reaction for fragmentation, and adenylation, products were purified (AMPure XP system) and the library quality was assessed on an Agilent Bioanalyzer 2100 system. Before library preparations were sequenced on the Illumina Hiseq 4000 platform in Novogene (Beijing, China), the cluster generation using TruSeq PE Cluster Kit v3-cBot-HS (Illumina) was performed. Then, 150 bp paired-end reads were generated.

The circRNA identification and differential expression analysis

Before circRNA identification, quality control was carried out by calculating Q20, Q30, and GC content. Subsequently, paired-end clean reads were aligned to the reference genome using Bowtie [22]. Then, the circRNA were identified by find_circ [23] and CIRI2 [24]. The Circos figure was drawn using Circos software [25].

The raw counts were first normalized using TPM. circRNAs exhibiting fold changes ≥2 with P-adjust ≤ 0.05 were classified as circRNAs with significant differential expression.

Bioinformatics functions analysis

Gene ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were implemented using clusterProfiler package [26] in R software for the host genes of DE circRNA and predicted target genes. MicroRNA target sites in exons of circRNA loci were identified using miRanda (http://www.microrna.org/microrna/home.do) with a score of 140 or higher and free energy − 10 or lower. Target genes of miRNAs were predicted by TargetScan [27] and miRDB [28]. The circRNA-miRNA network was constructed according to the prediction of miRNA binding sites. Cytoscape software was used to construct circRNA-miRNA networks [29].

Validation of circRNAs by sanger sequencing and qRT-PCR

The circRNAs were validated with convergent and divergent primers according to a previous study [18]. Details of divergent and convergent primers are in Additional file 1: Table S1. PCR products of divergent and convergent primers for cDNA and genomic DNA were analyzed by agarose gel electrophoresis. Back-splicing sites of circRNAs were confirmed by Sanger sequencing at Sango Biotech Co. Ltd. (Shanghai, China).

The expression of PCR products from divergent primers for each GC were validated using qRT-PCR. The qRT-PCR program was implemented using ABI7500 (Life Technologies, USA), with the SYBR Green (TaKaRa, cat. #RR820A) method in a final volume of 20 μL. Each assay was employed in triplicate using the following cycling conditions: 95 °C/30 s and 40 cycles of 95 °C/5 s and 60 °C/34 s. The △△Ct method was used to compare gene expression, with β-actin as a reference gene.

In addition, RNA samples from hierarchal follicles from the F2 to F5 GCs and theca cells, LWF, O, Po, Ut, and Ov were analyzed by qRT-PCR. The expression of circRNAs at different ovary tissues was calculated via the 2- ΔΔ Ct method.

Results

Overview of circRNAs

A total of 864 million raw reads were obtained from RNA-seq data Clean data obtained by removing adapter and low quality sequences were mapped to the chicken reference genome (Gallus-gallus-5.0/galGal5). Table 1 summarizes the total number of reads generated from nine samples, and each sample yielded more than 72 million raw reads data. The average GC content was 60.34% (Table 1). Circos figure for circRNAs in GCs is displayed in Fig. 1, which showed that the numbers of circRNAs in GCs were distributed almost on all chromosome, with the largest on chromosome 1 (chr1) and the smallest on the chr16.

Table 1 Summary of sequencing results in chicken follicle granulosa cells at three stages
Full size table
Fig. 1
figure 1

Overview of sequencing data of circRNA by Circos plots. Note: Plot inside are chromosomes from 1 to 28 and sex chromosome converged across granulosa cells (GCs) from SYF (G1, G2, G3), F6 (G1, G2, G3) to F1 (G1, G2, G3), SYF: small yellow follicle; F6: the smallest hierarchal follicles; F1: the largest hierarchal follicle; G: granulosa cell

Full size image

After ribosomal-depleted RNA analysis, a total of 11,642 circRNAs were detected by find_circ and CIRI2 (Additional file 1: Table S2). The length of most circRNAs was about 200–300 nt (Fig. 2a) with genomic distances of ~ 1000–3000 bp (Fig. 2b). We grouped circRNAs based on their genomic region including intron, CDS, 5′UTR, intergenic, 3′UTR, and lincRNA. After this step was performed, circRNAs were mainly produced from intron and CDS regions (Fig. 2c). Multiple circRNAs were produced by one gene and the majority of circRNAs have two exons, followed by three exons and one exon (Fig. 2d). Among these CDS sourced circRNAs, about 2371 host gene lengths were more than 8000 nt (Fig. 2e). Approximately, 9165 and 9149 flank introns (42.17%) of circRNAs were above 104~105nt, respectively (Fig. 2f, Additional file 1: Table S3).

Fig. 2
figure 2

Profiling of circular RNAs in chicken follicle granulosa cells. Note: (a), splice length of circRNAs; (b), genomic distance between circRNAs; (c), source of circRNAs; (d), exon numbers; (e), length of host genes; (f), flank intron length distributed in granulosa cells

Full size image

Differentially expressed analysis

From the nine GC samples, a total of 6337 circRNAs were in SYF, 8189 in F6, 7323 in F1, and 3742 prevailed in all three follicle GCs (Fig. 3a, Additional file 1: Table S4). Volcano plots of circRNAs expressed differentially in GCs are displayed in Fig. 3b, c, d. In a comparison of SYF and F6 GC libraries, 219 circRNAs were significantly differentially expressed, including 77 upregulated and 142 downregulated (Fig. 3b), while 687 circRNAs were significantly differentially expressed between the SYF and F1 GC libraries (Fig. 3c). In the hierarchal follicle GC libraries, a total of 412 circRNAs with 191 upregulated and 221 downregulated in F6 GCs were significantly differentially expressed (Fig. 3d).

Fig. 3
figure 3

Features of differential expression circRNAs. Note: (a), Venn analysis of differentially expressed circRNAs in three stages. b, volcano map of differentially expressed circRNAs (DEcircs) between SYF and F6 granulosa cells. c, volcano map of DEcircs between SYF and F1 granulosa cells. d, volcano map of DEcircs between F6 and F1 granulosa cells. SYF: small yellow follicle; F6: the smallest hierarchal follicles; F1: the largest hierarchal follicle; G: granulosa cell

Full size image

A heatmap of all differentially expressed circRNAs based on transcript per million (TPM) values is shown in Fig. 4. Samples at the same stages are clustered together, the F6_GCs and SYF_GCs are clustered together. The circRNAs changed their expression levels according to the heatmap. In particular, the expression between F1_GCs and SYF_GCs showed the opposite tendency.

Fig. 4
figure 4

Heat map representing the relative expression levels of all differentially expressed circRNAs in the three stages (Fold change > 2, P < 0.05). Note: SYF: small yellow follicle; F6: the smallest hierarchal follicles; F1: the largest hierarchal follicle; G: granulosa cell.

Full size image

Functional analyses of host genes of differentially expressed circRNAs

The terms from GO of all differentially expressed circRNAs (DEcirc) host genes included GTPase related process and cell development (Additional file 1: Table S5, Fig. 5). The KEGG pathway enrichment of all DEcirc host genes displayed several significant pathways (P-adj < 0.05), including the Mucin type O-glycan biosynthesis (4.35%), MAPK signaling pathway (13.04%), Adherens junction (6.09%), FoxO signaling pathway (8.70%), regulation of actin cytoskeleton (10.43%), and progesterone-mediated oocyte maturation (6.09%) (Additional file 1: Table S6).

Fig. 5
figure 5

GO and KEGG analysis for the host genes of all differentially expressed circRNAs. Note: (a), The top 10 GO enrichment term in BP, CC, and MF for the host genes of all differentially expressed circRNAs; (b), the top 20 KEGG enrichment term for the host genes of all differentially expressed circRNAs

Full size image

The DE circRNAs between three stages according to Venn analysis are shown in Fig. 6. Notably, 36 pervasive circRNAs existed in three types of GC (Fig. 6a). We clustered the 36 circRNAs based on the relative expression value (Fig. 6b), and the deciphered heat map showed that the expression pattern between SYF and F1 GCs displayed the opposite trend.

Fig. 6
figure 6

Differentially expressed circRNAs in three different stages of granulosa cells. Note: a, Venn diagram of differentially expressed circRNAs (n = 891, fold change > 2, q-value < 0.05); b, Heat map of pervasive circRNAs shows expression value at three stages. SYF: small yellow follicle; F6: the smallest hierarchal follicles; F1: the largest hierarchal follicle; G: granulosa cells.

Full size image

Validation circRNAs and expression pattern

We randomly confirmed back-site junctions of all chosen circRNAs by agarose gel electrophoresis, Sanger sequencing, and quantitative real-time PCR (qRT-PCR) (Fig. 7). Divergent and convergent primers were designed for each circRNAs, and both genomic DNA (gDNA) and complementary DNA (cDNA) were used as templates for PCR.

Fig. 7
figure 7

Validation of circRNAs by experimental and sequencing. Note: (a), agarose gel electrophoresis test for divergent primers and convergent primers amplify circRNAs; The brightest DNA marker is 200 bp. cc: convergent primers for cDNA; cg: convergent primers for genomic DNA; dc: divergent primers for cDNA; dg: divergent primers for genomic DNA. b, Sanger sequencing confirmed the back-splicing junction of circRNAs; (c), RT-qPCR validation of nine differentially expressed circRNAs in three types of granulosa cell. TPM: Transcripts Per Million; RQ: relative expression; Graph number of a, b, c, d, e, f, g, h, i represent circRNAs: circRNA_8: 6369673|6,402,248, circRNA_8: 6369673|6,422,097, circRNA_8: 6384248|6,402,248, circRNA_26: 4587446|4,664,529, circRNA_3: 66399375|66,402,723, circRNA_20: 2426006|2,445,902, circRNA_14: 84666291|8,458,456, circRNA_12: 9000575|9,001,598, circRNA_23: 5084518|5,086,817, respectively. SYF: small yellow follicle; F6: the smallest hierarchal follicles; F1: the largest hierarchal follicle; G: granulosa cell.

Full size image

Divergent primers from each circRNA amplified the expected fragments by cDNA as a template, while there were no PCR products for the gDNA template, which suggested the presence of back-site junctions (Fig. 7a), which were validated by Sanger sequencing (Fig. 7b). The qRT-PCR validation for the nine circRNAs was consistent with the trends obtained from circRNA sequencing data (Fig. 7c). Notably, three circRNAs, chr8:6369673|6,402,248 (circRalGPS2_1), chr8:6369673|6,422,097 (circRalGPS2_2), and chr8:6384248|6,402,248 (circRalGPS2_3), were from the same gene, RalGPS2, and showed pervasive differential expression in three types of GCs.

Moreover, the expression pattern of circRalGPS2 in different follicle cells were detected by qRT-PCR, and results are shown in Fig. 8. These three types of circRNAs were all detected in all reproduction tissues. Interestingly, the expression level in the ovary stroma was the highest. The expression level was lower in the hierarchal follicle, GCs, or theca cells compared to the prehierarchal follicle and reproduction tact.

Fig. 8
figure 8

Expression pattern of circRLAGPS2 during ovarian development

Note: F1, F2, F3, F4, F5, F6 represent hierarchal follicles from larger to small; SYF, small yellow follicle; O, ovary; LWF, large white follicle, POF, post-ovulation follicle; Ut, uterus; Ov, oviduct.

Full size image

Target miRNA or gene prediction, networks, and pathways

As demonstrated by a previous study regarding circRNA main functions as a microRNA sponge, we predicted the 36 circRNAs by miRanda based on its sequence. The top five predicted microRNAs are displayed in Table 2. It was observed that five miRNAs (miR-1625-3p, miR-1552-3p, miR-16-2-3p, miR-18b-3p, and miR-200a-3p) were common among the three types of circRNAs sourced from the same gene RalGPS2 (Fig. 9). A number of overlap genes (538) were predicted by TargetScan and miRDB for these five miRNAs (Additional file 1: Table S6). GO and KEGG analysis was performed for the putative overlap target genes to identify the potential pathways (Additional file 1: Table S8 and S9).

Table 2 Thirty-six pervasive DE circRNAs with miRNA binding sites
Full size table
Fig. 9
figure 9

Predicted biomathematical circRNA–miRNA network for the same soured circRNAs. Note: Pink square represents circRNAs, blue circle represents miRNAs.

Full size image

The results showed that a total of 145 predicted target genes were annotated to 92 pathways (Additional file 1: Table S9). Among these, pathways involving tight junctions and oocyte meiosis were significantly enriched (P-adj < 0.05). Other pathways involving in GnRH signaling, Wnt signaling, MAPK signaling pathway, progesterone-mediated oocyte maturation, TGF-beta signaling pathway, and FoxO signaling were related to animal reproduction. The enriched genes related to reproduction pathways included MAP3K5, ADCY3, PLCB1, CAMK2G, PIK3R1, and another 37genes (Additional file 1: Table S8).

Discussion

The field of chicken reproductive function biology is still in its infancy and many questions remain. If hens can produce longer layer periods, they could be taken full advantage of and economic benefits would improve. Investigations into the GCs of follicles could provide information underlying follicle development and make further improvements in egg production. Moreover, other applications could include a model to study the molecular events that are the most common fate of ovarian follicles [30]. Folliculogenesis in birds is a very complex process, which includes recruitment of the primitive follicle, development or atresia of small follicles and 6–8 mm follicle selection, hierarchal follicle rapid growth, and ovulation. Knowledge of developmental changes during the folliculogenesis process in chicken can be incorporated into egg production and has important ramifications for understanding ovarian function with chicken as a model organism [30, 31]. The circRNA and polyA-RNAs have rich and extremely important functions on a variety of biological processes [23]. Though many studies have reported that circRNAs were widespread in mammals [32] and plants [14], little is known about circRNAs in chickens, especially in the folliculogenesis process.

In the current study, a strategy of second generation sequencing for circRNAs in chicken follicle GCs was performed. The circRNAs were widely distributed in chromosomes 1–28, 30, 32–33 Z, W, and two linkage groups (Fig. 1, Additional file 1: Table S1), which suggested that circRNAs had complexity and functional diversity. The circRNAs distributed on chr1 were the largest, while those on chr16 were the smallest, which was in line with the chromosome length. GGA1 is the longest chromosome in the entire genome of the chicken and chr16 is the shortest [33]. We predicted a total of 11,642 circRNAs in three types of follicle GC, a number which is less than that detected in the muscle of embryonic chicken (13,377) [34] and more than that found in the liver of ALVJ-resistant and susceptible chicken (1800) [35], which suggested that the expression of circRNAs were tissue-specific and assisted in the hypothesis that the circularization is a tightly regulated event in different tissues [32]. In our study, the spliced length of the majority of circRNAs was 200~300 nt, which was longer than those obtained in a previous study from pre-implantation embryos (124–227 nt) [36]. In total, 94.83% of spliced circRNAs were less than 500 nt, which was higher than in other studies of circRNA expression in rat liver (45.56%) [37]. The difference of length in different animals may be due to a species difference in the function performed by circRNAs. With respect to the source of circRNAs in GCs, it was shown that circRNAs mainly came from introns, which was unlike in previous research where most circRNAs came from coding exons [34]. Other aspects about the production rate or exon number in circRNAs were similar with previous studies, for instance, the production rate of circRNAs is determined by intronic sequence [38] and a single gene can produce multiple circRNAs [39]. When all the above features on circRNA are combined, it was concluded that not only was there a species difference in circRNA expression, which was tissue-specific, but common features were also shared.

With respect to differential expression in circRNAs, Fig. 3 shows that GC circRNAs exhibited stage-specific expression, which is in line with reports in human and mouse [32]. The most abundant expression was in F6_GCs with 2095 stage-specific expression circRNAs, followed by F1_GCs. The fluctuation of circRNA abundance in GCs may be related to their specific roles during the development of follicles.

We drew a heat map for all differential expression circRNAs according to its TPM (Fig. 4), which indicated that these circRNAs have a dynamic expression pattern during GC growth. Notably, the expression level between F6_GCs and SYF_GCs showed an opposite tendency, and these two GCs have different cell morphologies during follicle development; small follicles contain several layers of GCs whereas mature follicles have a single layer of GCs [30]. It was concluded that circRNAs may exert different functions on different cell morphology of the same tissue.

Moreover, GO and KEGG pathway analyses for these differentially expressed circRNAs were performed. Most GO terms were enriched in GTPase related processes and cell developmental processes. GTPases of the Ras superfamily are omnipresent signaling proteins that play important roles in a wide range of vital cellular processes [40], including ovarian follicle development [41, 42]. A study showed that the activation of GTPases is critical for the regulation of cell junctions and cell-cell adhesion [43] and it was shown to be responsible for the accumulation of cell junction proteins between oocytes and GCs [41]. Tight junctions between oocytes and GCs are dynamically regulated in a developmental stage-specific manner [44]. Some GTPase proteins also regulate GC proliferation, differentiation and follicle selection in chicken [42]. Taken together, these observations and our results indicate that those circRNAs are involved in GCs progression and communications between oocyte and GCs. Moreover, several other pathways related to reproduction are also significantly enriched including pathways involved in MAPK signaling, FoxO signaling, and progesterone-mediated oocyte maturation. Similar results were observed in noncoding RNAs including long non-coding RNA and microRNA enrichment analysis in previous studies [45, 46], which suggest that these circRNAs may act as competing endogenous RNAs (ceRNAs) to modulate gene transcription. Findings obtained from our study and previous reports suggest that follicle development is a complex process that includes numerous events, and circRNAs distributed on three types of follicle GCs are mainly involved in tight junctions to participate in follicle growth, reproduction events, and other biological processes.

As demonstrated in Fig. 6, a total of 36 differentially expressed circRNAs showed pervasive expression patterns in three types of GCs. Interestingly, expression levels presented an increased or decreased trend with the development of follicle GCs. Previous work has reported that the abundance of circRNAs were negatively correlated with proliferation in cancer tissue compared to normal tissue [47] and circularization of circRNAs can compete with canonical pre-mRNA splicing [38]. The correlation of circRNA expression level and the quality of embryo numbers were reported to be negative in Quan et al. [48]. It is tempting to speculate that these upregulated circRNAs may inhibit GC proliferation and downregulated circRNAs may affect follicle number by participating in the regulation of follicle recruitment.

We observed that three circRNA isoforms (circRalGPS2_1, circRalGPS2_2, and circRalGPS2_3) from the same host gene RalGPS2 (known as a hot spot gene [49]) widely distributed on three GCs and displayed similar differential expression trend, which is unlike previous reports where they displayed dissimilar expression profiles [49]. Presumably, these circRNAs are generated from one particular exon that is common to RalGPS2. It remains to be validated whether circRNAs from the same back-splicing can perform a similar expression level to modulate the transcription process. In addition, previous research showed that RalGPS2 may play a role in cytoskeleton reorganization [50] and is responsible for survival and the cell cycle in lung cancer cells [51]. To date, few data have been published regarding the role of RalGPS2 in follicle GCs, whether it is GC-type-specific as with mRNA. Interestingly, RalGPS2 was enriched in many GTPase processes in our study (Additional file 1: Table S4). We examined the expression changes of circRalGPS2s during follicle growth. The results showed that the expression of all circRalGPS2s were predominant in the ovary stroma and a decreased level was found in follicle growth, suggesting that circularization of circRNAs may not accompany GCs development. It has been proven that circRNAs such as circRNA_103827 and circRNA_104816 expression level is negatively related to embryonic development [20]. With the combined analysis of the abovementioned reports and current results, it is plausible to hypothesize that the lower the level of circRNAs, the more abundant are the display of biological activities, and circRNA from different host genes may exert different functions in different tissues, organs, or stages.

The main mechanism of circRNAs may act as a miRNA sponge [16] to modulate post-transcriptional regulation [23]. The higher level of circRNAs, the stronger is the efficiency of the sponge [16]. The circRNA-miRNA network (Fig. 9) displayed that circRNAs generated from the same gene may regulate common target miRNA that bound to the same genes to perform its function. A study about human GCs showed that the predicted target genes of circRNA_103827 and circRNA_104816 are mainly enriched in glucose metabolism, the mitotic cell cycle, and ovarian steroidogenesis [20]; these pathways play a key role in follicle development. GO and KEGG analysis for the predicted target genes of circRalGPS2 showed that the target genes were enriched in the regulation of signaling, secretion by the cell, and cell motility (Fig. 10a), suggesting that circRNAs may be potentially responsible for these functional areas. KEGG bioinformatics analysis results indicated that circRalGPS2s were significantly involved in pathways of tight junctions and oocyte meiosis. Tight junctions were also enriched in DE circRNA annotations, which further explain the role of tight junctions as a key pathway in follicle development. Other pathways related to reproduction were also enriched, such as GnRH signaling, the Wnt signaling pathway, mTOR signaling, MAPK signaling, progesterone-mediated oocyte maturation, and FoxO signaling. Genes enriched in many pathways include ADCY3 (ENSGALG00000016608), PIK3R1 (ENSGALG00000014786), PAK3 (ENSGALG00000008058), and CAMK2G (ENSGALG00000005088), which have been reported to play important roles during reproduction processes. For example, knock-out of ADCY3 leads to the resumption of gamete mitosis [52]. PIK3R1 is involved in the PI3K/AKT pathway that is responsible for chicken GC proliferation [53]. CAMK2A (the paralog gene of CAMK2G), and PAK3 are located in a putative QTL region that is required for increasing broiler follicle number [54]. Although the function of these genes in chicken granulosa development is unclear, previous reports indicate that they play important roles in regulating reproduction. Furthermore, on the basis of the ceRNA hypothesis [55], circRalGPS2 can compete with miRNA to modulate mRNA expression that regulates downstream target genes. Follicle development is a complex process, with the susceptibility of multiple genes or transcriptional factors as well as the interactions between them. Thus, an assessment of the function of circRNAs as a miRNA sponge to modulate gene expression during follicle development needs tentative overexpression and knockdown experimental designs.

Fig. 10
figure 10

GO and KEGG analysis of circRalGPS2-miRNA target genes. Note: (a), Gene ontology of target genes. b, KEGG pathways of target genes.

Full size image

Conclusions

In summary, our study is the first to elucidate spatiotemporal comprehensive circRNA expression profile in follicle GCs of chicken. These results from RNA-seq for circRNAs suggested that circRNA are abundant, with tissue-specific and stage-specific expression in chicken follicle GCs. In particular, we found a hot spot gene RalGPS2 that produced three isoforms of circRNAs which exhibited a similar expression pattern during follicle development. Furthermore, the function of these circRNAs needs to be investigated by experimental study. However, these findings could be helpful to decipher the regulatory mechanisms underlying follicle growth and engineer practical breeding programs.

Abbreviations

BP:

Biological process

CC:

Cell component

cDNA:

Complementary DNA

ceRNA:

Competing endogenous RNA

DE:

Differential expression

F1:

The largest hierarchal follicle of ~ 40 mm in diameter

F6:

The smallest hierarchal follicle of 9–12 mm in diameter

GC:

Granulosa cell

gDNA:

Genomic DNA

GO:

Biological process in gene ontology

KEGG:

Kyoto Encyclopedia of Genes and Genomes

MF:

Molecular function

O:

Ovary

Ov:

Oviduct

Po:

Postovulatory

qRT-PCR:

Quantitative real-time PCR

SYF:

Small yellow follicle

TPM:

Transcripts Per Million

Ut:

Uterus

References

  1. Bahr JM. The chicken ovary as a model of follicular development. Semin Reprod Med. 1991;9(04):352–9.

    Article  Google Scholar 

  2. Johnson PA. Follicle selection in the avian ovary. Reprod Domest Anim. 2012;47(4):283–7.

    Article  Google Scholar 

  3. Hocking PM. Biology of breeding poultry; 2009.

    Book  Google Scholar 

  4. Palmer SS, Bahr JM. Follicle stimulating hormone increases serum oestradiol-17 beta concentrations, number of growing follicles and yolk deposition in aging hens (Gallus gallus domesticus) with decreased egg production. Br Poult Sci. 1992;33(2):403–14.

    Article  CAS  Google Scholar 

  5. Johnson AL, Bridgham JT, Wagner B. Characterization of a chicken luteinizing hormone receptor (cLH-R) complementary deoxyribonucleic acid, and expression of cLH-R messenger ribonucleic acid in the ovary. Biol Reprod. 1996;55(2):304–9.

    Article  CAS  Google Scholar 

  6. Johnson AL. Ovarian follicle selection and granulosa cell differentiation. Poult Sci. 2015;94(4):781–5.

    Article  CAS  Google Scholar 

  7. Zhu G, Mao Y, Zhou W, Jiang Y. Dynamic changes in the follicular transcriptome and promoter DNA methylation pattern of steroidogenic genes in chicken follicles throughout the ovulation cycle. PLoS One. 2015;10(12):e0146028.

    Article  Google Scholar 

  8. Kang L, Yang C, Wu H, Chen Q, Huang L, Li X, Tang H, Jiang Y. miR-26a-5p Regulates TNRC6A Expression and Facilitates Theca Cell Proliferation in Chicken Ovarian Follicles. DNA Cell Biol. 2017;36(11):922–9.

    Article  CAS  Google Scholar 

  9. Shen M, Sun H, Qu L, Ma M, Dou T, Lu J, Guo J, Hu Y, Wang X, Li Y, et al. Genetic architecture and candidate genes identified for follicle number in chicken. Sci Rep. 2017;7(1):16412.

    Article  Google Scholar 

  10. Hernandez A, Bahr J. Role of FSH and epidermal growth factor (EGF) in the initiation of steroidogenesis in granulosa cells associated with follicular selection in chicken ovaries. Reproduction. 2003;125(5):683–91.

    Article  CAS  Google Scholar 

  11. Wang J, Zhao C, Li J, Feng Y, Gong Y. Transcriptome analysis of the potential roles of FOXL2 in chicken pre-hierarchical and pre-ovulatory granulosa cells. Comp Biochem Physiol Part D Genomics Proteomics. 2017;21:56–66.

    Article  CAS  Google Scholar 

  12. Salzman J, Gawad C, Wang PL, Lacayo N, Brown PO. Circular RNAs are the predominant transcript isoform from hundreds of human genes in diverse cell types. PLoS One. 2012;7(2):e30733.

    Article  CAS  Google Scholar 

  13. Lu T, Cui L, Zhou Y, Zhu C, Fan D, Gong H, Zhao Q, Zhou C, Zhao Y, Lu D, et al. Transcriptome-wide investigation of circular RNAs in rice. RNA. 2015;21(12):2076–87.

    Article  CAS  Google Scholar 

  14. Ye CY, Chen L, Liu C, Zhu QH, Fan L. Widespread noncoding circular RNAs in plants. New Phytol. 2015;208(1):88–95.

    Article  CAS  Google Scholar 

  15. Hansen TB, Kjems J, Damgaard CK. Circular RNA and miR-7 in cancer. Cancer Res. 2013;73(18):5609–12.

    Article  CAS  Google Scholar 

  16. Hansen TB, Jensen TI, Clausen BH, Bramsen JB, Finsen B, Damgaard CK, Kjems J. Natural RNA circles function as efficient microRNA sponges. Nature. 2013;495(7441):384–8.

    Article  CAS  Google Scholar 

  17. Xu T, Wu J, Han P, Zhao Z, Song X. Circular RNA expression profiles and features in human tissues: a study using RNA-seq data. BMC Genomics. 2017;18(Suppl 6):680.

    Article  Google Scholar 

  18. Jeck WR, Sorrentino JA, Wang K, Slevin MK, Burd CE, Liu J, Marzluff WF, Sharpless NE. Circular RNAs are abundant, conserved, and associated with ALU repeats. RNA. 2013;19(2):141–57.

    Article  CAS  Google Scholar 

  19. Tao H, Xiong Q, Zhang F, Zhang N, Liu Y, Suo X, Li X, Yang Q, Chen M. Circular RNA profiling reveals chi_circ_0008219 function as microRNA sponges in pre-ovulatory ovarian follicles of goats (Capra hircus). Genomics. 2017;S0888-7543(17):30129–5.

    Google Scholar 

  20. Cheng J, Huang J, Yuan S, Zhou S, Yan W, Shen W, Chen Y, Xia X, Luo A, Zhu D, et al. Circular RNA expression profiling of human granulosa cells during maternal aging reveals novel transcripts associated with assisted reproductive technology outcomes. PLoS One. 2017;12(6):e0177888.

    Article  Google Scholar 

  21. Gilbert AB, Evans AJ, Perry MM, Davidson MH. A method for separating the granulosa cells, the basal lamina and the theca of the preovulatory ovarian follicle of the domestic fowl (Gallus domesticus). J Reprod Fertil. 1977;50(1):179–81.

    Article  CAS  Google Scholar 

  22. Langmead B, Salzberg SL. Fast gapped-read alignment with bowtie 2. Nat Methods. 2012;9(4):357–9.

    Article  CAS  Google Scholar 

  23. Memczak S, Jens M, Elefsinioti A, Torti F, Krueger J, Rybak A, Maier L, Mackowiak SD, Gregersen LH, Munschauer M, et al. Circular RNAs are a large class of animal RNAs with regulatory potency. Nature. 2013;495(7441):333–8.

    Article  CAS  Google Scholar 

  24. Gao Y, Zhang J, Zhao F. Circular RNA identification based on multiple seed matching. Brief Bioinform. 2017;19(5):803–10.

    Article  Google Scholar 

  25. Krzywinski M, Schein J, Birol I, Connors J, Gascoyne R, Horsman D, Jones SJ, Marra MA. Circos: an information aesthetic for comparative genomics. Genome Res. 2009;19(9):1639–45.

    Article  CAS  Google Scholar 

  26. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics. 2012;16(5):284–7.

    Article  CAS  Google Scholar 

  27. Agarwal V, Bell GW, Nam JW, Bartel DP. Predicting effective microRNA target sites in mammalian mRNAs. Elife. 2015;4:e05005.

    Article  Google Scholar 

  28. Wong N, Wang X. miRDB: an online resource for microRNA target prediction and functional annotations. Nucleic Acids Res. 2015;43(Database issue):D146–52.

    Article  CAS  Google Scholar 

  29. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504.

    Article  CAS  Google Scholar 

  30. Etches RJ, Petitte JN. Reptilian and avian follicular hierarchies: models for the study of ovarian development. J Exp Zool Suppl. 1990;4:112–22.

    Article  CAS  Google Scholar 

  31. Cogburn LA, Porter TE, Duclos MJ, Simon J, Burgess SC, Zhu JJ, Cheng HH, Dodgson JB, Burnside J. Functional genomics of the chicken--a model organism. Poult Sci. 2007;86(10):2059–94.

    Article  CAS  Google Scholar 

  32. Rybak-Wolf A, Stottmeister C, Glazar P, Jens M, Pino N, Giusti S, Hanan M, Behm M, Bartok O, Ashwal-Fluss R, et al. Circular RNAs in the mammalian brain are highly abundant, conserved, and dynamically expressed. Mol Cell. 2015;58(5):870–85.

    Article  CAS  Google Scholar 

  33. Sequencing ICG. Sequence and comparative analysis of the chicken genome provide unique perspectives on vertebrate evolution. Nature. 2004;432(7018):695–716.

    Article  Google Scholar 

  34. Ouyang H, Chen X, Wang Z, Yu J, Jia X, Li Z, Luo W, Abdalla BA, Jebessa E, Nie Q, et al. Circular RNAs are abundant and dynamically expressed during embryonic muscle development in chickens. DNA Res. 2017;25(1):71–86.

    Article  Google Scholar 

  35. Zhang X, Yan Y, Lei X, Li A, Zhang H, Dai Z, Li X, Chen W, Lin W, Chen F, et al. Circular RNA alterations are involved in resistance to avian leukosis virus subgroup-J-induced tumor formation in chickens. Oncotarget. 2017;8(21):34961–70.

    PubMed  PubMed Central  Google Scholar 

  36. Dang Y, Yan L, Hu B, Fan X, Ren Y, Li R, Lian Y, Yan J, Li Q, Zhang Y, et al. Tracing the expression of circular RNAs in human pre-implantation embryos. Genome Biol. 2016;17(1):130.

    Article  Google Scholar 

  37. Li L, Guo J, Chen Y, Chang C, Xu C. Comprehensive CircRNA expression profile and selection of key CircRNAs during priming phase of rat liver regeneration. BMC Genomics. 2017;18(1):80.

    Article  Google Scholar 

  38. Ashwal-Fluss R, Meyer M, Pamudurti NR, Ivanov A, Bartok O, Hanan M, Evantal N, Memczak S, Rajewsky N, Kadener S. circRNA biogenesis competes with pre-mRNA splicing. Mol Cell. 2014;56(1):55–66.

    Article  CAS  Google Scholar 

  39. Guo JU, Agarwal V, Guo H, Bartel DP. Expanded identification and characterization of mammalian circular RNAs. Genome Biol. 2014;15(7):409.

    Article  Google Scholar 

  40. Goitre L, Trapani E, Trabalzini L, Retta SF. The Ras superfamily of small GTPases: the unlocked secrets. Methods Mol Biol. 2014;1120:1–18.

    Article  CAS  Google Scholar 

  41. Jiang C, Diao F, Sang YJ, Xu N, Zhu RL, Wang XX, Chen Z, Tao WW, Yao B, Sun HX, et al. GGPP-mediated protein Geranylgeranylation in oocyte is essential for the establishment of oocyte-granulosa cell communication and primary-secondary follicle transition in mouse ovary. PLoS Genet. 2017;13(1):e1006535.

    Article  Google Scholar 

  42. Xu R, Qin N, Xu X, Sun X, Chen X, Zhao J. Inhibitory effect of SLIT2 on granulosa cell proliferation mediated by the CDC42-PAKs-ERK1/2 MAPK pathway in the prehierarchical follicles of the chicken ovary. Sci Rep. 2018;8(1):9168.

    Article  Google Scholar 

  43. Fukata M, Kaibuchi K. Rho-family GTPases in cadherin-mediated cell-cell adhesion. Nat Rev Mol Cell Biol. 2001;2(12):887–97.

    Article  CAS  Google Scholar 

  44. Schuster MK, Schmierer B, Shkumatava A, Kuchler K. Activin a and follicle-stimulating hormone control tight junctions in avian granulosa cells by regulating occludin expression. Biol Reprod. 2004;70(5):1493–9.

    Article  CAS  Google Scholar 

  45. Peng Y, Chang L, Wang Y, Wang R, Hu L, Zhao Z, Geng L, Liu Z, Gong Y, Li J, et al. Genome-wide differential expression of long noncoding RNAs and mRNAs in ovarian follicles of two different chicken breeds. Genomics. 2018;S0888-7543(18):30268–4.

    Google Scholar 

  46. Kang L, Cui X, Zhang Y, Yang C, Jiang Y. Identification of miRNAs associated with sexual maturity in chicken ovary by Illumina small RNA deep sequencing. BMC Genomics. 2013;14(1):352.

    Article  CAS  Google Scholar 

  47. Bachmayr-Heyda A, Reiner AT, Auer K, Sukhbaatar N, Aust S, Bachleitner-Hofmann T, Mesteri I, Grunt TW, Zeillinger R, Pils D. Correlation of circular RNA abundance with proliferation--exemplified with colorectal and ovarian cancer, idiopathic lung fibrosis, and normal human tissues. Sci Rep. 2015;5:8057.

    Article  CAS  Google Scholar 

  48. Quan G, Li J. Circular RNAs: biogenesis, expression and their potential roles in reproduction. J Ovarian Res. 2018;11(1):9.

    Article  Google Scholar 

  49. Veno MT, Hansen TB, Veno ST, Clausen BH, Grebing M, Finsen B, Holm IE, Kjems J. Spatio-temporal regulation of circular RNA expression during porcine embryonic brain development. Genome Biol. 2015;16:245.

    Article  Google Scholar 

  50. Ceriani M, Scandiuzzi C, Amigoni L, Tisi R, Berruti G, Martegani E. Functional analysis of RalGPS2, a murine guanine nucleotide exchange factor for RalA GTPase. Exp Cell Res. 2007;313(11):2293–307.

    Article  CAS  Google Scholar 

  51. OS A, Parrini MC, Camonis J. RalGPS2 Is Essential for Survival and Cell Cycle Progression of Lung Cancer Cells Independently of Its Established Substrates Ral GTPases. PLoS One. 2016;11(5):e0154840.

    Article  Google Scholar 

  52. Liu L, Kong N, Xia G, Zhang M. Molecular control of oocyte meiotic arrest and resumption. Reprod Fertil Dev. 2013;25(3):463–71.

    Article  CAS  Google Scholar 

  53. Johnson AL, Bridgham JT. Regulation of steroidogenic acute regulatory protein and luteinizing hormone receptor messenger ribonucleic acid in hen granulosa cells. Endocrinology. 2001;142(7):3116–24.

    Article  CAS  Google Scholar 

  54. McDerment NA, Wilson PW, Waddington D, Dunn IC, Hocking PM. Identification of novel candidate genes for follicle selection in the broiler breeder ovary. BMC Genomics. 2012;13(1):494.

    Article  CAS  Google Scholar 

  55. Chiu HS, Martinez MR, Bansal M, Subramanian A, Golub TR, Yang X, Sumazin P, Califano A. High-throughput validation of ceRNA regulatory networks. BMC Genomics. 2017;18(1):418.

    Article  Google Scholar 

Download references

Acknowledgments

We would like to thank Huiqiang Shi and Yongjuan Wang from Jiangsu Jinghai Poultry Industry Co. Ltd. for assistance with the recording of animal data.

Funding

This work was supported by the National Broiler Industrial and Technology System (nycytx-42-G1–05), Priority Academic Program Development of Jiangsu Higher Education Institutions, Natural Science Foundation of Jiangsu Province (BK20181217), and Institutes of Agricultural Science and Technology Development, PR. China. The funders had no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Availability of data and materials

The datasets generated and/or analyzed during the current study are available in the NCBI Sequence Read Archive (SRA) repository with the accessions SRP153771, BioProject accession number PRJNA481176 (https://www.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA481176). Additional dataset(s) supporting the conclusions of this article are included within the article and its additional files.

Author information

Author notes

  1. Manman Shen and Tingting Li contributed equally to this work.

Authors and Affiliations

  1. College of Animal Science and Technology, Yangzhou University, Yangzhou, 225009, China

    Manman Shen, Tingting Li, Genxi Zhang, Pengfei Wu, Fuxiang Chen, Qiuhong Lou, Lan Chen, Xuemei Yin, Tao Zhang & Jinyu Wang

  2. Jiangsu Institute of Poultry Science, Chinese Academy of Agricultural Science, Yangzhou, 225216, China

    Manman Shen

Authors
  1. Manman Shen
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Tingting Li
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Genxi Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Pengfei Wu
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Fuxiang Chen
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Qiuhong Lou
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Lan Chen
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Xuemei Yin
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Tao Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. Jinyu Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

JW, GZ and MS conceived and designed the work. MS and TL performed the transcriptome data, did bioinformatics analyses and wrote the paper. PW, FC and LC prepared the Figs. QL, XY and TZ collected samples and performed qRT-PCR experiment. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Genxi Zhang or Jinyu Wang.

Ethics declarations

Ethics approval and consent to participate

The experimental chickens were from Jiangsu Jinghai Yellow Chick Co. Ltd., located in Nangtong, Jiangsu Province. All experimental procedures were conducted in compliance with Experimental Animals Regulations and all efforts were made to alleviate the suffering of the birds. Ethics approval for this study was granted by the Animal Care Committee of Yangzhou University (Yangzhou, China) with permit number SYXK (Su) 2012–0029 from the Chinese Ministry of Science and Technology.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1:

Table S1. Detailed information of primers used in this study. Table S2. All novel circRNAs detected in granulosa cells. Table S3. Flank intron length statistics. Table S4. Differentially expressed circRNAs (Fold change > 2, q-value < 0.05) in different granulosa cells. Table S5. Gene ontology for host genes of differential expressed circRNAs. Table S6. KEGG pathway for host genes of all differentially expressed circRNAs in three granulosa cells stage. Table S7. Overlapped prediction genes by TargetScan and miRDB. Table S8. Gene ontology of predicted target genes based on circRALGPS2. Table S9. KEGG pathway of predicted target genes based on circRALGPS2. (XLSX 1280 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Shen, M., Li, T., Zhang, G. et al. Dynamic expression and functional analysis of circRNA in granulosa cells during follicular development in chicken. BMC Genomics 20, 96 (2019). https://doi.org/10.1186/s12864-019-5462-2

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-019-5462-2

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us