The previously frozen bovine mammary epithelial cell lines MAC-T (20) in liquid nitrogen were thawed and cultured in three 6-well plates supplemented with Dubecco's modified Eagle medium/ F12 containing 10% fetal bovine serum (System Biosciences, Mountain View, CA). When bMEC monolayers reached 80% confluence, one 6-well plate cells (0 h) were collected as the control group, and the remaining two 6-well plates were treated with 50 ng/μL LPS to stimulate an inflammatory response (21). The respective 6-well plates were collected at 6 and 12 hpi as the experimental groups.

Total RNA was extracted from bMECs induced with LPS at 0, 6, and 12 h using TRizol reagent (Takara Biomedical Technology Co., Ltd, Beijing, China) according to the manufacturer's instructions. The quality of RNA was evaluated by Agarose Gel Electrophoresis, NanoDrop (NanoDrop Technologies, Wilmington, DE, USA), and Agilent 2,100 Bioanalyzer (Agilent, Santa Clara, CA, USA) to ensure the use of qualified samples for sequencing. Subsequently, ribosomal RNA (rRNA) was depleted from total RNA using the Epicentre Ribo-zero™ (rRNA) Removal Kit (Epicentre, Madsion, WI, USA). Preparation of the RNA sequencing library construction was prepared as previously reported (22), that is, RNA library for RNA-seq was prepared as rRNA depletion and stranded method. Then the library concentration was preliminarily tested using Qubit 2.0, and was adjusted to 1.0 ng/μL. Afterwards, the Agilent 2,100 Bioanalyzer (Agilent, Santa Clara, CA, USA)was used to measure insert sizes of the acquired library. To ensure library quality, quantitative real-time PCR (qRT-PCR) was used to accurately determine the library effective concentration (>3 nM). After qualification of the library, the Illumina PE150 platform was used for sequencing according to the library effective concentration and the pooling of data output requirements.

The raw data (raw reads) containing a small number of reads with sequencing adapters or with low sequencing quality was obtained by RNA-seq. For the reliability of the data analysis, it was necessary to filter the raw data using the following parameters: (1) reads with adapter were removed; (2) unidentifiable base information reads with a ratio of > 0.002 were removed; (3) the paired reads were removed when the number of low-mass bases in single-end read was more than 50% of the read length ratio. The obtained clean reads were mapped to the internal reference genome (Bos taurus Ensembl 97) by using Hisat2 (23).

In this study, all the transcripts were merged using Cuffmerge software. LncRNAs were identified according to their structural and functional characteristics. The criteria used to identify the lncRNAs were as follows: (1) transcripts with exon numbers ≥ 2 and lengths > 200 bp were selected; (2) transcript fragments per kilobase per million reads (FPKM) > 0.5 were selected; (3) transcripts that overlap with the database annotation exon regions were identified using the Cuffcompare software (v2.2.1); (4) candidate data sets of novel lncRNAs were obtained through potential coding identification. Then, the final identification and naming were performed to obtain the novel lncRNAs of this analysis according to their loci relative to the coding genes (24).

Quantification of lncRNAs was performed using the StringTie software (v1.3.3). The expression level of each transcript was calculated according to the frequency of clean reads. The data were subsequently normalized to the FPKM (25), using the following formula:

Differential expression analyses between two groups were conducted using the Cuffdiff software (v2.1.1). The DElncRNAs were identified using |log2(FoldChange)| > 1 and P < 0.05 as standard.

To verify the sequencing results, 10 DElncRNAs were randomly selected to detect their relative expression levels by qRT-PCR. Total RNA extracted by using the Trizol reagent was reverse-transcribed into cDNA using the PrimeScriptTM RT reagent Kit with gDNA Eraser (Takara Biomedical Technology Co., Ltd, Beijing, China) according to the manufacturer's instructions. The qRT-PCR was performed using 2 × M5 HiPer SYBR Premix Es Taq (with Til RNaseH) kit (Mei5 Biotechnology, Co., Ltd, Beijing, China) on a Bio-Rad CFX96 qRT-PCR system (Bio-Rad, Hercules, CA, United States) to analyze the relative expression level of DElncRNAs. The qRT-PCR reaction consisted of the following components: 10.0 μL 2 × M5 HiPer SYBR Premix Es Taq (with Til RNaseH), 0.4 μL Forward Primer, 0.4 μL Reverse Prime, 1.5 μL cDNA template and 7.7 μL RNase-Free ddH2O. The qRT-PCR amplification conditions were as follows: at 95 °C for 30 s, followed by 40 cycles of denaturation at 95 °C for 5 s followed by annealing/extension at 60 °C for 30 s. The qRT-PCR primers for DElncRNAs were designed with Primer Premier 5 (PREMIER Biosoft international, Palo Alto, CA, USA). The primer sequences are presented in Supplementary Table 1. For each sample, 3 biological replicates were amplified. Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) and ribosomal protein S18 (RPS18) genes were used as internal references (26), and the relative expression of DElncRNAs was calculated using the 2^−ΔΔCT method.

The qRT-PCR and RNA-seq data were compared with each other. When the expression trends of qRT-PCR and RNA-seq were consistent, the DElncRNA identification was considered to be reliable.

To explore the potential function of DElncRNAs, target genes of the DElncRNAs were predicted as being either cis-acting and trans-acting. For cis-acting target (co-location) prediction, genes that co-localized to within 100 kb up-or down-stream of a lncRNA were considered to be a cis-acting target (27). For trans-acting target (co-expression) prediction, the Pearson's correlation coefficients between the coding genes and lncRNAs were calculated and analyzed for the identification of trans-acting regulatory elements. This analysis required the sample size was > 5, and a Pearson's correlation coefficient > 0.95 for a gene to be considered a target of the lncRNA (28). Then, Gene Ontology (GO) functional enrichment analysis of target genes of DElncRNAs was implemented by the cluster Profiler R package (R-3.2.4), in which gene length bias was corrected (29).

All experiments were repeated at least 3 times. Significant differences were analyzed using Student's t-test, and P < 0.05 was considered to be a significant difference and P < 0.01 was extremely significant. The experimental data were presented as mean ± SEM.

Nine lncRNA libraries were sequenced using the Illumina PE150 to acquire raw data for quality control, and to extract high-quality clean reads (Table 1). In addition, the Hisat2 software (v2.0.5) was used to compare the clean reads to the reference genome for statistical analysis. For each sample, the clean bases (G) should be ≥ 11.29 G and the Q30 (%) should be ≥ 92.78% to ensure the reliability of results obtained from subsequent analyses. Additional details are provided in Table 1.

In this study, the correlation of expression levels between samples was determined to ensure the confidence and rationality of the experiment (Figure 1A). Through analysis of the sequencing data, a total of 2,257 lncRNAs were identified, including 210 annotated and 2,047 novel. Novel lncRNAs were classified as intergenic lncRNA (lincRNA, 29.8%), antisense lncRNA (antisense, 28.9%), and sense overlapping (41.3%) (Figure 1B).

To further detect the expression of lncRNA in LPS treated bMECs, cells were incubated for 0, 6, and 12 h following treatment for high-throughput sequencing. In each sample, box-plot diagrams were constructed based on FPKM to observe the expression symmetry and distribution of all lncRNAs. The results showed that the distribution of each box-plot was flat (Figure 2), which indicated that the expression distribution of lncRNA in each sample was correct. At the same time, it can be seen that the distribution of lncRNA expression levels in samples between groups were similar, but small differences were observed.

A total of 112 DElncRNAs were observed between the 3 experimental groups. Relative to the 0 h group, there were 22 up- and 25 down-regulated DElncRNAs in the 6 hpi group (Supplementary Table 2; Figure 3A), and 27 up- and 22 down-regulated DElncRNAs in the 12 hpi group (Supplementary Table 3; Figure 3B). Compared with the 6 hpi group, 32 up- and 25 down-regulated DElncRNAs were observed in the 12 hpi group (Supplementary Table 4; Figure 3C). Moreover, among the 3 groups, a total of 17, 23 and 31 specific DElncRNAs were observed in 6 vs. 0 hpi, 12 vs. 0 hpi, and 12 vs. 6 hpi groups by Venn diagram analyses, respectively (Figure 4). It is worth noting that lncRNA TCONS_00039271 was highly expressed and significantly down-regulated (log2(FoldChange) = −1.9) at 12 hpi relative to 0 h. Meanwhile, lncRNA TCONS_00139850 was not detected in the 0 h and 6 hpi group and was significantly up-regulated (log2(FoldChange) = 15.94) in the 0 h and 12 hpi group. The results indicated that these DElncRNA may play an important role in the regulation of the inflammatory response following LPS treatment of bMECs.

To verify the accuracy and reliability of the RNA-seq data, a total of 10 DElncRNAs were randomly selected for expression level analysis by qRT-PCR. Comparing the qRT-PCR results with the RNA-seq results, as shown in Figure 5, the data obtained by the two methods yielded similar patterns, indicating that the screened DElncRNAs based on RNA-seq were reliable.

To further investigate the physiological and molecular functions of DElncRNAs on inflammatory reactions in bMECs, or more generally, the pathogenesis of bovine mastitis, possible target genes of DElncRNAs were predicted via co-location and co-expression methods (27, 28). The prediction results of target genes showed that 112 DElncRNAs may have 1,042 target genes. The results of GO enrichment showed that these target genes were involved in the regulation of various biological processes such as inflammatory response, immune response, stress response and regulation of cytokine production (Figure 6). Hence, it can be speculated that these DElncRNAs may participate in the inflammatory response by regulating their target genes during LPS-induced inflammatory response of bMECs.

It is worth noting that TCONS_00039271 and TCONS_00139850 were stably down- and up-regulated, respectively. Constructing the network between lncRNAs and target genes indicated that each DElncRNA possibly targets approximately14 genes (such as NOTCH2, TNFSF14, CDC42, GNG2, RHEB and ITGA5) (Figure 7; Supplementary Table 5). Previous studies have shown that these target genes involved in inflammatory pathways including the Notch (30), NF-κB (31). MAPK (32), chemokine (33), mTOR (34) and PI3K-Akt signaling pathways (34, 35). In addition, Li et al. (30) and Zhao et al. (36) respectively reported that lncRNA DLEU2 and lncRNA LINC01410 can elevate the expression of NOTCH2 to stimulate Notch signaling pathway, thereby boosting cell proliferation and hindering cell apoptosis. Taken together, the data suggest that lncRNA TCONS_00139850 and TCONS_00039271 might play a role in regulating the inflammatory response of bMECs via their target genes (Figure 7; Supplementary Table 5). As such, these lncRNAs and their target genes will be the focus of future research on the molecular mechanism of mammary inflammation-related lncRNAs.