Selective packaged circular RNAs in milk extracellular vesicles during Staphylococcus aureus infection may have potential against bacterial infection

ABSTRACT

Extracellular vesicles (EVs) provide a novel intercellular communication mechanism to transfer biologically important molecules to target cells. Although several pieces of evidence have shown that EVs have potential to respond to bacterial infections, our knowledge about the role of circular RNA (circRNA), an important cargo of EV, behind this process remains poor. In particular, the mechanism by which circRNAs are packaged into EVs remains elusive during bacterial infection. In the present study, EVs from bovine milk samples with or without Staphylococcus aureus (S. aureus) infection were isolated. The presence of circRNAs in milk-derived EVs (MEVs) was validated for the first time by PCR amplification with convergent and divergent primers and the RNase R resistance test. Through high-throughput sequencing, the expression profile of circRNAs in EVs was found to be changed during S. aureus infection. Moreover, we demonstrated that circRNAs were selectively packaged into EVs. Finally, bioinformatic analyses predicted the involvement of differentially expressed circRNAs in immune functions. In summary, our findings offer an insight into the packaging mechanism of EV circRNAs and underscore the potential by which host used the EV circRNAs in response to pathogenic bacterial infections.

KEYWORDS:

• Circular RNA
• extracellular vesicles
• Staphylococcus aureus
• expression profile
• packaging mechanism

Introduction

Infectious disease is the leading cause of mortality and morbidity worldwide, claiming about 15 million lives each year [1]. Among the pathogenic bacteria, Staphylococcus aureus (S. aureus) is a major cause of community and health care associated infections [2], seriously threatening human health and causing huge economic losses. For example, an estimated 119,247 S. aureus bloodstream infections occurred in 2017 with 19,832 associated deaths [3]. Additionally, S. aureus is a predominant pathogen causing bovine mastitis, which causes the greatest economic losses to the dairy industry [4]. Because of its pathogenicity, contagiousness, colonization of skin or mucosal epithelia, and poor cure rates associated with current therapies, S. aureus induced bovine mastitis is hard to be eliminated [5]. Furthermore, numerous types of toxins and enzymes in the milk produced by S. aureus can lead to severe food-borne diseases [6].

In order to protect a host from bacterial infections, vertebrates have evolved a number of mechanisms for the immune system to recognize and respond to pathogenic organisms. With the deepening of research, it has become increasingly clear that extracellular vesicles (EVs) can be used as a vehicle for intercellular communication to modulate immune functions during bacterial infection [7,8]. EVs, including exosomes with a diameter of ~100 nm, are surrounded by a lipid bilayer and generated inside multivesicular endosomes or multivesicular bodies and are secreted when these compartments fuse with the plasma membrane [9]. After being released by donor cells, EVs are present in perhaps all biological fluids including blood, saliva, urine, amniotic fluid, bronchoalveolar lavage fluid, and milk [10]. Although initially thought to be cellular debris, and thus under-appreciated, EVs are now increasingly recognized as important vehicles of intercellular communication and circulating biomarkers for disease diagnosis and prognosis [9]. Recent studies have shown that immune cells such as macrophages, dendritic cells, T cells, and natural killer cells use EVs to perform various functions in the immune response [11]. Membrane-associated proteins of EVs released from human macrophages were significantly changed by infection of Mycobacterium tuberculosis (M. tuberculosis) [12]. Moreover, infection of monocyte-derived macrophages with Mycobacterium bovis Bacillus Calmette-Guérin led to release of several specific EV miRNAs [13]. Bacteria-infected non-immune cells can also use EVs to induce an immune response. For example, human cervical epithelial cells released EVs containing elevated levels of a cellular inhibitor of apoptosis 2 (cIAP2) after Neisseria gonorrhoeae infection and suppressed apoptosis of the epithelial cells [14]. Our previous study also showed that milk-derived EVs, presumably released by mammary epithelial cells, contained differentially expressed miRNAs after S. aureus infection [15]. However, little is known about the roles of other EV cargos, such as circular RNAs (circRNAs), in bacterial infections.

CircRNAs are a type of single-stranded RNAs with a covalent 5ʹ to 3ʹ linkage loop structure formed by a process called back-splicing [16]. The proposed functions of circRNAs include sequestration of microRNAs (miRNAs) or proteins, modulation of transcription and interference with splicing, and even translation to produce polypeptides [17]. Thus, circRNAs may play considerable roles in bacterial infection. For example, circRNA_101128 was involved in peripheral blood mononuclear cells’ (PBMCs) response to M. tuberculosis infection [18], while circ_0005836 and circ_0009128 were significantly down-regulated in PBMCs of active pulmonary tuberculosis [19].

Although several studies have shown that circRNAs were enriched and stable in EVs derived from a number of bio-fluids, such as serum [20], urine [21], cerebrospinal fluid [22] and saliva [23], very little is known about EV circRNAs in milk. A few questions remain unsolved, such as whether circRNAs exist in MEVs. If so, how are circRNAs packaged into MEVs? Finally, what are the roles of MEV circRNAs in defending against bacterial infections?

In the present study, we aimed to address these questions. We showed the presence of circRNAs in MEVs for the first time. In addition, involvement of EV circRNAs in bacterial infection was investigated by comparing expression of circRNAs in EVs derived from milk samples with and without S. aureus infection and linking differentially expressed circRNAs to immune functions. Finally, we revealed that circRNAs were actively packaged into MEVs. These results suggest that hosts may use EVs as a tool to transport functional circRNAs to fight bacterial infection.

Methods and materials

Collection of milk samples

The milk samples were the same as the samples used in our previous research [15]. The details for these milk samples were described in the Supplementary Table S1 [15], including the Cow ID, the udder, the number of colonies in BPA after 24 h and 48 h culture and the test results of the S. aureus-specific thermonuclease gene (nuc). The milk samples from all four udders of 76 cows were aseptically collected as previously described [15]. The animal use protocol was approved by the Animal Use and Care Committee of Northwest A&F University (NWAFAC3751). The presence of S. aureus in milk samples was used for identifying S. aureus infections. S. aureus from milk samples were isolated, using a selective medium, Baird-Parker agar (Oxoid, Basingstoke, Hampshire, UK) and confirmed by the presence of nuc gene using PCR and bacterial 16S rRNA sequencing. Six samples from the S. aureus infected (the SA group) quarter and three from the non-infected (the control group) quarter were randomly selected and used for further studies.

Preparation and identification of MEVs

MEVs from these nine samples were prepared as previously described [15]. Due to the high complexity of milk components, ultracentrifugation was used to isolate mEVs. Briefly, milk samples (50 mL) were centrifuged at 5,000 × g (Sorvall ST 16 R, Thermo Scientific, USA) for 20 min at 4°C to remove milk fat and milk somatic cells. For removal of caseins and cell debris, skimmed milk samples were subjected to three successive centrifugations at 4°C for 1 h each at 12,000 × g, 35,000 × g (Avanti J-26 XP, Beckman Coulter, USA) and 70,000 × g (Optima™XPN, Beckman Coulter). The whey was collected and centrifuged at 135,000 × g at 4°C for 90 min (Optima™XPN, Beckman Coulter) to remove large particles. The supernatant was carefully collected and filtered through a 0.22 μm syringe-driven filter unit (Merck KGaA, Darmstadt, Germany). The percolate was collected and centrifuged at 150,000 × g for 90 min at 4°C (Optima™XPN, Beckman Coulter) to obtain the EV precipitate. The pellet was re-suspended in 1 mL of sterile PBS and filtered through a 0.22 μm syringe-driven filter unit (Merck KGaA). Finally, MEVs were stored in aliquots of 200 μL at −80°C until being used.

An aliquot of MEV preparation was used for protein estimation using the BCA kit (Thermo Scientific, Rockford, IL). MEV preparations, usually diluted by 10-fold, were compared in triplicate against serially diluted BSA as standard according to manufacturer’s instructions. Values were extrapolated from the standard curve, using a third-order polynomial equation, with r² > 0.98 for each assay.

To identify the MEVs, nanoparticle tracking analysis (NTA), TEM and Western blot were performed as previously described [15]. Briefly, nanoparticle sizes of MEVs were detected by NTA. Aliquot of 200 μL stored MEVs was diluted to 1 mL volume with sterile PBS stored on ice. The EVs solution was slowly injected into the sample cell of the Malvern Zetasizer Nano ZS90 (Malvern Panalytical Ltd., United Kingdom) system and measurements were taken according to manufacturer’s instructions. For TEM, MEVs were fixed in 3% (w/v) glutaraldehyde and 2% paraformaldehyde in a cacodylate buffer, pH 7.3. The fixed EVs were then applied to a continuous carbon grid and negatively stained with 2% uranyl acetate. The samples were examined with an HT7700 transmission electron microscope (HITACHI, Japan). Western blot assay was performed to assess the purity of exosomes. For protein extraction, MEVs were lysed in RIPA buffer and protease inhibitor cocktail (R0010, Solarbio, Beijing, China), and equal amounts of proteins (~30 μg) were subjected to Western blot assay. Total protein samples were separated by SDS-PAGE, followed by immune blot analyses with appropriate antibodies and detection with the Super ECL Plus substrate (P1050, Applygen, Beijing, China). The following antibodies were used: CD63 (1:1,000; System Biosciences, Inc., Palo Alto, CA, USA), CD81 (1:800, Santa Cruz, Santa Cruz, CA, USA), TSG101 (1:500, 14497-1-AP, Proteintech, Rosemont, IL, USA) and CD9 (1:1,000; System Biosciences, Inc.). Then, incubation with secondary antibodies was performed using either horseradish peroxidase (HRP)-conjugated goat anti-mouse (1:2,000; Dako Denmark A/S, Glostrup, Denmark) or anti-rabbit (1:20,000; System Biosciences, Inc.) antibodies. The immunoreactive bands were visualized using the Pierce ECL Western Blot Substrate (catalogue no. 32106, Thermo Scientific) according to the manufacturer’s instructions and exposed by Image Quant LAS4000 (GE Healthcare, Little Chalfont, UK). The protein band densities were assayed with the Image J2x 2.1.4.7 Analyser (Rawak Software Inc., Stuttgart, Germany).

Library construction and sequencing of circRNAs

To determine expression patterns of circRNAs in bovine MEVs, RNA in the EVs from milk in the control and SA groups were sequenced using Illumina HiSeq 2500 as described previously [24]. Briefly, total RNA in MEVs was rapidly extracted using TRIzol (Invitrogen, CA, USA). The concentration and integrity of RNA were determined using the NanoDrop 2000 (ThermoFisher Scientific, Waltham, MA, USA) and Bioanalyzer 2100 systems (Agilent Technologies, CA, USA), respectively. In order to construct cDNA libraries, the Epicentre Ribo-Zero ™ rRNA removal kit (Epicentre, Madison, WI, USA) was used to remove ribosomal RNAs from total RNA, and then the RNase R (3 U per μg of RNA) (Epicentre) was used to digest linear RNAs. Strand-specific RNA-seq libraries were prepared using the NEBNext Ultra Directional RNA Library Prep Kit for Illumina (New England BioLabs, Beverly, MA, USA) following the manufacturer’s instructions. Finally, the libraries were subjected to paired-end sequencing on an Illumina HiSeq sequencer (Illumina, San Diego, CA, USA). The datasets generated and analysed during the current study are available in the NCBI Sequence Read Archive (BioProject No.: PRJNA671561).

Bioinformatics analysis

Bioinformatic analyses were performed as previously described with some modifications [25]. Briefly, the FastQC software (v0.11.2) [26] was used to assess sequencing quality of raw data in fastq format. Low-quality data were filtered out using the NGSQC software (v2.3.2) [27]. The clean reads with high quality were then aligned to the bovine reference genome (Bos_taurus_UMD_3.1.1, http://www.ensembl.org/Bos_taurus/Info/Index) using the Tophat2 software (v2.0.13) [28]. Twenty-mers from both ends of the unmapped reads were extracted and aligned according to the reference genome to find unique anchor positions within the splice site. Anchor reads that aligned in the reversed orientation (head-to-tail) indicated circRNA splicing and were then subjected to find_circ [29] and CIRI [30] to identify circRNAs. The identified circRNAs were subjected to statistical analyses for type, chromosome distribution and length distribution.

Analysis of differentially expressed circRNAs

To quantify overall expression levels of circRNAs in MEVs after S. aureus infection, the total number of reads spanning back-spliced junctions was used as an absolute measure of circRNA abundance. To estimate the relative expression of a circRNA, SRPBM (Spliced Reads per Billion Mapped Reads) was calculated as (Spliced reads/Total mapped reads) × 10⁹ for each back-spliced junction. The DESeq2 package of R software was used to identify differentially expressed circRNAs between the two groups [31]. CircRNAs with total counts more than 10 and a fold change greater than 2 and a p-value <0.05 in the comparison between the two groups were identified as differentially expressed circRNAs.

The Gene Ontology (GO) analysis was performed using the online enrichment tool from the GO database powered by the PANTHER classification system (http://pantherdb.org/) [32]. Lists of genes corresponding to linear RNA or circRNA parental genes were provided as input for the PANTHER significant enrichment detection tool. The option of Bonferroni correction for multiple hypothesis testing was used for all samples. Molecular Function, Biological Process, and Cellular Component were analysed. Pathway-based analyses further elucidated the biological function of genes. Kyoto Encyclopaedia of Genes and Genomes (KEGG) enrichment analyses of the differentially expressed circRNAs was performed with the DAVID software [33]. The cut-off for detecting significant enrichment was set to p < 0.05. The interactions between circRNAs and their target miRNAs were predicted using the RegRNA 2.0 (http://regrna2.mbc.nctu.edu.tw/) by the miRanda algorithm [34,35].

Isolation of genomic DNA and total RNA from milk somatic cells

In order to validate the existence of expressed circRNAs, genomic DNA (gDNA) and total RNA were isolated from the somatic cells of nine milk samples. Milk somatic cells were isolated as described earlier [36]. The gDNA of milk somatic cells was isolated using the Easy-DNA gDNA Purification Kit (Invitrogen), while total RNA was extracted using TRIzol (Invitrogen) according to the manufacturer’s protocol. The total RNA from milk somatic cells was divided into three parts. One part of the total RNA was reverse transcribed to cDNA using an iScript RT SuperMix kit (Bio-Rad, Hercules, CA, USA). The cDNA and gDNA were used for PCR amplification with both divergent and convergent primers. The other two parts of the total RNA were stored at −80°C until the RNase R treatment.

PCR amplification using divergent and convergent primers

Using cDNA and gDNA as templates, PCR amplification was performed with the divergent and convergent primers (Supplementary Table 1). Divergent primers were used to flank the circRNA-specific junction (circular isoforms), while the convergent primers were used for the linear splice junction of downstream exons (linear isoforms). These primers were designed by an NCBI primer design tool and synthesized from Sangon Biotech (Shanghai, China), and 35 cycles of PCR were performed under standard conditions. The PCR products were then subjected to electrophoresis in 2% ethidium bromide-stained agarose gel. To confirm the PCR results, the PCR products were purified using an EasyPure Quick Gel Extraction Kit (Transgen Biotech, Beijing, China) and Sanger sequenced by Sangon Biotech (Shanghai, China).

RNase R treatment and real-time PCR analysis

The RNase R treatment of total RNA from milk somatic cells was performed as previously reported [37]. In brief, two parts of the total RNA samples were thawed quickly in a 37°C water bath. One part was treated for 15 min at 37°C with 3 U/μg of RNase R (Epicentre), while the other was added to an equal amount of RNase-free water under the same conditions as the control. Both RNase R-treated and untreated samples were subsequently purified by phenol-chloroform extraction, and immediately reverse transcribed to cDNA using an iScript RT SuperMix kit (Bio-Rad).

Real-time PCR was carried out using the TranStart Tip Green qPCR SuperMix (Transgen Biotech) on a Bio-Rad CFX96 Touch Real-Time PCR Detection System by using the cDNA as a template with the specific primers. The relative levels of circRNAs were determined by 2^−ΔΔCt and normalized to the housekeeping gene GAPDH.

RNase A treatment

The RNase A was used to treat MEVs according to previous research [29]. Briefly, milk EVs were incubated with RNase A (final concentration 0.5 µg/µL; RNase A/T1 mix, EN0551, ThermoFisher Scientific, Burlington, ON, CA) for 20 min at 37°C. Then, the extracellular vesicle RNA was extracted with TRIzol (Invitrogen, Carlsbad, CA, USA). Control samples were subjected to the same treatments, but by replacing RNase A with PBS. The RNA concentrations were measured by NanoDrop 2000 (ThermoFisher Scientific).

Packing of circRNAs into EVs

In order to reveal the packaging mechanism of EV circRNAs, the ratio of the relative levels of circular/linear isoforms in EVs from milk somatic cells was calculated. The relative expression of circular isoforms and linear isoforms was normalized to the housekeeping gene GAPDH. For RT-qPCR, eight biological replicates were used to calculate standard deviations (SDs).

Results

Isolation and biological characterization of milk-derived EVs

MEVs from bovine milk samples were isolated by differential centrifugation. To characterize the isolated MEVs biophysically, transmission electron microscopy (TEM) and nanoparticle tracking analysis (NTA) were performed. Most MEVs had diameters ranging from 50 nm to 200 nm by TEM (Fig. 1A). The mean diameters of bovine MEVs isolated from the control group and the SA group samples were 155.1 nm and 165.7 nm as measured by NTA (Fig. 1B). The average yields of the MEVs from the control group and the SA group were 1.40 × 10⁹ ± 5.14 × 10⁷ and 1.51 × 10⁹ ± 8.70 × 10⁷ particles/mL, respectively (Fig. 1B). Western blot analysis indicated the presence of surface marker proteins CD63, CD81, TSG101 and CD9 in bovine MEVs isolated from both the control group and the SA group samples (Fig. 1C). The concentrations of total protein and total RNA in MEVs from both the control group and the SA group samples were measured. Additionally, the protein to particle ratio and the RNA to particle ratio were calculated. Although the S. aureus infection decreased the protein (from 169.31 ± 21.45 to 114.90 ± 17.89 μg/mL) and RNA (from 32.53 ± 0.97 to 30.58 ± 1.42 ng/mL) concentrations of bovine MEVs, neither the protein to particle ratio (120.93 ± 15.32 μg/10⁹ particles for the control group and 114.90 ± 17.89 μg/10⁹ particles for the SA group) nor the RNA to particle ratio (23.23 ± 0.69 ng/10⁹ particles for the control group and 20.25 ± 0.94 ng/10⁹ particles for the SA group) changed significantly (Fig. 1D).

Figure 1. Isolation and characterization of MEVs. (A) Electron microscopy images of EVs isolated from bovine milk. The field image was obtained at 100,000× magnification. (B) Comparison of particle size distribution of MEVs isolated from the control group (green, n = 3) and the SA group (red, n = 6) samples. The dots represent the mean value, and the shading represents the SD. (C) Western blot analysis of surface marker proteins on isolated MEVs from the control group and the SA group samples. (D) RNA to particle ratio and protein to particle ratio for MEVs isolated from the control group and the SA group samples. The standard errors of the means (SEMs) are indicated by error bars

Characteristics of bovine circRNAs in MEVs

In order to explore the influence of S. aureus infection on the characteristics of bovine circRNAs in MEVs, the raw data were filtered by removing adapter and low-quality sequences. A total of 850 million clean reads were obtained in the control and SA groups (Supplementary Table 2). Most of the clean data were aligned to the bovine reference genome, ranging from 87.93% to 91.40% (Supplementary Table 2). The average numbers of circRNAs in the MEVs from the SA group were significantly lower than those from the control group (p < 0.05) (Fig. 2A). Further analyses revealed a total of 39,276 circRNAs. Among them, 15,011 circRNAs were shared in both the control and SA groups. However, the number of circRNAs only expressed in the SA group was much higher than that in the control group (Fig. 2B). These results indicate that S. aureus infection decreased the total amount of circRNAs in MEVs but increased the variety of MEV circRNAs.

Figure 2. Characteristics of the identified circRNAs in MEVs by RNA sequencing. (A) The counts of circRNAs in MEVs from the control (n = 3) and SA (n = 6) groups. (B) The Venn diagram demonstrating the intersection of circRNAs in the control and the SA groups. (C) Violin plots showing the length distribution of circRNAs in MEVs. The vertical axis represents the length of circRNAs in MEVs. The white dot and the box within each violin stand for the median and the interquartile range. (D) Distribution of the circRNAs in MEVs on the bovine chromosomes (the outermost circle, one unit of the scale stands for one million base-pairs). The middle circle (blue lines) and the innermost circle (orange lines) represent the circRNAs in the control group and the SA group, respectively. The height of the column is proportional to the expression level and the position of the column corresponds to the location of circRNAs on chromosomes. (E) Genomic origin of the circRNAs in MEVs. * p < 0.05; chr: chromosome; CDS: coding sequence; UTR: untranslated regions; TSS: transcription start sites; TES: transcription end sites

The length of bovine MEV circRNAs ranged from under 100 nt to over 10,000 nt, with most being about 1,000 nt (Fig. 2C). Furthermore, these circRNAs originated from every bovine chromosome and the most abundant MEV circRNAs (2,258 different types) were enriched on bovine chr1 (Fig. 2D). Among these MEV circRNAs, 92.82% of them were identified in exons, with 73.72%, 5.23% and 13.87% on CDS exons, 5ʹ-UTR exons and 3ʹ-UTR exons, respectively (Fig. 2E).

The relationship between differentially expressed circRNAs in MEVs and immune functions

In order to investigate how bacterial infection affected circRNAs in MEVs, differentially expressed circRNAs in EVs from milk samples without or with S. aureus infection were profiled. The results showed that 290 EV circRNAs were differentially expressed after S. aureus infection, based on fold change >2 and p-value <0.05 (Supplementary Table 3). Among them, 121 were up-regulated and 169 were down-regulated by S. aureus infection (Fig. 3A). To understand the biological function of these differentially expressed circRNAs, the KEGG (Kyoto Encyclopaedia of Genes and Genomes) pathway analysis and GO (Gene Ontology) annotation were performed. Three (‘Epstein-Barr virus infection’, ‘Salivary secretion’ and ‘Ubiquitin mediated proteolysis’) and five (‘Antigen processing and presentation’, ‘Phagosome’, ‘Protein processing in endoplasmic reticulum’, ‘Endocytosis’ and ‘Regulation of actin cytoskeleton’) KEGG pathways (p < 0.05) were identified among the up-regulated and down-regulated circRNAs, respectively (Fig. 3B). Furthermore, 80 enriched GO processes (p < 0.05) were revealed. Among them, 38, 23, and 19 belonged to categories of biological process, cell components, and molecular function, respectively (Supplementary Table 4). The top five enriched terms of each category were ranked by -log₁₀(p-value), and the number of genes in these GO terms was also counted (Fig. 3C). The top three enriched GO terms of biological process were ‘the response to dehydroepiandrosterone’, ‘the response to 11-deoxycorticosterone’, and ‘the microtubule nucleation’. The ‘Golgi lumen’ and ‘antigen binding’ were the top-enriched terms in categories of cell components and molecular function, respectively.

Figure 3. Analysis of immune function of differentially expressed circRNAs in EVs derived from bovine milk infected by S. aureus. (A) Heat map of differentially expressed circRNAs. Colour from blue (low) to red (high) represents expression levels of circRNAs. Each row in the map represents a sample, and each column represents a gene. (B) KEGG pathway analysis of up-regulated and down-regulated circRNAs. The sizes of the dots represent the counts of genes. The gene ratio indicates the ratio between the number of target genes associated with a KEGG term and the total number of genes in the KEGG term. (C) GO enrichment of the genes that produce differently expressed circRNAs. GO annotations of the linear counterparts of differently expressed circRNAs covering the domains of biological processes, cellular components, and molecular functions. The five most significantly enriched terms are shown for each category of circRNAs, ranked according to -log₁₀(p-value). The grey histogram represents number of circRNAs in certain GO terms. (D) Schematic diagram of the known relationship between the significantly upregulated circRNAs in MEVs and their targeted miRNAs as well as immune-related genes. The detailed information for this diagram is presented in Supplementary Table 6. KEGG: Kyoto Encyclopaedia of Genes and Genomes; GO: Gene Ontology

Figure 3. (Continued)

One of the functions of circRNAs is sequestration of microRNAs (miRNAs) [24]. miRNAs are a family of small non-coding RNAs that are important post-transcriptional regulators of gene expression [38]. We hypothesized that EV circRNAs (especially significantly up-regulated circRNAs) had the possibility to affect immune-related genes through interacting with miRNAs. Therefore, interactions between the eight most significantly upregulated EV circRNAs (circCSN1S2-Up, circMEMO1, circXDH, circCYFIP1, circRNF4, circTAF15, circCOPA, and circDENND1A) and their target miRNAs were predicted. A total of 94 miRNAs were predicted to be associated with these eight EV circRNAs (Supplementary Table 5). Thirty-five of these miRNAs were reported to have immune functions targeting at least 136 immune-related genes (Supplementary Table 6) in the NCBI PubMed database (https://www.ncbi.nlm.nih.gov/pubmed/). The reported relationship among eight EV circRNAs, their target miRNAs and the target genes for the miRNAs is presented in Fig. 3D.

Validation of circRNAs in milk somatic cells

Three groups of circRNAs (22 kinds) were selected from the sequencing data for validation of their existence in milk somatic cells: eight significantly up-regulated circRNAs (circCSN1S2-Up, circMEMO1, circXDH, circCYFIP1, circRNF4, circTAF15, circCOPA, and circDENND1A), six significantly down-regulated circRNAs (circCSN1S2-Down, circCSN1S1-1, circCSN1S1-2, circCSN1S1-3, circPICALM, and circCSN2) and eight circRNAs (circGNPTAB, circSLC39A8, circTNFRSF21, circARID1A, circSACM1L, circTXNDC11, circMAP7, and circHNRNPLL) with the highest expression levels but without significant changes between the two groups. The significantly up-regulated circRNAs were chosen with threshold log₂(Fold Change) >3 and p-value <0.05, while significantly down-regulated circRNAs were selected based on log₂(Fold Change) < −4 and p-value <0.05 (Fig. 4A). These 22 circRNAs ranged in length from 234 nt to 1161 nt, distributed on 11 bovine chromosomes (chr 2, 3, 5, 6, 9, 11, 19, 22, 23, 25, and 29), and cyclizing within 1 to 11 exons (Fig. 4B and Supplementary Table 7). It is worth noting that although both circCSN1S2-Up and circCSN1S2-Down originated from gene CSN1S2, circCSN1S2-Up was cyclized from Exons 14 to Exons 17, while circCSN1S2-Down was from Exons 6 to Exons 11. More importantly, compared to the control group, circCSN1S2-Up was significantly up-regulated in MEVs, while circCSN1S2-Down was significantly down-regulated.

Figure 4. Validation of circRNAs. (A) Volcano plots comparing the expression of EV circRNAs in the SA group with those in the control group. The horizontal axis represents the normalized circRNA expression difference (the SA group/the control group), while the vertical axis is the normalized p-value. The horizontal line represents p-value = 0.05. Red triangles and green-inverted triangles represent circRNAs, which were the most significantly up-regulated (log₂(Fold Change > 3), p-value < 0.05) and down-regulated (log₂(Fold Changes < −4), p-value < 0.05), respectively. Gray dots represented circRNAs that were filtered out by Fold Change>2 and/or p-value < 0.05. (B) Schematic illustration of 22 selected circRNAs genomic region. Arc lines indicate exonic regions that generate circRNAs. The colours of the arc lines represent different types of circRNAs. Red: the significantly up-regulated circRNAs; Green: the significantly down-regulated circRNAs; Orange: the circRNAs without significant changes. (C) PCR products of selected candidate circRNAs were analysed by gel electrophoresis with cDNA and gDNA from milk somatic cells as templates. Divergent primers amplify circRNAs in cDNA but not in gDNA. GAPDH, linear control, size marker in base pairs. ◄ ►stands for divergent primers, ► ◄ while stands for convergent primers. (D) The RNase R digestion resistance test of circRNAs by qPCR. Total RNA samples of milk somatic cells were digested with or without RNase R before reverse transcription (RT) reactions (right panel). The -ΔCt values were calculated for each of 9 samples (left panel) with at least four replicates. GAPDH was used as a linear control. The error bars represent ±SD. **, p-value < 0.01

Figure 4. (Continued)

In order to validate the presence of the 22 circRNAs, PCR amplifications with convergent primers and divergent primers were performed. The convergent primers amplified products from both the cDNA and gDNA, whereas divergent primers only amplified circRNAs from cDNA. As shown in Fig. 4C, all 22 circRNAs were detected in the samples. These results were then confirmed by Sanger sequencing, which showed the head-to-tail circularization splicing of these circRNAs (Fig S1).

The existence of the selected 22 candidate circRNAs was further confirmed by the RNase R treatment and real-time PCR analysis. The RNase R is an exonucleolytic enzyme that preferentially digests linear mRNA from the 3ʹ tails but not circular RNA. The qPCR results indicated that all the 22 circRNAs showed varying degrees of resistance to the RNase R digestion in comparison with the linear housekeeping gene GAPDH, which was significantly down-regulated to 13.68% (2.87 cycles) (Fig. 4D).

Selective package of circRNAs into EVs during S. aureus infection

The packaging mechanism of EV circRNAs was investigated. In order to eliminate the possibility that circRNAs might be attached to the surface of MEVs, MEVs were digested with RNase A and the RNA concentration was measured. The result showed that the RNA concentrations of MEVs did not change significantly by RNase A treatment (Fig S2). Then, expression of 22 circRNAs in the MEVs was detected by RT-qPCR with the circRNA-specific primers (Supplementary Table 1), and the results were consistent with the sequencing data generally (Fig S3). Compared with the control group (milk samples without S. aureus), expression levels of the circRNAs in milk somatic cells infected with S. aureus, regardless of the groups (the up-regulated circRNAs, the down-regulated circRNAs, or the circRNAs without significant changes), were significantly down-regulated, except for circGNPTAB and circHNRNPLL (Fig. 5A). The last two circRNAs (circGNPTAB and circHNRNPLL) did not show any difference between the control milk somatic cells and the milk somatic cells infected with S. aureus. These results indicated that the S. aureus infection could down-regulate the expression levels of most circRNAs in milk cells. More importantly, compared with the control group, the relative expressions of candidate circRNAs in milk somatic cells and MEVs of the SA group were clustered by the circRNA groups (the significantly up-regulated circRNAs, the significantly down-regulated circRNAs, and the circRNAs without significant changes) (Fig. 5B). There was no linear relationship of expression levels between milk somatic cells and MEVs of the SA group (Fig. 5C) and this refutes existence of a clearance mechanism by which cells secreted internal circRNAs into EVs indiscriminately.

Figure 5. The selective package of circRNAs in MEVs were regulated by S. aureus infection. (A) The relative expressions of the significantly up-regulated circRNAs (top), the significantly down-regulated circRNAs (middle) and the circRNAs without significant changes (bottom) in milk somatic cells (blue squares) and MEVs (orange dots). The expression level of circRNAs in the absence of S. aureus was set to zero. The red up and the green down arrows represent significantly up- and down-regulated circRNAs, respectively. * and ** were calculated by the relative expression of circRNAs in the SA group vs. circRNAs in the control group. # and ## were calculated by the relative expression of circRNAs in the MEVs vs. circRNAs in the milk cells. * and #, p-value < 0.05; ** and ##, p-value < 0.01. (B, C) The relative levels of candidate circular (B) and linear (C) isoforms in the presence of S. aureus infection compared with the absence of S. aureus infection. The relative expression in bovine MEVs (Y-axis), against those in milk somatic cells (X-axis). Ellipses with different colours represent clusters of candidates. (D, E) CircRNAs preferentially released into EVs under the normal and S. aureus infection, respectively. The relative circ/lin levels are plotted for MEVs (Y-axis), against those in milk somatic cells (X-axis). The upper and lower dashed lines respectively represent the relative expression up-regulated or down-regulated four-folds. CircRNAs preferentially packaged into EVs or retained in milk somatic cells are plotted in red or blue, respectively

Figure 5. (Continued)

Unlike the canonical splicing of linear RNAs, the biogenesis of circRNAs occurs via back-splicing [39]. Therefore, the packaging of circRNAs into EVs could be related to the relative levels of circular/linear isoforms (circ/lin). Thus, we calculated the relative levels of circ/lin in milk somatic cells and MEVs. In the absence of S. aureus infection, most of the 22 circRNAs (circTAF15, circMAP7, circTNFRSF21, circCSN1S1-3, circHNRNPLL, circMEMO1, circCSN2, circPICALM, circGNPTAB, circCOPA, circCSN1S1-2, circXDH, circARID1A, circCSN1S1-1, and circCSN1S2-Up) were preferentially retained in milk somatic cells and only one (circCSN1S2-Down) was preferentially packaged into EVs (Fig. 5D). On the other hand, during the S. aureus infection, five circRNAs (circDENND1A, circTAF15, circRNF4, circSACM1L, and circCYFIP1) were preferentially packaged into EVs (Fig. 5E). Meanwhile, 12 circRNAs (circCSN1S1-3, circHNRNPLL, circCSN1S2-Up, circXDH, circMAP7, circGNPTAB, circTNFRSF21, circARID1A, circPICALM, circCSN1S1-2, circCSN2, and circCSN1S1-1) were preferentially retained in milk somatic cells during S. aureus infection (Fig. 5E). Taken together, these data indicate that circRNAs were selectively packaged into MEVs during S. aureus infection.

Discussion

Bacterial infection remains as one of the major threats to human and animal health. It has been suggested that EVs could be used as potential tools to fight against bacterial infections [40]. This hypothesis has been supported by results from an increasing number of studies, which showed a close relationship between circRNAs and diseases caused by bacterial infections, such as active pulmonary tuberculosis, community‐acquired pneumonia and bacterial meningitis [41–43]. However, solid evidence is still lacking. In this study, we confirmed the presence of circRNAs in MEVs for the first time and revealed that the profiles of circRNAs in EVs were significantly altered by S. aureus infection. Such changes were induced by an active sorting process into milk somatic cells. These up-regulated circRNAs were linked to immune-related pathways. Consequently, EV circRNAs have the potential to play regulatory roles in response to pathogenic bacterial infection.

In the present study, MEVs were chosen as an appropriate model to study the involvement of circRNAs in bacterial infections as bovine milk is easily obtained in large quantities and large amounts of EVs in milk can be isolated and purified [44]. In addition, it has been reported that miRNAs of MEVs have the potential to respond to bacterial infection [15]. Finally, milk could represent a unique means of cell-to-cell communication between individuals within a species or across species.

Our research demonstrated that circRNAs were abundant in MEVs for the first time. Milk is a complex and dynamic bio-fluid comprised of hundreds to thousands of distinct bioactive structures, among which EVs are one of the most important substances. MEVs have the potential to be internalized by human intestinal crypt-like cells, human colon carcinoma Caco-2 cells, and rat small intestinal IEC-6 cells [45,46]. In addition, meaningful amounts of cargos in MEVs, such as miRNAs, can be detected in human blood after milk consumption [47]. Therefore, MEVs have the potential to act as the ‘bridge’ through which circRNAs can play regulatory roles between a mother and her offspring, or even across species (Fig. 3D).

The results from the current study suggest the involvement of EV circRNAs in response to bacterial infections. The circRNAs have been reported to be involved in regulating transcription in the nucleus [48], functioning as efficient microRNA sponges [49], competing with pre-mRNA splicing [50], and serving in circRNA–protein interactions [51]. In the current study, 290 circRNAs in MEVs were dysregulated significantly with S. aureus infection. Further, the KEGG analysis showed that the differentially expressed circRNAs in MEVs were enriched in pathways related to immunity (Epstein-Barr virus infection and antigen processing) and secretion (salivary secretion, protein processing in endoplasmic reticulum, endocytosis, and regulation of actin cytoskeleton). More importantly, the eight most significantly upregulated circRNAs in MEVs were predicted to affect immune-related genes through miRNAs. These results suggest the involvement of circRNAs in MEVs in fighting bacterial infections.

Although the presence of circRNAs in EVs has long been recognized, the sorting mechanism by which cells package circRNAs into EVs remains unclear. After being transferred from nuclei into cytoplasm, circRNAs can be packaged into EVs along with an abundant cargo of other nucleic acids, lipids, and proteins [52]. There are two main hypotheses about the secretion of circRNAs into EVs: clearance mechanism and selective release mechanism [53,54]. In the current study, we showed that the relative levels of circular/linear isoforms in milk somatic cells were much different from those in MEVs (Fig. 5B). More importantly, S. aureus infection selectively released certain circRNAs into EVs in comparison with their corresponding linear isoforms. These results provide additional support for the selective release mechanism.

Author contributions

XZ contributed to the experiment design, manuscript draft and data analysis. SM and CT contributed to the experiment implementation, manuscript draft and data analysis. SM and MN participated in designing the experiment. SM, MN, ZH and ML performed the experiment and analyzed the data. XZ and SM wrote the manuscript. All authors read and approved the final manuscript.

Acknowledgments

We would like to thank Genergy Bio Co., Ltd. (Shanghai, China) for assistance with bioinformatics analysis. This work was supported by the National Key Research and Development Program of China [2016YFD0500507].

Disclosure statement

No potential conflict of interest was reported by the authors.

Supplementary material

Supplemental data for this article can be accessed here.

Additional information

Funding

This work was supported by the National Key Research and Development Program of China [2016YFD0500507].