Expression Patterns of Circular RNAs from Primary Kinase Transcripts in the Mammary Glands of Lactating Rats.

PURPOSE: Circular RNAs (circRNAs), a novel class of RNAs, perform important functions in biological processes. However, the role of circRNAs in the mammary gland remains unknown. The present study is aimed at identifying and characterizing the circRNAs expressed in the mammary gland of lactating rats.
METHODS: Deep sequencing of RNase R-enriched rat lactating mammary gland samples was performed and circRNAs were predicted using a previously reported computational pipeline. Gene ontology terms of circRNA-producing genes were also analyzed.
RESULTS: A total of 6,824 and 4,523 circRNAs were identified from rat mammary glands at two different lactation stages. Numerous circRNAs were specifically expressed at different lactation stages, and only 1,314 circRNAs were detected at both lactation stages. The majority of the candidate circRNAs map to noncoding intronic and intergenic regions. The results demonstrate a circular preference or specificity of some genes. DAVID analysis revealed an enrichment of protein kinases and related proteins among the set of genes encoding circRNAs. Interestingly, four protein-coding genes (Rev3l, IGSF11, MAML2, and LPP) that also transcribe high levels of circRNAs have been reported to be involved in cancer.
CONCLUSION: Our findings provide the basis for comparison between breast cancer profiles and for selecting representative circRNA candidates for future functional characterization in breast development and breast cancer.

INTRODUCTION

RNA functions as an intermediate molecule in the transfer of genetic information from DNA to protein. Noncoding RNAs (ncRNAs) are defined as RNAs that do not code proteins; they include all traditional classes of RNAs (except for mRNAs) such as rRNAs, tRNAs, and snRNAs, as well as the more recently discovered microRNAs (miRNAs) and long noncoding RNAs (lncRNAs). With the advent of robust high throughput sequencing technologies and advances in bioinformatics, thousands of long and short ncRNAs have been discovered in the past decade in prokaryotes, eukaryotes, and viruses. These ncRNAs encompass transcripts of different sizes, from processed small RNAs ~20-30 bases in length to lncRNAs that are hundreds or even thousands of bases long [1]. These recently discovered ncRNAs represent a new paradigm of gene regulation. In particular, many miRNAs and lncRNAs have been reported to regulate both normal development and tumorigenesis in the mammary gland [2].

Circular RNA (circRNA) is a recently discovered ncRNA. CircRNA isoforms have largely gone unnoticed since their discovery in 1979, with some notable exceptions [3]. Recent studies by two independent teams demonstrated that a human circRNA running antisense to the cerebellar degenerationrelated protein 1 (CDR1) locus functions as a sponge for miR-7 [45]. Interestingly, miR-7 inhibits the epithelial-mesenchymal transition and metastasis of breast cancer, and sensitizes cells to DNA damage [67], which indicates that circRNAs may play profound roles in breast cancer. A circRNA transcribed from the Sry gene was also reported to act as a miR-138 sponge [8]. These two examples represent a strong paradigm for circRNA function, and have therefore attracted attention towards circRNAs. Growing evidence indicates that circRNAs, which are expressed from a large fraction of human, mouse, zebrafish, Caenorhabditis elegans, and fruit fly genes, constitute the major RNA isoform from these genes [9101112]. However, no information on mammary gland circRNAs has been reported.

Breastfeeding is not only beneficial to infants, but also provides many health benefits to mothers [13]. The maternal health benefits include decreased postpartum bleeding, more rapid uterine involution, decreased menstrual blood loss, increased child spacing, earlier return to prepregnancy weight, and possible decreased risk of hip fractures and osteoporosis in the postmenopausal period [13]. One of the most notable maternal health benefits is a decreased risk of breast cancer, which has been demonstrated by several groups using different methods [14]. These findings suggest that breastfeeding affects biological processes in the mammary gland. However, the molecular mechanisms underlying the role of breastfeeding in breast cancer prevention remain unclear. Accordingly, an understanding of the molecular landscape of the normal breast during lactation is critical for breast cancer prevention.

In the present study, we investigate the circRNA expression pattern in the lactating mammary gland to provide a basis for comparison between breast cancer profiles, and for selection of representative circRNA candidates for future functional characterization in breast development and breast cancer.

METHODS

Animal material

Mammary gland tissue samples from ten Sprague Dawley rats (Rattus norvegicus) fifth parturition were obtained from Vital River Laboratories (Beijing, China). Lactating mammary glands were harvested from five rats at day 1 and 7 postpartum and frozen in liquid nitrogen for subsequent extraction of RNA. All samples were obtained in accordance with the experimental procedures guideline of Jiangsu Normal University, and all experimental protocols were approved by the Animal Ethics Committee of Jiangsu Normal University.

RNA isolation and RNase R enrichment

Total RNA was extracted using Trizol (Life Technologies, Carlsbad, USA), and treated with DNase I. The RNase R enrichment protocol was modified from a previous report [11]. Each total RNA sample (20 µg) was depleted for ribosomal RNAs using the RiboMinus kit (Life Technologies), and then subjected to RNase R treatment. RiboMinus RNA was incubated for 15 minutes at 37℃, with or without (control) 3 U/µg of RNase R (Epicentre Biotechnologies). RNA was subsequently purified by phenol-chloroform extraction, and divided into two parts, one of which was used to prepare libraries for sequencing, using Illumina's whole-transcriptome RNA sequencing protocols (Illumina, San Diego, USA). The Illumina Hiseq2000 System was used to perform sequencing runs for RNA sequencing. The other was reverse transcribed using the SuperScript III FirstStrand Synthesis System (Life Technologies) with random hexamers, according to the manufacturer's instructions.

Computational pipeline for prediction of circRNAs

The rat reference genome was downloaded from the National Center for Biotechnology Information (NCBI) genome browser (ftp://ftp.ncbi.nlm.nih.gov/genomes/R_norvegicus/). Sequencing reads from each pooled sample were mapped independently using the short-read mapper Bowtie 2 [15]. Reads aligning contiguously and full length to the genome were discarded. We extracted 20-mers from both ends and aligned them independently to find unique anchor positions within spliced exons from the remaining reads. Anchors that aligned in the reverse orientation (head-to-tail) indicated circRNA splicing. We extended the anchor alignments such that the complete read aligns and the breakpoints were flanked by GU/AG splice sites. Ambiguous breakpoints were discarded. More information can be found in a previous report [5].

Gene ontology terms of circRNA-producing genes

The DAVID tool (National Institute of Allergy and Infectious Diseases, Bethesda, USA) was used for enrichment of gene ontology (GO) terms of circRNA-producing genes.

Detecting putative miRNA binding sites

FASTA files of rat miRNAs were obtained from miRBase release 20.0 (http://www.mirbase.org/). Only mature miRNAs were considered for seed analysis. The miRNAs were aligned with circRNAs. A putative target site of a miRNA is a 6-nucleotide-long sequence in the genome that represents the reverse complement of nucleotides 2-7 of the mature miRNA sequence.

Modeling of linear to circular expression

In order to test the global relationship between the circRNAs and linear RNAs, we fit a Poisson model per gene modeling circular expression by reads per kilobase of exon model per million mapped reads (RPKM) to poly(A) gene expression using SPSS version 19.0 (IBM Corp., Armonk, USA).

Quantitative reverse-transcription polymerase chain reaction analysis

cDNAs of mock and RNase R-treated RNA were diluted with water and used as a template for polymerase chain reaction (PCR). In order to validate the unpolyadenylated status of circRNAs, the total RNAs were reverse-transcribed with oligo-dTs and random hexamer primers, and the cDNAs were subjected to qPCR. Standard PCR was performed using Taq DNA polymerase (New England Biolabs, Beverly, USA), and quantitative PCR was performed using SYBR Green Mix (Roche, Basel, Switzerland) on a StepOne plus Real-Time PCR System. qPCR Ct values were calculated automatically using the manufacturer's software. The PCR primers are listed in Supplementary Table 1 (available online). PCR products were directly Sanger-sequenced in both directions using amplification primers.

RESULTS

CircRNA transcriptome composition

We obtained 11,853,653 and 11,970,307 clean sequence reads of 86 nucleotides from the two RNase R-treated RiboMinus mammary gland RNA pools. Nearly half (45%) of the clean reads were predicted as circRNAs, and most of the circRNAs were transcribed from gene-containing regions covering 1,569 genes. These data were consistent with findings from previous studies that identified circRNAs from mainly exons and introns of gene-containing regions (Figure 1) [34591011].

Figure 1

The distribution of circular RNAs in the genome. (A) Rat mammary gland library at day 1 postpartum. (B) Rat mammary gland library at day 7 postpartum.

CDS=coding sequences.

CircRNA expression patterns

We detected 6,824 and 4,523 circRNAs from lactating mammary glands at day 1 (D1) and 7 (D7) postpartum, respectively. The length of circRNAs varied from 66 nt to 1,955 nt. Numerous circRNAs seemed to be specifically expressed at different lactation stages, as only 1,314 circRNAs were detected at both lactation stages.

The circularization of circRNAs and their relative expression levels, as obtained from RNA-sequencing, were validated. Inward and outward facing primers were designed against 11 randomly selected representative transcripts. Each primer pair amplified a single, distinct product of the expected size and sequence. RNase R resistance of all 11 novel backsplice events was apparent following RNase R treatment, whereas the abundance of linear RNAs was found to decrease (Table 1). Moreover, oligo-dT priming for reverse transcription resulted in significantly lower levels of backspliced products relative to poly(A)-containing transcript levels, indicating that these species were not polyadenylated (Table 1). The difference in the expression of the 11 circRNAs between the two time points was also confirmed by qPCR (Supplementary Figure 1).

Table 1

Quantitative polymerase chain reaction validation of circular RNAs by different reverse transcription primers and RNase R treatment

circRNAs	Reverse transcription treatment*		RNase R treatment†
	Circular	Linear	Circular	Linear
circRNA_555	-4.41 (0.29)	-0.08 (0.08)	-0.74 (0.11)	3.99 (0.04)
circRNA_127	-3.57 (2.46)	-2.34 (1.87)	-0.36 (0.56)	5.25 (0.43)
circRNA_3385	0.43 (1.45)	-0.90 (0.31)	0.06 (0.10)	6.16 (0.42)
circRNA_659	-4.55 (0.57)	1.06 (0.01)	0.61 (0.62)	9.10 (0.42)
circRNA_10008	-7.34 (1.51)	-2.16 (0.07)	0.44 (0.33)	4.16 (0.56)
circRNA_537	-5.45 (0.32)	-6.20 (0.22)	0.38 (0.53)	2.06 (0.11)
circRNA_388	-5.24 (1.40)	-2.93 (1.59)	-0.19 (0.56)	2.95 (0.11)
circRNA_15703	-3.84 (0.69)	-1.98 (0.05)	1.51 (0.54)	5.77 (0.39)
circRNA_7070	-2.79 (0.61)	-2.69 (0.24)	-2.48 (0.58)	3.19 (0.22)
circRNA_36159	-0.32 (0.68)	-1.99 (0.04)	0.23 (0.34)	7.98 (0.09)
circRNA_4815	-3.31 (1.51)	-4.31 (0.29)	-1.37 (0.74)	1.96 (0.52)

circRNAs=circular RNAs.

*The values were ΔCT values between random hexamer primer generated templates and oligo-dT generated templates with standard deviation in bracket; †The values were ΔCT values between templates from RNase R treated RNA and mock treated RNA with standard deviation in bracket.

In order to exclude the possibility that circRNAs might be the result of the background "noise" level of dysfunctional splicing, we analyzed the relationship between linear RNA isoform expression from a given gene and the probability of detecting a circRNA isoform from that gene. We did not find any correlation between the expressions of the two RNA isoforms (Figure 2).

Figure 2

Quantitative regulation of circular to linear isoform. No genome-wide trend between circular expression and linear transcript expression as measured by log PolyA reads per kilobase of exon model per million mapped reads (RPKM).

Y axis=RPKM for circular RNA; X axis=log PolyA RPKM of linear isoform.

Many circRNAs are derived from kinase genes

In order to investigate the function of circRNAs, the GO term enrichment of the set circRNA-producing genes was analyzed using the DAVID tool. The top five annotation clusters, summarized in Table 2, involved nucleotide and ATP binding, nuclear lumen, GTPase regulator activity, ion binding and protein kinase activity.

Table 2

Gene ontology term enrichment of circular RNA generating genes

Category (enrichment score)	Description	Count	Fold enrichment	False discovery rate
Nucleotide and adenosine-triphosphate binding (6.60)	Purine nucleotide binding	87	1.78	7.79E-05
	Purine ribonucleotide binding	84	1.80	8.05E-05
	Ribonucleotide binding	84	1.80	8.51E-05
	Nucleoside binding	77	1.86	9.54E-05
	Purine nucleoside binding	76	1.85	1.46E-04
	Nucleotide binding	98	1.68	1.53E-04
	Adenyl ribonucleotide binding	72	1.89	1.55E-04
	Adenyl nucleotide binding	75	1.86	1.56E-04
	ATP binding	71	1.91	1.57E-04
Nuclear lumen (4.88)	Nucleoplasm	40	2.23	4.18E-03
	Membrane-enclosed lumen	59	1.82	7.17E-03
	Nuclear lumen	47	1.97	1.38E-02
	Intracellular organelle lumen	55	1.81	2.08E-02
	Organelle lumen	56	1.78	2.82E-02
	Nucleoplasm part	26	2.34	1.62E-01
GTPase regulator activity (3.63)	Small GTPase regulator activity	19	3.53	1.05E-02
	Regulation of small GTPase mediated signal transduction	19	3.37	2.35E-02
	GTPase regulator activity	22	2.76	7.08E-02
	Nucleoside-triphosphatase regulator activity	22	2.70	9.83E-02
	Regulation of Ras protein signal transduction	15	3.38	2.39E-01
	Enzyme activator activity	16	2.86	7.12E-01
	GTPase activator activity	12	3.37	1.24E+00
	Regulation of GTPase activity	9	3.83	3.93E+00
	Regulation of Ras GTPase activity	8	4.26	4.38E+00
Ion binding (3.53)	Zinc ion binding	71	1.63	4.83E-02
	Transition metal ion binding	84	1.50	1.36E-01
	Metal ion binding	116	1.36	2.46E-01
	Cation binding	116	1.34	4.35E-01
	Ion binding	116	1.32	8.71E-01
Protein kinase activity (3.49)	Phosphate metabolic process	54	2.17	1.92E-04
	Phosphorus metabolic process	54	2.16	2.04E-04
	Phosphorylation	46	2.19	1.54E-03
	Protein amino acid phosphorylation	41	2.26	3.32E-03
	Protein kinase activity	38	2.26	7.63E-03
	Protein serine/threonine kinase activity	30	2.55	9.75E-03

ATP=adenosine triphosphatase; GTPase=guanosine triphosphatase.

Docking sites for miRNAs in circRNAs

In the present study, 716 out of 728 mature rat miRNAs were predicted to have docking sites in the circRNAs detected in the rat mammary gland. Only a few circRNAs were predicted to have many miRNA-binding sites, which indicated that circularization was not a global feature under selection in the sequence of the thousands of circRNAs profiled.

DISCUSSION

Previous studies have identified circRNAs mainly from exons and introns in gene-containing regions [34591011]. In the present study, circRNAs were not only identified from gene containing regions, but also from intergenic regions. The circRNAs distributed in noncoding regions accounted for ~70% of the total circRNAs, including those from intronic and intergenic regions, of which intronic circRNAs accounted for more than 40%. This is similar to the expression pattern of other ncRNAs [16].

The top five annotation clusters involved nucleotide and ATP binding, nuclear lumen, GTPase regulator activity, ion binding and protein kinase activity. As the activation of virtually all cell-surface receptors leads directly or indirectly to changes in protein phosphorylation via activation of protein kinases, the circRNAs derived from protein kinase-encoding genes may also be functional [17]. We also analyzed mouse and human circRNA data from several tissues [51011]. Notably, DAVID analysis revealed an enrichment of protein kinases and related proteins among the set of genes producing circRNAs in the tissues of various species. This finding further verified that circRNAs show a high degree of evolutionary conservation, suggesting that they might have important biological functions.

In the present study, most circRNAs were transcribed from protein-coding gene regions, which indicated a circular preference or specificity of some genes. The set of genes producing circRNAs in the tissues of various mammal species were all enriched for protein kinases and related proteins [511]. As shown in Table 3, seven of the 10 most highly expressed circRNA families were transcribed from seven different protein-coding gene regions. Interestingly, four of the protein-coding genes (Rev3l, IGSF11, MAML2, and LPP) were reported to be involved in carcinogenesis. Rev3l, which encodes an error-prone DNA polymerase Pol ζ subunit, is critical for embryonic development and maintenance of genome stability [1819]. Loss of Rev3l enhances spontaneous mammary tumorigenesis, and genetic variation in the Rev3l gene is reported to affect the development and progression of breast cancer [2021]. The expression of the immunoglobulin superfamily 11 gene (IGSF11) was elevated in colorectal cancers and hepatocellular carcinomas as well as in intestinal-type gastric cancers, and IGSF11 is suggested to be a potential target for cancer immunotherapy [22]. MAML2, a member of the mastermind-like (MAML) family, functions as an essential coactivator of NOTCH in regulating a variety of malignancies [23]. There is increasing evidence that lipoma preferred partner (LPP) plays an important role in the development of a variety of human cancers [24]. A recent report revealed that LPP, which was normally operative in cells of mesenchymal origin, could be co-opted by breast cancer cells during an epithelial-mesenchymal transition to promote their migration and invasion [25]. Expression levels of the circRNA families transcribed from the above-mentioned four genes were all remarkably higher than those in their linear counterparts. The abundance of circRNA family from Igsf11 was 1,700-fold higher than that of Igsf11 mRNA expressed in the mammary gland. These data indicate that although the functions of these circRNAs remain unclear, they are unlikely to be the result of background "noise."

Table 3

The 10 most abundance co-expressed circular RNA families in rat mammary gland at day 1 and day 7

Genomic region	Day 1		Day 7		Gene
	circRNA*	linRNA†	circRNA*	linRNA†
NW_003812807: 1828679...1829413	2,055.8 (1)	87.9 (1,019)	1,562.6 (1)	71.0 (1,310)	Nip30
NW_003813123: 1588160...1589235	569.9 (2)	3.9 (9,366)	556.2 (2)	3.4 (9,162)	Rev3l
NW_003811247: 99932...100069	515.5 (3)	ND	0	ND	Intergenic
NW_003811537: 153458...165077	172.7 (20)	2,408.2 (31)	365.8 (3)	1,919.2 (32)	Rn45s
NW_003810754: 1053550...1055309	513.0 (4)	0.3 (14,430)	342.0 (4)	0.2 (14,802)	Igsf11
NW_003807897: 6882661...6883299	423.6 (5)	6.9 (7,668)	65.7 (37)	8.6 (6,511)	Phtf2
NW_003809175: 1605166...1605873	373.7 (6)	ND	26.1 (72)	ND	Intergenic
NW_003811373: 1201130...1201246	362.6 (7)	ND	40.0 (51)	ND	Intergenic
NW_003810582: 2592261...2592708	320.6 (8)	8.1 (7,137)	5.8 (234)	7.0 (7,209)	LOC100361417
NW_003809684: 1179751...1180471	306.6 (9)	14.8 (5,128)	190.8 (5)	36.0 (2,520)	Maml2
NW_003810792: 4017988...4018799	299.9 (10)	7.1 (7,540)	196.3 (6)	9.0 (6,400)	Lpp
NW_003812337: 587010...587645	230.4 (14)	ND	232.0 (7)	ND	Intergenic
NW_003809106: 939227...940296	151.0 (22)	ND	217.4 (8)	ND	Intergenic
NW_003806649: 2252840...2253928	20.2 (150)	11.5 (2,073)	205.0 (9)	19.8 (2,924)	Polr3g
NW_003806649: 2161597...2161746	272.8 (11)	ND	200.2 (10)	ND	Intergenic

circRNA=circular RNA; linRNA=mRNA of the circular RNA generating gene; ND=not determined.

*The number in bracket is the abundance ranking of a circRNA among the identified circular RNAs; †The number in bracket is the abundance ranking of a mRNA among the expressed mRNAs.

Recent reports have shown that, in mammals, circRNAs may function as miRNA sponges or compete against endogenous RNAs [456]. However, the enrichment of miRNA binding sites was not a global feature under selection in the sequence of the thousands of circRNAs profiled in the present study, which is consistent with the findings of a recent report [10]. However, we did identify a highly expressed circRNA that had several miRNA-binding sites. CircRNA_1093, with a short length of 108 bp, had 4 miR-342-3p binding sites. Notably, miR-342 was recently reported to regulate BRCA1 expression in breast cancer [26]. Therefore, we speculate that these circRNAs may be involved in the regulation of oncogenesis. In addition to their general regulatory functions in normal development, miR-200 and let-7 families are relatively well characterized in breast cancer initiation and self-renewal [2728]. A total of 3,100 circRNAs possess putative binding sites for members of the miR-200 family, of which 599 circRNAs have more than 10 binding sites, and a total of 2,539 circRNAs have putative binding sites for members of the let-7 family, of which 932 circRNAs have more than 10 binding sites. The high stability of circRNA in cells [11] suggests that they might be involved in the regulation of breast cancer via miRNAs. However, the potential genome-wide interplay between miRNAs and circRNAs needs further experimental and computational investigation.

Although the function of circRNAs remains largely unclear, an association between circRNAs and cancers has been reported. There is a global reduction in circRNA abundance in colorectal cancer cell lines and cancer compared with that in normal tissues, revealing a negative correlation between global circRNA abundance and proliferation and cancer [29]. Additionally, typical circRNAs were found to be show a significantly lower expression in gastric cancer tissues than in the paired adjacent nontumorous tissues [30]. These findings indicate that circRNA is a promising biomarker for cancer diagnosis. In the present study, a total of 6,824 and 4,523 circRNAs were identified by deep RNA sequencing from normal rat mammary glands at two lactation stages. These results provide the basis for comparison of breast cancer profiles and for selection of representative circRNA candidates for future functional characterization in breast development and breast cancer.

Supplementary Table 1

The polymerase chain reaction primers for circular RNA validation

circRNAs	Gene	Outward primer	Inward primer
circRNA_555	Tc2n	GGCACAATACCTTCCAAGTTCA	AAGAAGACACGCTTACTGAAGG
		CCTTCAGTAAGCGTGTCTTCTT	TGAACTTGGAAGGTATTGTGCC
circRNA_127	Arhgap10	GAACCTTATCGCTGCCACCAA	GCACAGTTGGATAAGATGGGAT
		GCACCGTTCATCATCAGTCTCA	ACTGCTTCCTTTCCTCTTCTGA
circRNA_3385	Hnrnpm	ATTGATGCCAACCACCTGAACA	CAGGATGAGAGGGCTTTACCAA
		GCTCCACCTCCTCCTTGCTT	ATGTCTGCTCGGACCACCA
circRNA_659	Cox4i1	TATGTTGATCGGCGTGACTACC	TTCAGTGTACCGCATCCAGTT
		CCACCACTGTCTTCCACTCATT	AGTCACGCCGATCAACATAAGA
circRNA_10008	Zmynd11	AGGCATATTGACTCCACTGACA	GAACCTTATCGCTGCCACCAA
		AACTGGCTGGCTCGCTTCT	GCACCGTTCATCATCAGTCTCA
circRNA_537	Tead1	GCTACCGCCATCCACAACAA	TCACATTCAGGTTCTTGCCAGA
		GGCTGTCCTGTCTGTATCATCC	TTGTCTTTCCCGTCCTCAGTTT
circRNA_388	Gnaq	CCTTCCTATCTGCCTACACAAC	GGCAAGAGTACCTTCATTAAGC
		CTCTCCTCCATTCGGTTCTCAT	AGCACATCAGTCACCTACAAC
circRNA_15703	Nars	ACAGGTTACGAAGGCAAGGAA	AAAGAGCCGTTTCCCACCATTT
		ACGCTACTCTCAGTGGACAAG	AACAACGCCATGTGCTCAGAA
circRNA_7070	Hdac4	AAGCACAAGCAACAGTTCCAG	GGGAATTTGGAGCGTGTTGGA
		GCCTCTACTTCTTCCTCCTCAC	CCACTGGCAGCGTCTTGTAA
circRNA_36159	Eef1a1	CTGATTGTTGCTGCTGGTGTTG	GCTTCTCTGACTACCCTCCACT
		TGTATGGTGGCTCGGTGGAA	GCAGCATCACCAGACTTCAAGA
circRNA_4815	Clmn	GCACCGACTCAGATTCTTCCTA	ACCGACTCAGATTCTTCCTACC
		GACCGCCACACCATACTTCC	GATGCTAACCAGCTTGACTTCC

circRNAs=circular RNAs.