High-throughput analysis of CircRNA in cows with naturally infected Staphylococcus aureus mammary gland

Abstract Circular RNAs (CircRNA) are a special type of non-coding RNA molecule with a closed ring structure and are not affected by RNA exonucases. It has stable expression, is not easy to degrade, and exists in most eukaryotes. However, circRNA regulation of cow mastitis has not been widely recognized. Mammary epithelial tissues were collected from healthy Holstein cows (HCN) and mastitis Holstein cows (HCU). RNA sequencing (RNA SEQ) was performed for the differentially expressed circRNAs, and analysis results showed that 19 differentially expressed circRNAs were identified in HCN and HCU, among which 6 circRNAs were up-regulated and 13 circRNAs were down-regulated. We randomly selected nine circRNAs for Q-PCR verification, and the results showed consistent expression. Three circRNAs: circRNA2860, circRNA5323 and circRNA4027 were confirmed to be significantly differentially expressed circRNAs in cow mastitis. Also, their host genes TRPS1, SLC12A2 and MYH11 might be directly or indirectly play a role in cow mastitis. Furthermore, RNA polymerase transcription factor binding and tight junction are most enriched in GO and KEGG pathways, respectively. In addition, the regulatory network of circRNA-miRNA has been inferred from a bioinformatics perspective, which may help to understand the underlying molecular mechanism of circRNAs involved in regulating mastitis in cows.


Introduction
Mastitis in dairy cattle and buffalo is a clinical condition that causes significant economic losses and is considered one of the enormous constraints to the dairy industry worldwide. 1 Staphylococcus aureus is recognized worldwide as one of the cattle's main contagious mastitis agents.It can express a set of antimicrobial resistance genes and virulence-associated genes that explain the wide range of outcomes of intramammary infections. 2 Staphylococcus aureus is a pathogen that is the causative agent of several human and veterinary infections and plays a critical role in cattle's clinical and subclinical mastitis.Staphylococcus aureus survival in cells is a significant cause of chronic persistent mastitis infection.However, it is unclear whether Staphylococcus aureus can escape autophagy in innate immune cells. 35][6] Furthermore, numerous types of toxins and enzymes in the milk produced by Staphylococcus aureus can lead to severe food-borne diseases. 7In addition, several essential pathways have been demonstrated to be related to the formation of cow mastitis, for example, TLR4/NF-jB pathway, 8 PI3K/Akt/mTOR signaling pathway, 9 and other signaling pathways.However, there is no systematic study on the molecular regulation of cow mastitis in mammary epithelium.
CircRNAs are a class of long, non-coding RNA molecules that form covalently closed continuous rings with a relatively stable framework and have a high tissue-specific expression in the eukaryotic transcriptome. 10Earlier, circRNAs were found and considered to have no biological function. 117][18][19] Many studies have demonstrated that circRNAs contribute to the generation of cancer, 20,21 regulate gene expression in many biological processes, and participate in the occurrence and development of various diseases. 22n the present study, we aimed to study the expression of circRNA in dairy cow mastitis.RNA-seq identification of differentially expressed circRNA in HCN and HCU.Functional analysis of differentially expressed circRNA was performed by gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment.In addition, the circRNA-miRNA network diagram was constructed to explore dairy mastitis further, but the molecular mechanism of mastitis has not been fully determined.This study provides a valuable reference for further understanding the relationship between circRNA and dairy cow mastitis and helps inhibit the development of dairy cow mastitis.

Ethics statement
All experiments in this study were approved and conducted according to the Animal Experimental Committee of Shenyang Agricultural University, Shenyang, China (201606005).

Sample preparation
We collected scapular mammary tissues from six Holstein cows (three HCN and three HCU).The cows we choose are 1 to 2 years old with 1 to 4 parities.The sample size was 2 cm 2 of female adults within the growth period, which had the same feed conditions, same body conditions, and no other diseases.We used procaine for local anesthesia to reduce animal pain.Mammary tissues were frozen in liquid nitrogen and stored at À80 C until RNA was isolated.

Total RNA isolation, library construction and sequencing
The total RNA of each group was extracted by Takarakit, the concentration and mass of RNA were determined by NanoDrop2000 spectrophotometer, and the integrity of RNA was determined by agar-gel electrophoresis.Total RNA was stored at À80 C for later use.An equal amount of RNA was taken from each sample for use.Ribosomal RNA was removed using the Epicenter Riboo-Zero rRNA removal kit (Epicenter, Madison,WI, USA) and rRNA was removed by ethanol precipitation.The NEBNext ultrafine RNA library kit prepared with Illumina (NEB, Ipswich,MA, USA) was used to generate two libraries from rRNA-deficient RNAs.The preparation of two libraries was classified on Illumina Hiseq 2000 platform.Generate a 150 bp end-to-end read.The quality of the original data was evaluated by FastQC.More than 10% of unknown bases are read, including adapters, and low-quality base reads are deleted by FastQC.

Identification of circRNAs and analysis of differentially expressed circRNAs
The cDNA libraries were performed the paired-end sequencing on an Illumina Hiseq 4000 (LC Bio, China) following the vendor's recommended protocol.Firstly, low-quality reads and adapters were removed by Cutadapt v1.10, quality controlled by FastQC v0.10.1, and then obtained the high-quality clean reads.TopHat v2.0.4 was utilized to map the clean reads to the reference genome from National Center for Biotechnology Information (NCBI) (https://www.ncbi.nlm.nih.gov/genome/?term=Capra+hircus). 23,24 Also, StringTie v1.3.0 was used to assemble and quantify expressed genes and transcripts (https://ccb.jhu.edu/software/stringtie/index.shtml). 25CIRCExplorer2 v2.2.6 software and the following criteria were used to identify candidate circRNAs: mismatch 2, back-spliced junction reads ! 1, and distances of two splice sites of less than 100 kb in the genome. 26Then, the back-spliced reads with at least two supporting reads were annotated as circRNAs.The differential expression of circRNAs between the two groups was assessed using the Ballgown package.A p-value < 0.05 and jlog2 (fc)j > 1 were set as the threshold for differential expression, 27,28 and R's phetmap package was used to draw the heatmap.

Gene ontology (GO) analysis and KEGG analysis of host genes
GO analysis (http://www.geneontology.org) was applied to differentially expressed circRNA-hosting genes.Similarly, pathway analysis uncovered the significant pathways related to differentially expressed circRNAs according to the annotation of the Kyoto Encyclopedia of Genes and Genomes (KEGG) (http:// www.kegg.jp/kegg). 29A threshold of p < 0.05 was used as a criterion for the determination of whether the enrichment analysis was significant. 30

Quantitative real-time PCR validation
We randomly detected 6 differentially expressed circrnas for Q-PCR.To demonstrate the digestibility of circrnas to RNase R, we treated total RNA with RNase R prior to cDNA synthesis.To verify the differentially expressed circrnas, total RNA was synthesized directly into cDNA using RT-PCR kits.According to the manufacturer's instructions, realtime PCR was performed using SYBR Green (TaKaRa Biotech, Dalian).The glyceraldehyde À3phosphate dehydrogenase (GAPDH) gene was used as internal control to make circRNAs expression level normal. 33hree independent experiments were conducted on HCN and HCU mammary tissue samples.Six pairs of primers were designed by primer 5 software (www.premierbiosoft.com)and are listed in Supplementary Table S1.Different analysis of the relative expression level circRNAs 2 -DDCt method qPCR data. 34Data were presented as mean ± standard deviation (n ¼ 3).SPSS 22.0, Chicago, IL, USA was used for statistical analysis of the two groups of data.P < 0.05 was considered statistically significant.In addition, circRNAs were verified by Q-PCR under the same experimental conditions to find the differential circRNAs in HCN and HCU.

Identification of circRNAs in HCN and HCU
A total of 13,679 circRNAs were identified from the RNA-seq data, including 705 circular intronic RNAs (ciRNAs) (Fig. 1a).After deleting the low-quality original read data, the total read data was obtained.The mapping rates of HCN and HCU were 94.93% and 94.76%, respectively.The HCN samples (45.87%) were compared with the HCU samples (50.6%), and the percentage of mapped sequence reads that could be aligned with the exon region was significantly lower (Fig. 1b).We analyzed the expression levels of HCN and HCU in circRNA, and the results showed that the expression of HCN was lower than that of HCN (Fig. 1c).

Identification of differentially expressed circRNAs in HCN and HCU
We identified 19 differentially expressed circRNAs in HCN and HCU.6 and 13 proteins were up-and down-regulated in HCN and HCU, respectively (Table 1).Volcano plots representing these circRNAs were prepared (Fig. 2a), with differentially expressed circRNA (!1.2-fold, P < 0.05) being located in the upper quadrant.Hierarchical clustering analysis was performed for the differentially expressed circRNAs to display circRNAs abundance differences between groups better (Fig. 2b).

Functional analysis of differentially expressed circRNAs
Different gene products coordinate with each other to perform biological functions, and pathway annotations of differentially expressed circRNA host genes provide a better understanding of the function of these genes.19 circRNAs were enriched in GO terms, and the top 25, top 15 and top 10 in biological processes, cellular components and molecular functions, respectively (Fig. 3a).Different gene products coordinate with each other to perform biological functions, and pathway annotations of protein genes provide a better understanding of the function of these genes.The RNA polymerase transcription factor binding is the most enriched in GO terms (Fig. 3b), and the tight junction is the most enriched in KEGG pathway (Fig. 3c).They are closely related to cow mammary tissue epithelial tissue and could be involved in the regulation of cow mastitis.

Validation of circRNAs by qRT-PCR
To investigate the expression of circRNA and determine that circRNA may play an important role in regulating cow mastitis, we used qPCR to verify the differential expression of certain circRNAs in HCN and HCU.Nine differentially expressed circRNAs were selected, and specific qPCR primers were designed in the circRNA junction region.RNA-SEQ results showed that circRNA2860 had the highest expression level in the up-regulated circRNA, while circRNA4027 had the highest expression level in the down-regulated circRNA.The qPCR experimental results of HCN and HCU are shown in Fig. 5.It proved that these  circRNAs existed and showed similar expression patterns in the epithelial tissue of a cow's breast, with the majority exhibiting a higher expression level in cow mastitis.The results of circRNA2860, circRNA5323 and circRNA4027 are significantly differentially expressed in RNA-seq and qPCR, suggesting that they might play a key role in cow mastitis.

Discussion
6][37][38][39] One of the most important pathogens causing mastitis is Staphylococcus aureus. 2 The immune response mechanism of Staphylococcus aureus infection is very complex and cannot be fully explained some mechanisms.Certain studies have shown that circRNA plays a vital role in Staphylococcus aureus infection 40 and breast cancer, [41][42][43] but seldom in cow mastitis.So we need to study it.In this work, high-throughput sequencing was used to revealed circRNA and gene expression in HCN and HCU.As a result,19 differentially expressed lncRNAs were identified by high-throughput sequencing and bioinformatics analysis.In addition, nine differentially expressed circRNAs were randomly selected for RT-QPCR validation to verify the defined circRNAs, and the results supported the reliability of RNA-SEQ data.
To further define the biological process that involves circRNAs, we performed GO annotation and KEGG pathway analysis.These processes could point to the roles of circRNAs in cow mastitis.In GO terms, differentially expressed circRNAs are enriched in Transcription from RNA polymerase II promoter and Adherens junction, including circRNA2584, circRNA2860, circRNA4513, circRNA8726 and their host gene TEAD1, TRPS1, CREBRF, AFDN, respectively.It has been revealed that differentially expressed circRNAs play a significant role in many vital biological pathways, including transcription from RNA polymerase II promoters.Sharifi showed that transcription from RNA polymerase II promoter, as a novel term, enriched by genes in his research of cow mastitis 44 ; Fu found host gene of circRNAs are involved in biological adhesion of biological processes in his study about bovine. 45The two terms are associated with circRNA2584, circRNA2860, circRNA4513 and circRNA8726 play crucial roles in the biological processes of cow mastitis.
In KEGG pathway enrichment, three significantly enriched pathways (p < 0.05) of differentially expressed circRNA were obtained, including Mucin type O-glycan biosynthesis, Glycine, serine and threonine metabolism, and Tight junction.The result shows the three pathways related to cow mastitis infected with Staphylococcus aureus.Tight junctions are intercellular junctions critical for building the epithelial barrier and maintaining epithelial polarity.Mastitis is a bacterial infection of mammary epithelial tissue, which destroys the protective barrier of Tight junctions. 46In our study, two host genes MYH11 of circRNA4027 and host gene AFDN of circRNA8726 were involved in Tight junction.MYH11 has been reported to be involved in breast cancer. 47urthermore, host gene TRPS1 is a highly sensitive and specific marker for breast carcinoma, suggesting it can be used as a great diagnostic tool 48 ; TEAD1, as a vector, interacts with genes in breast cancer. 49astitis is a chronic disease, an inflammatory state of the breast, and epidemiological studies have identified chronic inflammation as a factor in cancer risk.Hence, mastitis is a risk factor for breast cancer. 50herefore, host genes MYH11, TRPS1 and TEAD1 act in breast cancer and may also play a role in mastitis.Host genes are associated with their differentially expressed circRNA were identified in our data, so they are participating in cow mastitis to some extent.In addition, the enriched MAPK signaling pathway and Salivary secretion also can act in cow mastitis.It has been reported that multiple genes regulate dairy mastitis through MAPK signaling pathway, and MAPK signaling pathway has been involved in dairy mastitis as an important signaling pathway. 10,51,52The host gene ENSBTAG00000010619 of differentially expressed circRNA1119 was found on this pathway, so it might be playing a role in cow mastitis.Moreover, some reports have shown that Salivary SA (Sialic acids) was significantly higher and total protein was lower in breast cancer patients than in control. 53herefore, SA by Salivary secretions could be reduced to prevent breast cancer.Host gene SLC12A2 of circRNA5323 was involved in the pathway of Salivary secretion, and it might directly or indirectly affect cow mastitis.These results suggest that Tight junctions, MAPK signaling pathway and Salivary secretion are the key events in the regulation of cow mastitis, and circRNA2860, circRNA5323 and circRNA4027 with their host genes are playing vital roles in cow mastitis.
CircRNAs are known to function via regulating gene expression as miRNA sponges. 54Overexpression of miRNA target junction (miRNA sponge) results in loss of miRNA function and increased expression of endogenous target genes.Therefore, some circRNAs may inhibit or alleviate miRNA inhibition of translation. 16In our study, 121 miRNAs were matched with 9 differentially expressed circRNAs, in which miR-4530, miR-6844, miR-3198 and miR-34b-3p are shown involved in breast cancer, [55][56][57][58] so they might be potential miRNAs play roles in cow mastitis.Hence, the related circRNAs may act as miRNA sponges for regulating cow mastitis.

Conclusions
In this work, we performed an RNA-seq analysis that identified 13,678 circRNAs in HCN and HCU, of which 19 circRNAs were found to be differential expression.The most enriched pathways are functional analysis of differentially expressed circRNAs by GO and KEGG, RNA polymerase transcription factor binding and tight junction.The result of qRT-PCR confirmed that three circRNA, circRNA2860, circRNA5323 and circRNA4027 were significantly differentially expressed in HCN and HCU, they might act as key regulators of cow mastitis caused by Staphylococcus aureus.Meanwhile, their host genes TRPS1, SLC12A2 and MYH11 may be associated with cow mastitis.This study lays a foundation for further research on the molecular regulation of cow mastitis, and provides a reference for genetic traits desirable to breeding efforts in the control of mastitis in dairy cows.

Figure 1 .
Figure 1.Information on circRNAs from RNA-seq in Holstein cow normal (HCN) and I Holstein cow unhealthy (HCU) mammary tissue.(a) The types of circRNAs.(b) Distribution of exons, introns, and intergenic circRNAs.(c) The expression of circRNA.

Figure 2 .
Figure 2. Differentially expressed circRNAs in HCN and HCU.(a) Volcano map of differentially expressed circRNAs.Red dots indicate up-regulation and blue dots indicate down-regulation.(b) Cluster heatmap of differentially expressed circRNAs.The abscissa represents the sample and the log value of circRNA expression is regarded by the ordinate, which means that the heatmap is drawn from log10 of circRNA expression.The highly expressed circRNA is indicated by red; meanwhile, the lowly expressed circRNA is presented by blue.

Figure 3 .
Figure 3. Functional analysis of differentially expressed circRNAs.(a) Classification of GO terms.(b) GO analysis of differentially expressed circRNAs.The dot's color corresponds to different p-value ranges, and the size of the dot indicates the number of genes in the pathway.Rich factor denotes the number of differentially expressed circRNAs in the GO/the total number of circRNAs in the GO.(c) KEGG pathways of differentially expressed circRNAs.The dot's color corresponds to different p-value ranges, and the size of the dot indicates the number of genes in the pathway.Rich factor denotes the number of differentially expressed circRNAs in the KEGG/the total number of circRNAs in the KEGG.

Figure 4 .Figure 5 .
Figure 4.The circRNAs-miRNAs network.Network of 9 differentially expressed circRNAs.Red and green represent up and downregulation, and blue represent target miRNA.

Table 1 .
Up-regulated and down-regulated differentially expressed circRNAs in HCN and HCU.