A total of 60 healthy purebred Landrace sows (2–3 parity and without genetic relationships) were separated into two balanced cohorts of 30 animals each, and HS tests were conducted at a local thoroughbred farm during December 2016 and during August 2017 (WENS Shuitai Breeding Pig Farm, Guangdong, China). All sows were fed the same commercial formula diet and raised under the same management conditions. In the experimental stage, the ambient temperatures and relative humidity were measured at 14:00pm in everyday. In details, one cohort of 30 sows was selected in the winter months with a moderate average temperature, designated as the non-heat stress (NS) population; other cohort of 30 animals was selected in the summer months with a higher average temperature, designated as the HS population. Within each cohort, the number of piglets born alive was recorded, and litter birth weights of piglets were obtained within 24 hours after farrowing. Piglets were not offered creep feed, and sow milk was the only feed available to the piglets during lactation. On day 21, weaning survival and the weights of living piglets per litter were recorded and used to calculate average daily weight gain. Blood samples (10 ml) were collected at 10:00am from fasted sows using jugular venipuncture at weaning, and ELISA kits (Nanjing Jiancheng Bioengineering Institute, Nanjing, China) were used to determine serum LDH, IgG, SOD, HSP70, COR, and PRL concentrations. Sow milk samples (approximate 20 ml) were obtained on d 3, 15, and 20 of lactation from the last two pairs of sow nipples at 10:30–11:30am in each animal, and oxytocin was used to stimulate let-down. Three milk samples from each animal was pooled equally to evaluate the effect of environmental temperature on milk composition. In each environmental group, six animals that balanced for weaned backfat thickness and weight were chosen and humanely slaughtered at 21 day postpartum, and suckled mammary glands were split down the mid-line and tissues were excised from the center portion of four glands from the fourth and fifth pairs of nipples. Connective tissue and fat were removed. Mammary tissues were cut into small pieces and snap-frozen in liquid nitrogen prior to subsequent processing. In general, the collected mammary tissues contain primarily secretory epithelial cells, with a small amount of myoepithelial cells, endothelial cells, adipocytes, fibroblasts, and immune cells (Kensinger et al., 1986). All procedures were conducted under protocols approved by the Institutional Animal Care and Use Committee of South China Agricultural University, China.

Total RNA was extracted from mammary tissue and purified using Trizol reagent (Invitrogen, Carlsbad, CA), according to the manufacturer's instructions. Each RNA sample was treated with DNase I (Takara, Dalian, China) for 15 minutes at 37°C to remove residual genomic DNA. RNA quantity and purity were analyzed using a Bioanalyzer 2100 (Agilent, Palo Alto, CA), and RNA samples with Integrity Number (RIN) value ≥ 7.5 were used for further analysis. In each experimental condition, we randomly selected two RNA samples and pooled 5 μg of RNA from each sample. In total, six RNA pools were depleted of ribosomal RNA using an Ribo-Zero™ rRNA Removal Kit (Illumina, San Diego, USA), and the left poly-A^−/+ RNA fractions were then reverse-transcribed to create the final cDNA using a mRNA-Seq sample preparation kit (Illumina, San Diego, CA). Finally, we performed paired-end sequencing on an Illumina Hiseq 4000 (LC Bio, Hangzhou, China) to yield 2 × 150 nucleotide reads, following the manufacturer's recommended protocol.

Raw sequences quality was verified using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/), and the reads that contained adaptor contamination, low quality, and undetermined bases were removed by Cutadapt (Martin, 2011). Filtered data from each library was aligned to the Sscrofa11.1 reference genome downloaded from Ensembl genome website (ftp://ftp.ensembl.org/pub/release-94/fasta/sus_scrofa/dna/) with TopHat v2.1.1 (Kim et al., 2013), and transcript assembly and abundance estimation were performed using Cufflinks v2.2.1 (Trapnell et al., 2012). Each assembly was then merged using Cuffmerge to create a single transcriptome annotation with known porcine genes in gtf format (ftp://ftp.ensembl.org/pub/release-94/gtf/sus_scrofa) for subsequent protein-coding gene analysis. To predict circRNA candidates, five different algorithms including CIRCexplorer2 (Zhang X. et al., 2016), circRNA_Finder (Fu and Liu, 2014), CIRI (Gao et al., 2015), find_circ (Memczak et al., 2013), and MapSplice (Wang et al., 2010) were performed on each RNAseq library. Only circRNA candidates identified by all five approaches were considered for further evaluation. Following the above primary analysis, expression levels of all coding genes in each library were estimated from the TopHat alignments as fragments per kilobase of exon per million mapped reads (FPKM), and Cuffdiff, included in the Cufflinks package, was used to compare expression levels between NS and HS with a false discovery rate (FDR) value < 0.05. The abundance of circRNA candidates was calculated with back-spliced junction read count (Zhang et al., 2014). Then, the edgeR software package (Robinson et al., 2010) was used to examine the differential expression (DE) of circRNA candidates with P value < 0.05 and fold change ≥ 2. Finally, Biological Processes GO terms and KEGG pathway analysis of the DE genes were performed using DAVID gene functional classification (https://david.ncifcrf.gov/).

Putative interactions between the DE circRNAs and lactation-related coding genes that responded to HS in our paper were evaluated by competing to bind with shared miRNAs. Porcine mature miRNAs published in miRBase (http://www.mirbase.org/) were prepared for further analysis. In details, the construction of ceRNA networks included three steps: (1) the correlations between DE circRNAs and lactation-related genes were calculated by the Pearson test, and only nodes in positive circRNA-gene interactions were retained; (2) the circRNA-miRNA and mRNA-miRNA interactions were predicted by miRanda algorithm (Betel et al., 2010) with with energies of ≤ −20.0 kcal/mol and no mismatch in the seed region (positions 2–8 in the 5′ end); (3) potential circRNA-miRNA-gene interactions were established and visualized using Cytoscape V3.4 (http://cytoscape.org/).

Total RNA from the NS and HS samples were isolated with Trizol reagent (Invitrogen, Carlsbad, CA), and cDNA synthesis was conducted using PrimeScript RT reagent Kit (Takara, Dalian, China) with random hexamers. Quantitative PCR was used to analyze the expression changes of the chosen transcripts with SYBR Premix Ex Taq II (Takara). All primers are listed in Table S1, and final expression data were calculated using the 2^−ΔCt method using porcine GAPDH as a reference gene.

All sequencing raw datasets have been deposited into the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) database (https://www.ncbi.nlm.nih.gov/sra/) with the BioProject accession number PRJNA578241.

The statistical analysis was performed using by SPSS 17.0 Statistics Software (Chicago, IL, USA). The results of ELISA and RT-qPCR analysis between two groups were compared with independent t-test; the correlation analyses of DE circRNAs with lactation-related coding genes were tested by function cor (x, y, use = “p”), and illustrated with function labeledHeatmap (Matrix, xLabels, yLabels) in R package WGCNA (http://127.0.0.1:11153/library/stats/html/cor.html).

This study was performed during December 2016 and during August 2017. During the winter experimental stage, ambient temperatures and relative humidity averaged 14.3 ± 0.81°C and 65.0% ± 0.69%, while the corresponding values for the hot season were 32.7 ± 1.40°C and 76.1% ± 0.38%, respectively. The effects of ambient temperature on the serum stress-associated variables of lactating sows are presented in Figure S1. Blood heat shock protein 70 (HSP70), lactate dehydrogenase (LDH), and immunoglobulin G (IgG) levels were significantly higher in the HS cohort compared with the NS group (P < 0.05). In contrast, serum superoxide dismutase (SOD) and prolactin (PRL) concentrations were lower in the HS population (P < 0.05). The high temperature group had significantly lower serum cortisol (COR) concentration than the NS group (P < 0.05). In addition, there was no effect of thermal stress on the number of piglets born alive, litter birth weights, or piglets alive per litter at weaning (P > 0.05, Table 1). In contrast, litter weights at weaning were significantly higher in the NS cohort than the HS cohort (P < 0.05). And a reduced average daily weight gain of piglets (193.9 ± 2.19 vs. 218.0 ± 1.89 g, respectively, for the HS and NS cohorts) was associated with the environment with a high temperature. Thermal stress also altered the milk composition of lactating sows. In particular, sow milk had less butterfat, and lactose when sows lactated in the hot environment (P < 0.05), and milk in the HS group had a downward trend for total milk solids and milk protein levels (Figure 1). HS individuals tended to have much lower casein alpha s1 (CSNαs1) and s2 (CSNαs2) distributions (P < 0.01), while casein beta (CSNβ) and casein kappa (CSNκ) concentrations decreased from 343.59 ± 6.42 μg/ml and 7.47 ± 0.16 μg/ml in the HS group to 259.14 ± 7.96 μg/ml and 6.35 ± 0.12 μg/ml in the NS group (P < 0.01), respectively. Under HS, we found no overall differences in whey acidic protein (WAP) concentrations between the NS and HS groups (P = 0.067), although there was a slightly decreasing trend in the HS group (Table 2).

Six cDNA libraries were constructed from porcine mammary tissues exposed to thermal stress or to suitable temperatures. Each RNA-seq library generated approximately 94.52 ± 2.83 million raw reads of 110 nt in length, representing about 3.70 ± 0.11 fold coverage of the porcine genome. After quality control trimming, a total of 88.08 ± 2.65 million valid reads were obtained, accounting for 93.19% ± 0.23% of the raw reads in each library. We aligned all valid reads onto the porcine Sscrofa11.1 reference genome and found that over 80.55% ± 1.63% of the reads could be mapped successfully to the genome, including 76.13% ± 1.63% of the mapped reads with proper pair alignment in the six libraries (Table S2A). In addition, most valid reads were mapped to exons (77.51% ± 2.64%), 17.27% ± 2.02% were mapped to introns, and the rest to intergenic regions (5.22% ± 0.62%), indicating confidence in the quality of library construction and sequencing analysis.

Transcript assemblies from porcine mammary tissue with Cufflinks revealed a total of 133,145 isoforms across six samples, including approximately 37.33% identified candidates that completely matched Ensembl transcript regions (Table S2B). A comparison of known Ensembl transcripts revealed that 36,271 isoforms were expressed across all tissues; NS-specific units accounted for 82.08% of known Ensembl isoforms, while 81.37% of the known isoforms existed in the HS libraries (Table S2C). Raw Ensembl gene expression levels were quantified by the FPKM algorithm, and the 10 most prevalent functional isoforms accounted for 7.97% ± 0.62% of the total raw reads. In addition, seven gene candidates of PAEP, CSN1S1, CSN3, CSN2, JCHAIN, COX1, and NUPR1 were shared in the top 10 expressed genes in each sequencing library. These highly expressed isoforms are well-known as having key functions in the lactation process, consistent with the physiological roles of genes expected to be found in mammary gland tissues.

Several tools have been developed for identification of circRNAs based on back-spliced reads produced from high-throughput RNA sequencing datasets (Hansen et al., 2015). Due to the rearranged exon ordering, these back-spliced events usually cannot be mapped onto the reference genomes (Zhang et al., 2014). In the present study, we identified 17.07 ± 0.25 (19.75% ± 0.80% of the valid reads) and 17.07 ± 3.19 (19.16% ± 3.55%) million unmapped reads in the NS and HS libraries, respectively. Among them, we found a total of 948.00 ± 23.98 thousand back-spliced junction events (1.09% ± 0.04% of the valid reads) in the rRNA-depleted libraries of NS animals, as well as 892.49 ± 80.27 thousand candidates in the HS libraries. We then compared five different circRNA predicting algorithms and found a total of 31,031 unique circRNAs identified across six libraries. Of these, 19,642 were found by a single algorithm, accounting for 63.29% of all the circRNAs identified (Figure 2), indicating that the circRNA landscape differs quite dramatically depending on the algorithm used. In particular, find_circ and MapSplice exhibited the highest and lowest level of sensitivity; this is in part reflected in the total number of circRNAs predicted, where find_circ and MapSplice output the highest and lowest number of circRNA species (27,439 and 2,399, respectively) compared to CIRCexplorer2, circRNA_finder, and CIRI methods (7,865, 10,451, and 6,841 species, respectively; Table S3). To limit the level of false positive circRNAs, only circRNA candidates identified by all five approaches were considered for further evaluation. Of the 31,031 predicted circRNAs, only a modest overlap of 1,728 circRNAs (5.57%) was observed among the five prediction pipelines. These 1,728 circularization events were found to be produced from 1,157 hosting isoform loci, including 571 transcripts that generated more than one circRNA. For instance, we found that the porcine genes SEC24A and SLC5A10 had nine and eight predicted circRNAs, respectively, and there were seven circularization events predicted from the BAZ2B, PIAS1, and CCAR1 genes (Table S4A).

To identify dissimilarities between the tested individuals, principal component analysis (PCA) of the globally expressed transcript with the FPKM levels was performed (Figure 3). This analysis indicated that the differences in expression between the NS and HS groups were greater than the differences between pools from each particular group. Therefore, we employed the Cuffdiff algorithm to analyze differences in mammary gland gene expression between the NS and HS groups to identify candidate transcripts that are responsive to thermal stress. In our dataset there were over 40,000 unique Ensembl transcripts sequenced, most of which had a very small FPKM value in total across all libraries. In order to filter out false-positive results, we only kept confident transcripts that were expressed in at least three libraries, and finally 9,789 tags were identified in our study (Table S5A). Among these, we detected a total of 67 differentially expressed (DE) transcript events by a limited cut-off of FDR < 0.05, representing 63 protein-coding genes, with 42 transcripts increased and 25 transcripts decreased in the HS groups compared to the NS group (Table S5A). Analysis of gene ontology (GO) enrichment for DE genes, using identified porcine genes as background in the current experiment, revealed that these genes were significantly enriched in lactation-related functions or stress-inducible biological processes, including “lactation”, “defense response”, “inflammatory response”, “response to stimulus”, and “regulation of immune system process” (Table S5B). These 67 DE genes are significantly involved in only one KEGG signaling pathways, termed as toll-like receptor signaling pathway. Although the role of these genes needs to be validated experimentally, the GO and KEGG pathway analyses collectively illustrate some possible avenues to improve HS resistance of lactating sows. In addition, we also clustered differential porcine circRNA expression counts between NS and HS libraries, as determined by the CIRCexplorer2 algorithm, and normalized with trimmed mean of M-values (TMM) (Robinson et al., 2010). In total, only 50 out of 1,728 identified circRNAs were DE, including 29 up-regulated candidates and 21 down-regulated candidates in the NS samples compared with the HS samples (Table S4B).

To identify possible correlations between circRNA and mRNA expression, we first used the Pearson test and found a total of 1,423 significant interactions between DE circRNAs and DE genes in our study (Table S6). In Table S6 we identified 464 sponge modulators participating in 100 miRNA-mediated regulatory interactions, including 45 circRNAs and 36 unique mRNAs. In addition, the Pearson analysis also revealed 84 significant interactions between DE circRNAs and four lactation-related coding genes (CSN1S1, CSN1S2, CSN3, WAP) that were annotated by GO analyses, and these interactions included 42 positive significant interactions and 42 negative interactions (Table S6). We observed that the highest expressed circRNA, circCSN1S1_2, was significantly and positively associated with the expression of the CSN1S1, CSN1S2, CSN3, and WAP genes; interestingly, these gene products represent the main components of lactoprotein. Recent reports also showed that diverse RNA species can communicate with and co-regulate each other by competing to bind with shared miRNAs, acting as competing endogenous RNAs (ceRNAs) (Tay et al., 2014). We constructed circRNA-miRNA-mRNA networks by pairing the shared miRNA recognition sequences. In total, we generated a network that contained eight nodes and five connections formed between four circRNAs, three miRNAs, and CSN1S1 gene, including five circRNA-miRNA interactions and three mRNA-miRNA interactions (Table S7). Of these, the circCSN1S1_2-miR-204 (miR-670)-CSN1S1 ceRNA axis was established, corresponding well with the ceRNA hypothesis.

Based on ceRNA networks involved in mammopoiesis and lactogenesis under HS that we constructed, we identified a total of 10 interaction core genes (Figure S2). We validated expression of these core genes by RT-qPCR, including eight lactation-related coding genes, one heat-response gene, and circCSN1S1_2, which had the highest expression levels in our study. Divergent and convergent primers were designed for circRNA candidates according to the method described in previous study (Sun et al., 2017). All tested candidate genes showed consistent expression patterns between RT-qPCR and sequencing analysis, suggesting that our estimation of abundance was accurate. Briefly, the expression of coding genes for lactoprotein were significantly lower in HS sows (P < 0.05), except for CSN3, which decreased in HS sows, but the change was not significant. Usually, HSP family member proteins have important roles as molecular chaperones that help prevent apoptosis under various stress conditions, including HS (Sakatani et al., 2012). In agreement with the sequencing analysis, HSP90AA1 showed a strong induction in response to HS, and it had high expression levels in the HS group. In addition, expression of circCSN1S1_2 was significantly lower in HS sows, demonstrating the validity of our post-transcriptional ceRNA regulatory model.