Genome-Wide Identification of CircRNAs of Infective Larvae and Adult Worms of Parasitic Nematode, Haemonchus contortus

Introduction

Circular RNAs (circRNAs), a novel class of endogenous non-coding RNAs with a covalently closed loop structure without 3′ end poly-adenylated tail and the 5′ end cap structure, have been widely discovered in the eukaryotic transcriptome. Although circRNAs with circular structure in eukaryotic cells were observed by electron microscopy more than four decades ago (Hsu and Coca-Prados, 1979), and have been reported sporadically thereafter (Grabowski et al., 1981; Capel et al., 1993), circRNAs were deemed to be functionless side-products of transcription or artifacts of aberrant RNA splicing. Owing to dramatic advances of high-throughput sequencing technology and specialized bioinformatic analysis pipelines, many circRNAs have been frequently recognized in diverse species, such as mammals (Guo et al., 2014), Drosophila (Westholm et al., 2014), Caenorhabditis elegans (Memczak et al., 2013), plants (Zhang et al., 2020) and Archaea (Danan et al., 2012). Characteristics analysis revealed that circRNAs could be generated from exonic, intronic, intergenic regions, and exonic-intronic regions through non-canonical splicing event called back-splicing (Li et al., 2018a; Kristensen et al., 2019), and the expression patterns were often cell, tissue, and developmental stage-specific (Salzman et al., 2013; Liang et al., 2017; Nicolet et al., 2018). The research on circRNAs has become one of the amazing hotspots in life science and the knowledge about circRNAs has been renovated uninterruptedly and fleetly.

More recently, the biological roles of circRNAs were continuously illuminated. Firstly, one of the prominent roles is that abundant circRNAs located in the cytoplasm and harbored miRNA binding sites can serve as miRNA sponges or decoys to regulate gene expression. For instance, in the mouse brain, ciRS-7/CDR1as was demonstrated to be a miR-7 sponge and to modulate miR-7’s biological activities and functions (Hansen et al., 2013). Secondly, some circRNAs can also act as protein sponges or decoys to indirectly regulate protein’s functions. For example, in mouse embryo fibroblast and cardiac fibroblast, circ-Foxo3 retarded cell cycle progression through interacting with p21 and CDK2 (Du et al., 2016). Thirdly, a novel subclass of circRNAs, ElciRNAs, have been demonstrated to promote the transcription of their parent genes by interacting with U1 snRNP and RNA Polymerase II in the promoter region (Li et al., 2015). Except for regulating gene transcription and protein expression, some circRNAs can be directly translated into detectable proteins. For example, in murine and human myoblast, circ-ZNF609 could be translated into a protein to control myoblast proliferation (Legnini et al., 2017). However, for most circRNAs, their functions remain unknown.

In contrast to vertebrate and plants, to date, the research on circRNAs in nematode is still limited. Only in the free-living nematode C. elegans, circRNAs have been identified from different developmental stages including sperm, oocytes, one/two-cell embryos (Memczak et al., 2013), the fourth-stage larvae (L4) (Cortes-Lopez et al., 2018), and adults (Ivanov et al., 2015). However, up to date, there is no information available for circRNAs of parasitic nematodes including the blood-sucking nematode Haemonchus contortus, which infects small ruminants and causes significant economic losses to the livestock industry around the world. This parasite has a direct and rapid life cycle (Laing et al., 2013; Schwarz et al., 2013). The adult worms live in the mucosa of host abomasum and females produce eggs that are excreted from host feces. The eggs hatch and develop to the infective third-stage larvae (L3) through first- and second-stage larvae in around seven days. The L3s are transmitted orally from contaminated herbage to the grazing host, exsheath, and develop to dioecious adults through L4 in two to three weeks. L4s and adult worms feed on blood from capillaries in the ruminant abomasum mucosa.

In order to study whether circRNAs play a key role in the development of gastrointestinal parasitic nematodes, it is necessary to analyze the expression and differences of circRNAs in different developmental stages/sexes. In the current study, circRNAs were sequenced from L3s, female and male worms of H. contortus. Their genomic feature, length distribution, and the number of exons coding for circRNAs were analyzed. Bioinformatic functional analyses of circRNAs were also conducted. In addition, Sanger sequencing was employed to validate the head-to-tail back-splicing site, and qRT-PCR was used to verify the differentially expressed circRNAs and the resistance of circRNAs to RNase R digestion. This is the first identification of circRNAs in parasitic nematodes, the data provide a basis for further studying the functions of circRNAs in H. contortus and other related nematode species.

Materials and Methods

Parasite Sample Collection

Haemonchus contortus L3s and adult worms were collected as described previously (He et al., 2019). Briefly, experimental goats (3 months of age) purchased from Hubei Academy of Agricultural Sciences were confirmed helminth-free by fecal examination and maintained under a parasite-free condition. The goats were orally inoculated with 8000 infective L3 of H. contortus (Haecon-5). Four weeks post-infection, fecal samples were harvested and incubated at 27°C for 1 week, then iL3s were collected by the Baermann funnel method. Adult female and male worms were collected from the abomasum after euthanasia and necropsy of infected goats according to their distinctive appearances. All samples were washed thoroughly in sterile phosphate-buffered saline and then transferred to liquid nitrogen for storage until use.

Extracting RNA, Strand-Specific Library Preparation, and Sequencing

Total RNA samples were extracted, respectively, from three independent experimental replicates of L3s (40,000 larvae per replicate), male adults (Am), and female adults (Af) (30 worms per replicate) of H. contortus using TranZol (Sigmen, China) with the manufacturer’s instructions. The total RNA yields and quality were evaluated on the Agilent 2100 Bioanalyzer (Agilent, USA). Then, approximately 5 μg of total RNA from each sample was treated to remove ribosomal RNA (rRNA) using Ribo-zero™ rRNA Removal Kit (Illumina, USA) following the manufacturer’s procedure, and further incubated with RNase R (Epicentre, USA) to digest linear RNAs. The remaining RNAs were used as templates for generating circRNA-seq libraries. Finally, all eligible libraries were sequenced using the Illumina HiSeq X Ten platform to generate 100 bp pair-end reads at BGI Technology Co. Ltd (Shenzhen, China). The samples were named as Af_1, Af_2, Af_3; Am_1, Am_2, Am_3; L3_1, L3_2, L3_3, respectively.

Identification of Differently Expressed circRNAs

Firstly, all the raw data were removed with low-quality reads, unknown bases and adaptors, and then FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) was used to confirm the quality of reads. The Q20, Q30, and GC contents were also simultaneously calculated. Subsequently, Bowtie2 (Langmead et al., 2009) and BWA (Li and Durbin, 2009) were employed to map the clean reads to the publicly available H. contortus reference genomes, which was downloaded from WormBase Parasite website (https://parasite.wormbase.org/Haemonchus_contortus_prjeb506/). Lastly, find_circ (Memczak et al., 2013) and CIRI2 (Gao et al., 2018) tools with default parameters were performed to identify circRNAs. Only the candidate circRNAs were identified by both algorithm methods in all replicates for further calculation of expression level.

The expression level of circRNAs from each sample was quantified based on the read counts, normalized using transcripts per million (TPM) as previously described (Zhou et al., 2010) and calculated as the normalized expression level = (read count*1,000,000)/the sum of circRNAs read count. The differentially expressed circRNAs were identified by pairwise comparison between L3s, adult female and male worms using DESeq2 algorithm (Love et al., 2014), and the cut-off criterion was set as |log₂^FoldChange| ≥ 1 and padj < 0.05. Venn, Heat map clustering, circos, and volcano diagrams were drawn using correspondence R packages.

Bioinformatics Functional Analysis

Based on the Gene Ontology (GO) (Young et al., 2010) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al., 2017) databases, the potential functional analysis of all the parental genes of differentially expressed circRNA was carried out. GO annotation was performed by BLAST (E-value < 10^-5) using C. elegans orthologs, which was then assigned to terms using clusterProfiler R packages (Yu et al., 2012). The padj<0.05 denotes the threshold representing the significantly enriched GO terms. All enriched GO terms that possessed a padj < 0.05 were displayed. The KOBAS (http://kobas.cbi.pku.edu.cn/) was used to test the statistical enrichment of the circRNA-host genes in the KEGG pathways. The top 20 enriched pathways were selected based on the p-value and rich factor.

The Construction of circRNA-miRNA Network

To determine whether circRNAs in H. contortus could function as miRNA sponges to regulate their target mRNAs, 194 mature miRNA sequences of H. contortus were downloaded from miRbase 22 (https://mirbase.org/ftp.shtml). The miRNA targeting sites for circRNAs were identified by miRanda algorithm (Enright et al., 2003) with threshold parameters as follows: max free energy values <− 20 and a score of 150 or higher. Based on the predicted results, network of circRNA–miRNA was then constructed and visualized using Cytoscape 3.7.1 software (Shannon et al., 2003).

Validation of Back-Splicing Site of circRNAs

The total genomic DNA (gDNA) was extracted from H. contortus using the Blood/Cell/Tissue Genomic DNA Extraction Kit (Tiangen, China) with the manufacturer’s instructions. The total RNA was extracted from L3s, adult female or male worms of H. contortus, respectively and then cDNA was reverse transcribed using HiScript II Q RT SuperMix (Vazyme, China) following the manufacturer’s protocol.

CircRNAs were validated by PCR as previously described (Memczak et al., 2013). To do it, divergent and convergent primers were designed, respectively. Linear transcripts were detected by convergent primers, and the candidate circular templates were detected using divergent primers. If the band with expected size was amplified with convergent primers from both cDNA and gDNA templates, but only amplified with divergent primers from cDNA template, it would suggest the presence of the back-splicing junction. For each PCR amplification, cDNA or gDNA was used with 2×Taq Master Mix (Vazyme, China), and 40 cycles of PCR cycling condition were performed. PCR products were examined by 1% agarose gel electrophoresis, and those with expected bands and amplified with divergent primers only from cDNA templates were sent to TSING Ke Biotech Co., Ltd. (Wuhan, China) for Sanger sequencing. Both specific divergent and convergent primers were synthesized by TSING Ke Biotech Co., Ltd. (Wuhan, China), and the primer information was listed in Additional File 1: Table S1.

Quantification of circRNAs With Quantitative Real-Time PCR

To test accuracy of high-throughput sequencing, differentially expressed circRNAs with the junction site verified were selected in each comparison group to check by qRT-PCR assay. In addition, for confirmation of the resistance of circRNAs to the exoribonuclease RNase R, the qRT-PCR experiment was also performed using RNA samples with or without RNase R treatment. The total RNA was incubated with 4U/μg of RNase R (RNR07250, Epicentre) for 15 min at 37°C, 10 min at 70°C, and mock treatment was carried out in the same conditions without RNase R treatment. The qRT-PCR experiment was performed using 2× SYBR Green Master Mix (Takara, China), and the expression levels of circRNAs was determined with divergent primers using the ViiA™ 7 Real-Time PCR System (Applied Biosystems, USA). 2^-△△Ct method was used to calculate the relative expression level of each target circRNA and the β-tubulin gene was an internal control. The qRT-PCR experiment was carried out with the following reaction system: 10 μL of the 2× SYBR Green, 0.4 μL dye, 8 μL of the cDNA template and RNase-free ddH₂O, 0.8 μL of the forward and reverse primers, respectively. The divergent primers were synthesized by TSING Ke Biotech Co., Ltd. (Wuhan, China), and the primer information was listed in Additional File 1: Table S2.

Statistical Analysis

These works were created by GraphPad Prism 9.0 (GraphPad Company, USA). Student’s t-test was performed to compare the expression levels of target circRNAs in L3s, adult female and male worms and the analysis of the results was with the mean ± SEM. The results were statistically significant at *P<0.05, **P<0.01, ***P < 0.001, ****P < 0.0001; ns, no significant.

Results

General Properties of H. contortus circRNAs

In order to identify H. contortus circRNAs, nine libraries were constructed and sequenced from L3s, adult female (Af) and male worms (Am) of this parasite. Clean data were used in subsequent analysis after wiping off the adaptor and low-quality reads. Among the nine libraries, the proportion of clean reads mapped to the H. contortus genome ranged from 61.07% to 74.73% (Additional File 1: Table S3). A total of 20073 candidate circRNAs with more than two back-spliced junction reads were identified and annotated from L3s, Af and Am worms using find_circ and CIRI2 computational pipelines (Additional File 1: Table S4). Among these circRNAs, 3491, 2550, and 18883 were detected in Af, Am, and L3 stages, respectively, and only 1650 were prevailed in all three stages/sexes, suggesting that most circRNAs in H. contortus were stage-specific (Figure 1A).

FIGURE 1

Figure 1 Characteristics of the identified circRNAs in H. contortus. (A) Venn diagram presenting the amount of circRNAs in third-stage larvae (L3), adult female worms (Af), and adult male worms (Am). (B) The numbers of circRNAs derived from exonic, intronic, and intergenic regions in L3, Af, and Am, respectively. (C) Circos plot showing the distribution of circRNAs on six H. contoruts chromosomes. Light pink bar represented circRNAs from plus strand; pale green bar from negative strand. Each track from the periphery to the core represents: chromosomes; L3; Af; Am. (D) Distribution of amount of circRNAs derived from different number of exons. (E) Length range distribution of the identified circRNAs. (F) The number of circRNAs per gene.

The subsequent analysis on genomic origin of circRNAs revealed that 14314 circRNAs (71.31%) were generated from exonic regions and 4476 circRNAs (22.30%) were generated from intergenic regions, but only 1283 circRNAs (6.39%) were stemmed from intronic regions (Figure 1B). In addition, the distribution of genome region coding for circRNAs on each chromosome was further examined using circlize R package. The mapping results showed that genome regions coding for these circRNAs were widely and unevenly distributed across the whole H. contortus chromosomes and were more abundant on the chromosome X than on other chromosomes. There is no correlation between the distribution of circRNAs and chromosome length (Figure 1C).

Next, the amount of exon-derived circRNAs were counted and showed that most of them typically encompass less than five exons. More than 94.18% of circRNAs were formed by multiple exons (up to 10 exons) (Figure 1D). Considering the length distribution, most circRNAs (83.48%) ranged from 200 bp to 600 bp (Figure 1E) and the average length of circRNAs was 354 bp, while the maximum length was 1875 bp. Moreover, the analysis on alternative circularization suggested that most genes yielded one or two circRNAs, some genes generated multiple distinguishing circularized products. 4937 genes encoded at least one circRNA, among them, 41.74% of genes only produced one circRNA, and 94.75% produced no more than 10 circRNAs (Figure 1F).

Expression Analysis of circRNAs

To analyze the expression pattern of circRNAs during H. contortus growth process, the expression level was quantified and differentially expressed circRNAs were screened by pairwise comparison between L3, adult female and male worms using DEGseq2. The hierarchical clustering presented the expression changes of circRNAs during three stages/sexes of H. contoruts. Overall, the abundance of most circRNAs decreased during development from L3 to adult stage (Figure 2A). 15948 circRNAs were identified differentially expressed between Af_vs_L3 group, shown in the volcano plot. Among them, 1707 were up-regulated and 14241 were down-regulated (Figure 2B). In the comparison between Am_vs_L3, 16847 circRNAs were identified as differentially expressed, including 1356 up-regulated and 15491 down-regulated (Figure 2C). Among the differentially expressed circRNAs, 13409 commonly existed in both Af_vs_L3 and Am_vs_L3 comparisons. In addition, comparison between Af_vs_Am revealed 1119 differentially expressed circRNAs. Among them, 637 were up-regulated and 482 were down-regulated (Figure 2D).

FIGURE 2

Figure 2 Differential expression pattern of circRNAs in H. contortus. (A) Hierarchical clustering analysis of differentially expressed circRNAs across L3, Af, and Am samples. The expression values were represented in different colors. Red strip represented high relative expression and green strip indicated low relative expression. (B–D) The volcano plot for differentially expressed circRNAs in Af_vs_L3, Am_vs_L3, Af_vs_Am comparisons. The volcano plot was constructed using log2 (Fold Change) and padj. light blue points: up-regulated circRNAs; red points: down-regulated circRNAs; purple points: no difference. The cut-off criterion was set as Fold Change ≥ 2 and padj< 0.05.

Functional Annotation Analysis of Differentially Expressed circRNAs

Analysis of the function of parental genes could provide valuable clues about circRNA’s functions, thus, GO and KEGG annotation analyses were conducted. GO categories with padj less than 0.05 were assigned to the parental genes, and the significantly enriched genes of each GO term were statistically analyzed. In Af_vs_L3 group (Figure 3A), source genes of differentially expressed circRNAs were divided into thirty-two significant GO terms, nineteen of them were involved in the biological process category. The remarkably enriched GO terms include “signaling” (GO:0023052), “single organism signaling” (GO:0044700), and “signal transduction” (GO:0007165), “response to stimulus” (GO:0050896). Only three GO terms, “cell part” (GO:0044464), “cell” (GO:0005623), and “membrane” (GO:0016020), were annotated in cellular component subcategories. As for the molecular function, ten GO subcategories were enriched and mainly included “protein binding” (GO:0005515), “receptor activity” (GO:0004872), and “signal transducer activity” (GO:0004871). In Am_vs_L3 comparison group (Figure 3B), 14, two, and 11 terms were categorized into biological process, cellular component, and molecular function, respectively. The significant enrichment GO terms were like those of the Af_vs_L3 group and no significant enrichment GO terms were identified in the Af_vs_Am group.

FIGURE 3

Figure 3 GO and KEGG enrichment analysis for differentially expressed circRNAs of H. contortus. (A, B) The most enriched GO terms in Af_vs_L3 and Am_vs_L3. Y-axis represented GO terms; X-axis represented the value of -log10 (padj). Different colors represented biological_process, cellular_component, and molecular_function, respectively. (C) The top 20 enriched KEGG pathways. Y-axis represents pathways; X-axis represents rich factor; The color and size of each bubble represent pvalue and the number of genes enriched in a pathway, respectively. Af_vs_L3, Am_vs_L3, and Af_vs_Am are represented from left to right.

Accordingly, the KEGG pathway analysis was further carried out to yield significantly enriched pathways. The top 20 enriched pathways are displayed in an enriched scatter diagram (Figure 3C). In Af_vs_L3 and Am_vs_L3 groups, the significantly represented pathways were pyruvate metabolism, MAPK signaling pathway, glycerophospholipid metabolism, and phosphatidylinositol signaling system, implying that source genes in H. contortus were mainly involved in signal processing and associated with the longevity of worm. However, the parent genes in Af_vs_Am comparison group were primarily related to metabolisms, such as pyruvate metabolism, propanoate metabolism, purine metabolism, fatty acid biosynthesis, and steroid biosynthesis.

Construction of Complete Endogenous (ceRNA) Network

Mounting researches reported that circRNAs might act as miRNAs sponges and likely modulate their activities (Yang et al., 2019; Tang et al., 2021; Xin et al., 2021). Meanwhile, miRNAs could be involved in multiple physiological processes during nematode development (Zhang et al., 2018; Mani et al., 2021). Information on the potential interaction between circRNA and miRNA will be helpful for understanding the function and regulatory mechanism of circRNA in H. contortus. However, the expression profiles of miRNAs were not available in the databases for different developmental stages of H. contortus (only 194 miRNAs available from L3 and mixed-sex adult worms), the potential miRNA targets can only be predicted based on sequence complementarity via miRanda algorithm. The predicted network consisted of 12015 interacting pairs and 8268 nodes, and the nodes included 8077 differentially expressed circRNAs and 191 miRNAs except for hco-miR-5892a, hco-miR-5892b, and hco-miR-5898-5p (Figure 4A and Additional File 1: Table S5). This indicated that a single miRNA could be targeted by different circRNAs, and that single circRNA could also sponge more than one miRNA. For example, hco-lin-4, hco-miR-1, and hco-miR-124 could be targeted by 36, 19, and 102 circRNAs, respectively. 55.89% of circRNAs contained more than one different miRNA-binding site, and hco_circ_00016196 harbored the greatest number (nine) of different miRNA binding sites (Figure 4B). The vast majority of circRNAs had one or two binding sites for one particular miRNA, only hco_circ_00015227 had 10 binding sites for one miRNA (hco-miR-307) (Additional File 1:Table S6).

FIGURE 4

Figure 4 (A) The interaction network of differentially expressed circRNAs and miRNAs for H. contortus. Red squares represent miRNAs, and blue nodes represent circRNAs. (B) The circRNA, hco_circ_00016196, and nine miRNAs interaction network. Red nodes represent circRNAs, and light blue squares represent miRNAs.

Validation of circRNAs by Sanger Sequencing

For validating individual circRNA’s back-splicing junction site predicted from the analysis pipeline, divergent and convergent primers were designed for candidate circRNAs and used in PCR amplification with cDNA and gDNA as templates, respectively. From cDNA template, PCR products with expected sizes were amplified with both convergent primers (black opposing triangle pairs) and divergent primers (black back-to-back triangle pairs). However, from gDNA, PCR products with expected sizes were only exclusively amplified with convergent primers. The amplified products were then sequenced and the correctness of head-to-tail back splicing site was confirmed. In total, twenty-four circRNAs were verified. The original sequences of parental genes of 10 selected circRNAs and the sequencing chromatogram of PCR products amplified from cDNA or gDNA were shown (Figure 5 and Additional File 1: Figure S1).

FIGURE 5

Figure 5 Experimental validation of back-spliced junction sites of candidate circRNAs. Linear RNAs were successfully amplified by convergent primers from both cDNA and gDNA. CircRNAs were successfully amplified by divergent primers from cDNA but failed to be amplified from gDNA. PCR products were detected by 1% agarose electrophoresis. The head-to-tail back-splicing site was further confirmed by Sanger sequencing. Black inverted triangles represented the back-spliced junction loci, and β-tubulin gene was used for the control group. Marker: 2000 bp ladder DNA marker.

Validation of circRNA Quantification and RNase R Digestion

To verify the reliability of high throughput sequencing data, differentially expressed circRNAs with the back-splicing junction sites validated were selected in each comparison group to assess the differential expression level by qRT-PCR (Figure 6A). In Af_vs_L3 group (three up-regulated and two down-regulated in adult female worms) and Am_vs_L3 group (two up-regulated and three down-regulated in adult male worms), the expression patterns of the selected circRNA detected by qRT-PCR were consistent with those from high-throughput sequencing results. In Af_vs_Am group (three up-regulated and two down-regulated circRNAs in adult female worms), qRT-PCR results also indicated a similar expression tendency with high-throughput sequencing results except for hco_circ_00058255 for which no difference was detected by qRT-PCR. In general, the correlation between qRT-PCR results and RNA-Seq is 0.85, suggesting that the high throughput sequencing results were dependable.

FIGURE 6

Figure 6 Reliability verification of RNA-Seq and the resistance of circRNAs to RNase R digestion. (A) The validation of differentially expressed circRNAs in all three comparisons. All experiments were conducted in three independent biological replicates and presented as means ± SEM. The expression levels among different groups were compared by Student’s t-test. *P < 0.05; **P < 0.01; ***P < 0.001, ****P < 0.0001. ns, no significant. (B) The results indicated that circRNAs were detected in both RNase R treated and control RNA samples, whereas the linear control β-tubulin gene wasn’t detected in RNase R treated sample. All experiments were conducted in two independent biological replicates and presented as means ± SEM.

RNase R, a 3’ to 5’ exoribonuclease from the E. coli RNR superfamily, can digest almost all linear RNAs but not circRNAs. To assess the resistance of H. contortus circRNAs to RNase R digestion, total RNA was treated with RNase R and then used for cDNA synthesis followed by qRT-PCR detection. The transcription level of linear RNA control (β-tubulin) was significantly reduced after RNase R treatment, nevertheless, the transcription of circRNAs was only slightly decreased or increased (Figure 6B). This result confirmed that H. contortus circRNAs were endowed with a strong resistance to RNase R treatment.

Alternative Circularization of circRNAs

Different circRNAs can be generated from the same source gene consisting of different exons and/or introns. To do alternative splicing analysis, two-parent genes (HCON_00020120 and HCON_00002920) each generating two circRNA forms were selected. HCON_00020120 generated both hco_circ_00007037 and hco_circ_00007039, which had different acceptors but shared the same donor. In contrast, HCON_00002920 produced both co_circ_00001063 and hco_circ_00001064, which had different donors but shared the same acceptor. Divergent primers for these four circRNAs were designed as illustrated in Figure 7A. The PCR products were amplified using divergent primers from cDNA and then sequenced to validate the predicted back-splicing site. Primer-1 and primer-4 only amplified hco_circ_00007037 and hco_circ_00001064, respectively, but primer-2 and primer-3 each amplified two circRNAs (hco_circ_00007037 and hco_circ_00007039 by primer-2) and (hco_circ_00001063 and hco_circ_00001064 by primer-3), respectively (Figure 7B). The sizes of PCR products and sequences in the back-splicing site of each PCR product were confirmed by agarose gel electrophoresis and Sanger sequencing, respectively (Figures 7B–D). Primer sequence and length of PCR products were listed in Additional File 1: Table S7.

FIGURE 7

Figure 7 Visualization and validation of two alternative back-splicing circularization events. (A) A schematic diagram of the divergent primers designed for alternative splicing event. (B) Representative PCR products amplified with primers presented in (A) were examined by 3% agarose electrophoresis. M: 500 bp ladder DNA marker. (C, D) Sanger sequencing to validate the backsplice sites. Vertical Bar represents the backsplice junction site.

Discussion

With the advances of sequencing technology, abundant resources of genomics and transcriptome have been accumulated for nematodes, while the research on the regulation of gene expression mediated by non-coding RNA in parasitic nematodes is still very scarce. High-throughput sequencing technology and bioinformatics methods make it possible to identify novel types of transcripts, even those often with low expression levels. Previous studies have identified the mRNA and miRNA expression profiles for key growth stages/sexes of H. contortus (Winter et al., 2012; Laing et al., 2013; Ma et al., 2018; Marks et al., 2019). Nevertheless, circRNAs, a recently discovered novel class of non-coding RNA from various animals and plants, have yet been studied in any parasitic nematodes including H. contortus. How circRNAs are transcribed in H. contortus and whether their expressions participate in the regulation of developmental processes of this parasite remain totally unknown. The present study represents the first one to perform genome-wide identification and potential function analysis of circRNAs in parasitic nematodes.

In the present study, from L3s, adult female and male worms of H. contortus, a total of 20073 circRNA were identified, the number of which is notably much higher than those from C. elegans (1166) (Cortes-Lopez et al., 2018) and Plasmodium falciparum (1381) (Broadbent et al., 2015). One of the reasons causing this was likely that our study performed an experimental strategy of ribo-depleted RNA-Seq coupled with RNase R enrichment. A previous study showed that from RNase R treated samples, indeed more candidate circRNAs were identified with approximately equal sequencing depth and prediction tools (Zeng et al., 2017). In our study, 2550, 3491, and 18883 circRNAs were detected in adult male, adult female, and L3 stages, respectively, and only 1650 circRNAs existed in all developmental stages/sexes, indicating that most circRNAs are stage-specific. However, what remains unclear is why the amount of circRNAs in L3 was significantly higher than those in the adult worms. We have checked out the analyses many times and ruled out the possibility of contamination. The reasons causing the significant difference in the number of circRNAs between L3 and adult worms hopefully can be discovered in the future when more circRNAs are identified from a wide range of parasitic nematode species.

Notoriously, based on the derivation of circRNAs in the genome, they can be divided into exonic, intronic, and intergenic ones (Kristensen et al., 2019). Among the identified 20073 circRNAs, the exon-derived were predominant (71.31%), consistent with the results from other species, such as Arabidopsis thaliana (Chen et al., 2017), C. elegans (Cortes-Lopez et al., 2018), zebrafish (Liu et al., 2019), and human (Xu et al., 2017). However, in wheat (Wang et al., 2016) and kiwifruit (Wang et al., 2016; Wang et al., 2017), intergenic circRNAs were in the majority, suggesting that the way that these genomes encode circRNAs occurs in a species-specific manner and the molecular basis of circRNA biogenesis is highly different among different species. In addition, a single gene locus could selectively generate multiple circRNAs through alternative back-splicing, and the events of alternative circularization phenomenon have been experimentally validated in human (Feng et al., 2019), fruit fly (Gao et al., 2016), and rice (Lu et al., 2015). Here, the alternative splicing events were analyzed and experimentally verified, indicating that such splicing mechanisms were conserved in parasitic nematodes.

CircRNAs were usually produced from protein-coding genes so that they may function in the same pathway as the source genes (Wilusz, 2017). The comparative analyses between adult worms and L3 revealed that significantly enriched GO terms for source genes were related to signaling, signal transduction, and response to stimulus, indicating that circRNAs could play roles in parasite’s response to environmental signals (e.g., temperature, humidity) during its transition from the free-living stage to the parasitic stage in order to adapt to the host abomasum environment. The KEGG pathway analysis identified several important pathways which are mainly involved in H. contortus growth, including MAPK signaling pathway, phosphatidylinositol signaling system, FoxO signaling pathway. In addition, interestingly, an identified exonic circRNA (hco_circ_00047001) was generated from the serine/threonine-specific protein kinase encoding gene, akt, which was known to associate with larval development (Di et al., 2020), suggesting that this circRNA may also play a role in larval development. Furthermore, hco_circ_00006206 was transcribed from HCON_00017820 (Hc-aap-1, Li et al., 2014), and four circRNAs (hco_circ_00012209, hco_circ_00012210, hco_circ_00012216, and hco_circ_00012219) derived from HCON_00034230 (Hc-age-1, Li et al., 2014). Hc-aap-1 and Hc-age-1 are important genes in the insulin-like signaling pathway of H. contortus, which play significant roles in the development of this parasite, especially during its transitional change from the environment to the host abomasum (Li et al., 2014; Li et al., 2016; Ma et al., 2018; Di et al., 2020). These results suggested that circRNAs could have important biological functions in the development and environmental information processing of H. contortus.

Numerous researches have indicated that circRNAs containing miRNA binding sites act as miRNA decoys or sponges to competitively suppress miRNA activity on their targets (Hansen et al., 2013; Memczak et al., 2013). In the present study, a big number of differentially expressed circRNAs contained one or more potential miRNA binding sites, suggesting that these circRNAs may function as miRNA sponges in H. contortus. In C. elegans, cel-let-7 was reported to regulate several transcription factors during the developmental transition from larvae to adult worms (Grosshans et al., 2005). Based on the homologous relationship between hco-miR-5991 and cel-let-7 (Winter et al., 2012), and our analytic result that hco-miR-5991 could interact with 167 circRNAs, we speculate that some circRNAs might participate in H. contortus larval development via the let-7 pathway. In addition, the up-regulation of C. elegans miRNA lin-4 reduced the activity of its target, the transcription factor lin-14, which then significantly extended the lifespan of this nematode (Boehm and Slack, 2005). The lin-4-lin-14 pair conferred life span extension through the insulin/insulin-like growth factor-1 pathway (Boehm and Slack, 2005). The mature sequence of cel-lin-4 has high homology with hco-lin-4, which was predicted to be bound by 36 circRNAs, suggesting that some circRNAs could also affect nematode developmental processes through this pathway in H. contortus. Also in this parasite, a previous study reported that hco-miR-228 and hco-miR-235 were enriched in L3 (Marks et al., 2019) and our analysis found that hco-miR-228 and hco-miR-235 combined with 78 and 55 circRNAs, respectively, suggesting that some circRNAs may be involved in L3’s development.

Except for being involved in the developmental regulation of nematodes, miRNAs could play roles in drug resistance. The C. elegans cel-miR-1 negatively regulated the expression of two nicotinic acetylcholine receptor (nAChR) subunits, UNC-29 and UNC-63, which are targets of levamisole (Simon et al., 2008). Intriguingly, the expression levels of both subunits were increased in C. elegans miR-1 mutant strains, whereas the muscle sensitivity to acetylcholine and levamisole was decreased. As hco-miR-1 was homologous to cel-miR-1 and could bind to 19 circRNAs, we surmise that circRNAs may also play functions in drug resistance in H. contortus.

Although whether circRNAs can act as miRNA sponges has not been experimentally verified in H. contortus, there should be a high probability that this mechanism exists in nematodes, like miRNA sponges in animals that regulate gene expression at the epigenetic level (Luo et al., 2021; Zhou et al., 2021). Certainly, the potential interplay between circRNA and miRNA in parasitic nematodes needs to be further studied in defined biological validation experiments in the future (Li et al., 2018b). Due to the lack of transgenic methods and the restriction of parasitic nematode culture conditions in vitro, cross-linking immunoprecipitation (CLIP) approaches can be applied to recognize circRNA-miRNA interactions using the antibody to the Argonaute (AGO) protein. Certainly, functions of more and more circRNAs will be elucidated in the future with the advances of new technologies.

Conclusion

In summary, the general features and expression profiles of circRNAs at critical stages of H. contortus were characterized and bioinformatic analysis revealed that circRNAs could have multiple biological functions during larval developmental processes. Abundant circRNAs contained one or more miRNA potential binding sites, suggesting that they may act as new post-transcriptional regulators involved in the development of H. contortus. Therefore, this study provided a basis for elucidating the dynamic regulation of circRNAs in H. contortus developmental processes. In addition, parasitic worms had relatively complex gene expression regulation system, so it was of great significance to study non-coding RNAs for revealing the physiological mechanism of parasites and the molecular mechanism of parasitic diseases.

Data Availability Statement

The data presented in the study are deposited in the NCBI repository, the accession number: PRJNA756739.

Ethics Statement

Experimental animals used in the project were treated in strict accordance with the Rules for Animal Ethics and Experimentation in the People’s Republic of China. The maintenance of goats was in accordance with protocols approved by The Scientific Ethics Committee of Huazhong Agricultural University (permit HZAUGO-2016-007).

Author Contributions

MH and CZ contributed to conception and design of the project. YZ, SW, ZW, CW, and FL contributed to animal experiments. CZ and WT performed the experiments and analyzed data. CZ drafted the manuscript, and MH edited and revised the manuscript. All the authors contributed to the article and approved the submitted version.

Funding

This work was funded by the National Natural Science Foundation of China (NSFC) (Grant no. 31872462 and 321719037) to MH.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcimb.2021.764089/full#supplementary-material

Supplementary Figure 1 | Verification of head-to-tail back-splicing site. Black inverted triangles represented the back-spliced junction loci. Marker: 2000 bp ladder DNA marker.

Supplementary Table 1 | Detailed information of convergent and divergent primers used for validating the back-splicing junction sites of circRNAs in H. contortus.

Supplementary Table 2 | Information on divergent primers for validating the RNA-Seq reliability and the resistance to RNase R digestion of circRNAs in H. contortus.

Supplementary Table 3 | Summary of sequencing results in H. contortus generated from nine libraries.

Supplementary Table 4 | Information of all novel circRNAs identified in H. contortus.

Supplementary Table 5 | The prediction of interaction network for H. contortus circRNA-miRNA.

Supplementary Table 6 | Information on ten specific binding sites between hco_circ_00015227 and hco-miR-307.

Supplementary Table 7 | Primer sequences and product lengths of circRNAs for alternative event.

References

Boehm, M., Slack, F. (2005). A Developmental Timing microRNA and its Target Regulate Life Span in C. Elegans. Science 310 (5756), 1954–1957. doi: 10.1126/science.1115596

PubMed Abstract | CrossRef Full Text | Google Scholar

Broadbent, K. M., Broadbent, J. C., Ribacke, U., Wirth, D., Rinn, J. L., Sabeti, P. C. (2015). Strand-Specific RNA Sequencing in Plasmodium Falciparum Malaria Identifies Developmentally Regulated Long non-Coding RNA and Circular RNA. BMC Genomics 16, 454. doi: 10.1186/s12864-015-1603-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Capel, B., Swain, A., Nicolis, S., Hacker, A., Walter, M., Koopman, P., et al. (1993). Circular Transcripts of the Testis-Determining Gene Sry in Adult Mouse Testis. Cell 73 (5), 1019–1030. doi: 10.1016/0092-8674(93)90279-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, G., Cui, J., Wang, L., Zhu, Y., Lu, Z., Jin, B. (2017). Genome-Wide Identification of Circular RNAs in Arabidopsis Thaliana. Front. Plant Sci. 8:1678. doi: 10.3389/fpls.2017.01678

PubMed Abstract | CrossRef Full Text | Google Scholar

Cortes-Lopez, M., Gruner, M. R., Cooper, D. A., Gruner, H. N., Voda, A. I., van der Linden, A. M., et al. (2018). Global Accumulation of circRNAs During Aging in Caenorhabditis Elegans. BMC Genomics 19 (1):8. doi: 10.1186/s12864-017-4386-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Danan, M., Schwartz, S., Edelheit, S., Sorek, R. (2012). Transcriptome-Wide Discovery of Circular RNAs in Archaea. Nucleic Acids Res. 40 (7), 3131–3142. doi: 10.1093/nar/gkr1009

PubMed Abstract | CrossRef Full Text | Google Scholar

Di, W., Gasser, R. B., He, L., Li, F., Liu, X., Zhou, C., et al. (2020). A Serine/Threonine-Specific Protein Kinase of Haemonchus Contortus With a Role in the Development. FASEB J. 34 (2), 2075–2086. doi: 10.1096/fj.201900888RR

PubMed Abstract | CrossRef Full Text | Google Scholar

Du, W. W., Yang, W., Liu, E., Yang, Z., Dhaliwal, P., Yang, B. B. (2016). Foxo3 Circular RNA Retards Cell Cycle Progression via Forming Ternary Complexes With P21 and CDK2. Nucleic Acids Res. 44 (6), 2846–2858. doi: 10.1093/nar/gkw027

PubMed Abstract | CrossRef Full Text | Google Scholar

Enright, A. J., John, B., Gaul, U., Tuschl, T., Sander, C., Marks, D. S. (2003). MicroRNA Targets in Drosophila. Genome Biol. 5 (1), R1. doi: 10.1186/gb-2003-5-1-r1

PubMed Abstract | CrossRef Full Text | Google Scholar

Feng, J., Chen, K., Dong, X., Xu, X., Jin, Y., Zhang, X., et al. (2019). Genome-Wide Identification of Cancer-Specific Alternative Splicing in circRNA. Mol. Cancer 18 (1), 35. doi: 10.1186/s12943-019-0996-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, Y., Wang, J., Zheng, Y., Zhang, J., Chen, S., Zhao, F. (2016). Comprehensive Identification of Internal Structure and Alternative Splicing Events in Circular RNAs. Nat. Commun. 7, 12060. doi: 10.1038/ncomms12060

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, Y., Zhang, J., Zhao, F. (2018). Circular RNA Identification Based on Multiple Seed Matching. Brief Bioinform. 19 (5), 803–810. doi: 10.1093/bib/bbx014

PubMed Abstract | CrossRef Full Text | Google Scholar

Grabowski, P. J., Zaug, A. J., Cech, T. R. (1981). The Intervening Sequence of the Ribosomal RNA Precursor Is Converted to a Circular RNA in Isolated Nuclei of Tetrahymena. Cell 23 (2), 467–476. doi: 10.1016/0092-8674(81)90142-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Grosshans, H., Johnson, T., Reinert, K. L., Gerstein, M., Slack, F. J. (2005). The Temporal Patterning microRNA Let-7 Regulates Several Transcription Factors at the Larval to Adult Transition in C. Elegans. Dev. Cell 8 (3), 321–330. doi: 10.1016/j.devcel.2004.12.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, J. U., Agarwal, V., Guo, H., Bartel, D. P. (2014). Expanded Identification and Characterization of Mammalian Circular RNAs. Genome Biol. 15 (7):409. doi: 10.1186/s13059-014-0409-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Hansen, T. B., Jensen, T. I., Clausen, B. H., Bramsen, J. B., Finsen, B., Damgaard, C. K., et al. (2013). Natural RNA Circles Function as Efficient microRNA Sponges. Nature 495 (7441), 384–388. doi: 10.1038/nature11993

PubMed Abstract | CrossRef Full Text | Google Scholar

He, L., Gasser, R. B., Li, T., Di, W., Li, F., Zhang, H., et al. (2019). A TGF-β Type II Receptor That Associates With Developmental Transition in Haemonchus Contortus In Vitro. PloS Negl. Trop. Dis. 13 (12), e0007913. doi: 10.1371/journal.pntd.0007913

PubMed Abstract | CrossRef Full Text | Google Scholar

Hsu, M. T., Coca-Prados, M. (1979). Electron Microscopic Evidence for the Circular Form of RNA in the Cytoplasm of Eukaryotic Cells. Nature 280 (5720), 339–340. doi: 10.1038/280339a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Ivanov, A., Memczak, S., Wyler, E., Torti, F., Porath, H. T., Orejuela, M. R., et al. (2015). Analysis of Intron Sequences Reveals Hallmarks of Circular RNA Biogenesis in Animals. Cell Rep. 10 (2), 170–177. doi: 10.1016/j.celrep.2014.12.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Kanehisa, M., Furumichi, M., Tanabe, M., Sato, Y., Morishima, K. (2017). KEGG: New Perspectives on Genomes, Pathways, Diseases and Drugs. Nucleic Acids Res. 45 (D1), D353–d361. doi: 10.1093/nar/gkw1092

PubMed Abstract | CrossRef Full Text | Google Scholar

Kristensen, L. S., Andersen, M. S., Stagsted, L. V. W., Ebbesen, K. K., Hansen, T. B., Kjems, J. (2019). The Biogenesis, Biology and Characterization of Circular RNAs. Nat. Rev. Genet. 20 (11), 675–691. doi: 10.1038/s41576-019-0158-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Laing, R., Kikuchi, T., Martinelli, A., Tsai, I. J., Beech, R. N., Redman, E., et al. (2013). The Genome and Transcriptome of Haemonchus Contortus, a Key Model Parasite for Drug and Vaccine Discovery. Genome Biol. 14 (8), R88. doi: 10.1186/gb-2013-14-8-r88

PubMed Abstract | CrossRef Full Text | Google Scholar

Langmead, B., Trapnell, C., Pop, M., Salzberg, S. L. (2009). Ultrafast and Memory-Efficient Alignment of Short DNA Sequences to the Human Genome. Genome Biol. 10 (3), R25. doi: 10.1186/gb-2009-10-3-r25

PubMed Abstract | CrossRef Full Text | Google Scholar

Legnini, I., Di Timoteo, G., Rossi, F., Morlando, M., Briganti, F., Sthandier, O., et al. (2017). Circ-ZNF609 Is a Circular RNA That can be Translated and Functions in Myogenesis. Mol. Cell 66 (1), 22–37.e29. doi: 10.1016/j.molcel.2017.02.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Liang, G., Yang, Y., Niu, G., Tang, Z., Li, K. (2017). Genome-Wide Profiling of Sus Scrofa Circular RNAs Across Nine Organs and Three Developmental Stages. DNA Res. 24 (5), 523–535. doi: 10.1093/dnares/dsx022

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y., Chen, B., Huang, S. (2018b). Identification of circRNAs for miRNA Targets by Argonaute2 RNA Immunoprecipitation and Luciferase Screening Assays. Methods Mol. Biol. 1724, 209–218. doi: 10.1007/978-1-4939-7562-4_17

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, H., Durbin, R. (2009). Fast and Accurate Short Read Alignment With Burrows-Wheeler Transform. Bioinformatics 25 (14), 1754–1760. doi: 10.1093/bioinformatics/btp324

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, F. C., Gasser, R. B., Lok, J. B., Korhonen, P. K., He, L., Di, W. D., et al. (2016). Molecular Characterization of the Haemonchus Contortus Phosphoinositide-Dependent Protein Kinase-1 Gene (Hc-Pdk-1). Parasit Vectors 9, 65. doi: 10.1186/s13071-016-1351-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, F. C., Gasser, R. B., Lok, J. B., Korhonen, P. K., Wang, Y. F., Yin, F. Y., et al. (2014). Exploring the Role of Two Interacting Phosphoinositide 3-Kinases of Haemonchus Contortus. Parasit Vectors 7, 498. doi: 10.1186/s13071-014-0498-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Z., Huang, C., Bao, C., Chen, L., Lin, M., Wang, X., et al. (2015). Exon-Intron Circular RNAs Regulate Transcription in the Nucleus. Nat. Struct. Mol. Biol. 22 (3), 256–264. doi: 10.1038/nsmb.2959

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, H., Hu, Y., Yin, J., Yan, X. Y., Chen, W. J., Jiang, C. Y., et al. (2019). Profiles Analysis Reveals Circular RNAs Involving Zebrafish Physiological Development. J. Cell Physiol. 234, 15922–15933. doi: 10.1002/jcp.28250

CrossRef Full Text | Google Scholar

Li, X., Yang, L., Chen, L. L. (2018a). The Biogenesis, Functions, and Challenges of Circular RNAs. Mol. Cell 71 (3), 428–442. doi: 10.1016/j.molcel.2018.06.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Love, M. I., Huber, W., Anders, S. (2014). Moderated Estimation of Fold Change and Dispersion for RNA-Seq Data With Deseq2. Genome Biol. 15 (12), 550. doi: 10.1186/s13059-014-0550-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, T., Cui, L., Zhou, Y., Zhu, C., Fan, D., Gong, H., et al. (2015). Transcriptome-Wide Investigation of Circular RNAs in Rice. RNA 21 (12), 2076–2087. doi: 10.1261/rna.052282.115

PubMed Abstract | CrossRef Full Text | Google Scholar

Luo, L. P., Miao, P. P., Ming, Y., Tao, J., Shen, H. C. (2021). Circ-ZFR Promotes Progression of Bladder Cancer by Upregulating WNT5A Via Sponging miR-545 and miR-1270. Front. Oncol. 10, 596623. doi: 10.3389/fonc.2020.596623

PubMed Abstract | CrossRef Full Text | Google Scholar

Mani, V., Assefa, A. D., Hahn, B.-S. (2021). Transcriptome Analysis and miRNA Target Profiling at Various Stages of Root-Knot Nematode Development for Identification of Potential Regulatory Networks. Int. J. Mol. Sci. 22 (14), 7442. doi: 10.3390/ijms22147442

PubMed Abstract | CrossRef Full Text | Google Scholar

Marks, N. D., Winter, A. D., Gu, H. Y., Maitland, K., Gillan, V., Ambroz, M., et al. (2019). Profiling microRNAs Through Development of the Parasitic Nematode Haemonchus Identifies Nematode-Specific miRNAs That Suppress Larval Development. Sci. Rep. 9 (1), 17594. doi: 10.1038/s41598-019-54154-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, G., Wang, T., Korhonen, P. K., Ang, C. S., Williamson, N. A., Young, N. D., et al. (2018). Molecular Alterations During Larval Development of Haemonchus Contortus In Vitro Are Under Tight Post-Transcriptional Control. Int. J. Parasitol 48 (9-10), 763–772. doi: 10.1016/j.ijpara.2018.03.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Memczak, S., Jens, M., Elefsinioti, A., Torti, F., Krueger, J., Rybak, A., et al. (2013). Circular RNAs Are a Large Class of Animal RNAs With Regulatory Potency. Nature 495 (7441), 333–338. doi: 10.1038/nature11928

PubMed Abstract | CrossRef Full Text | Google Scholar

Nicolet, B. P., Engels, S., Aglialoro, F., van den Akker, E., von Lindern, M., Wolkers, M. C. (2018). Circular RNA Expression in Human Hematopoietic Cells Is Widespread and Cell-Type Specific. Nucleic Acids Res. 46 (16), 8168–8180. doi: 10.1093/nar/gky721

PubMed Abstract | CrossRef Full Text | Google Scholar

Salzman, J., Chen, R. E., Olsen, M. N., Wang, P. L., Brown, P. O. (2013). Cell-Type Specific Features of Circular RNA Expression. PloS Genet. 9 (9), e1003777. doi: 10.1371/journal.pgen.1003777

PubMed Abstract | CrossRef Full Text | Google Scholar

Schwarz, E. M., Korhonen, P. K., Campbell, B. E., Young, N. D., Jex, A. R., Jabbar, A., et al. (2013). The Genome and Developmental Transcriptome of the Strongylid Nematode Haemonchus Contortus. Genome Biol. 14 (8), R89. doi: 10.1186/gb-2013-14-8-r89

PubMed Abstract | CrossRef Full Text | Google Scholar

Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: A Software Environment for Integrated Models of Biomolecular Interaction Networks. Genome Res. 13 (11), 2498–2504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

Simon, D. J., Madison, J. M., Conery, A. L., Thompson-Peer, K. L., Soskis, M., Ruvkun, G. B., et al. (2008). The microRNA miR-1 Regulates a MEF-2-Dependent Retrograde Signal at Neuromuscular Junctions. Cell 133 (5), 903–915. doi: 10.1016/j.cell.2008.04.035

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang, X., Wen, X. M., Li, Z. Y., Wen, D. X., Lin, L., Liu, J. Q., et al. (2021). Hsa_circ_0102171 Aggravates the Progression of Cervical Cancer Through Targeting miR-4465/CREBRF Axis. J. Cell Physiol. 236 (7), 4973–4984. doi: 10.1002/jcp.30210

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Z. P., Liu, Y. F., Li, D. W., Li, L., Zhang, Q., Wang, S. B., et al. (2017). Identification of Circular RNAs in Kiwifruit and Their Species-Specific Response to Bacterial Canker Pathogen Invasion. Front. Plant Sci. 8, 413. doi: 10.3389/fpls.2017.00413

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Yang, M., Wei, S., Qin, F., Zhao, H., Suo, B. (2016). Identification of Circular RNAs and Their Targets in Leaves of L. Under Dehydration Stress. Front. Plant Sci. 7, 2024. doi: 10.3389/fpls.2016.02024

PubMed Abstract | CrossRef Full Text | Google Scholar

Westholm, J. O., Miura, P., Olson, S., Shenker, S., Joseph, B., Sanfilippo, P., et al. (2014). Genome-Wide Analysis of Drosophila Circular RNAs Reveals Their Structural and Sequence Properties and Age-Dependent Neural Accumulation. Cell Rep. 9 (5), 1966–1980. doi: 10.1016/j.celrep.2014.10.062

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilusz, J. E. (2017). Circular RNAs: Unexpected Outputs of Many Protein-Coding Genes. RNA Biol. 14 (8), 1007–1017. doi: 10.1080/15476286.2016.1227905

PubMed Abstract | CrossRef Full Text | Google Scholar

Winter, A. D., Weir, W., Hunt, M., Berriman, M., Gilleard, J. S., Devaney, E., et al. (2012). Diversity in Parasitic Nematode Genomes: The microRNAs of Brugia Pahangi and Haemonchus Contortus Are Largely Novel. BMC Genomics 13, 4. doi: 10.1186/1471-2164-13-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Xin, R., Qu, D. H., Xu, H. Y., Chen, D. W. (2021). Circ_001504 Promotes the Development of Renal Cell Carcinoma by Sponging microRNA-149 to Increase NUCB2. Cancer Gene Ther. 28 (6), 667–678. doi: 10.1038/s41417-020-00247-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, T., Wu, J., Han, P., Zhao, Z., Song, X. (2017). Circular RNA Expression Profiles and Features in Human Tissues: A Study Using RNA-Seq Data. BMC Genomics 18 (Suppl 6), 680. doi: 10.1186/s12864-017-4029-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, Y., Zhang, Y., Chen, B., Ding, L., Mu, Z., Li, Y. (2019). Elevation of Circular RNA Circ-POSTN Facilitates Cell Growth and Invasion by Sponging miR-1205 in Glioma. J. Cell Biochem. 120 (10), 16567–16574. doi: 10.1002/jcb.28916

PubMed Abstract | CrossRef Full Text | Google Scholar

Young, M. D., Wakefield, M. J., Smyth, G. K., Oshlack, A. (2010). Gene Ontology Analysis for RNA-Seq: Accounting for Selection Bias. Genome Biol. 11 (2), R14. doi: 10.1186/gb-2010-11-2-r14

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, G. C., Wang, L. G., Han, Y. Y., He, Q. Y. (2012). "Clusterprofiler: An R Package for Comparing Biological Themes Among Gene Clusters.". OMICS 16 (5), 284–287. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeng, X., Lin, W., Guo, M., Zou, Q. (2017). A Comprehensive Overview and Evaluation of Circular RNA Detection Tools. PloS Comput. Biol. 13 (6), e1005420. doi: 10.1371/journal.pcbi.1005420

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, P. J., Li, S. D., Chen, M. (2020). Characterization and Function of Circular RNAs in Plants. Front. Mol. Biosci. 7, 91. doi: 10.3389/fmolb.2020.00091

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Zhang, W., Dong, M. (2018). The miR-58 microRNA Family Is Regulated by Insulin Signaling and Contributes to Lifespan Regulation in Caenorhabditis Elegans. Sci. China Life Sci. 61 (9), 1060–1070. doi: 10.1007/s11427-018-9308-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, L., Chen, J., Li, Z., Li, X., Hu, X., Huang, Y., et al. (2010). Integrated Profiling of microRNAs and mRNAs: microRNAs Located on Xq27.3 Associate With Clear Cell Renal Cell Carcinoma. PloS One 5 (12), e15224. doi: 10.1371/journal.pone.0015224

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, Z. B., Ma, J., Lu, J. J., Chen, A. M., Zhu, L. (2021). Circular RNA CircCDH13 Contributes to the Pathogenesis of Osteoarthritis via CircCDH13/miR-296-3p/PTEN Axis. J. Cell Physiol. 236 (5), 3521–3535. doi: 10.1002/jcp.30091

PubMed Abstract | CrossRef Full Text | Google Scholar