PropertyValue
is nif:broaderContext of
nif:broaderContext
is schema:hasPart of
schema:isPartOf
nif:isString
  • Holstein dairy cows were selected from two herds, one MAP free (verified by ELISA and absence of clinical cases) and the other positive for MAP (verified by ELISA, fecal culture and presence of clinical cases). From the positive herd, 5 positive (PP) cows based on ELISA test results and 5 ELISA negative, potentially exposed (NP) animals were chosen. All animals were 4 to 5 years old, had a body condition score (BCS) of 3 and were at 170 to 190 days in milk (DIM) when sampled. Positive subjects were confirmed by fecal culture. Five negative non-exposed control animals (NN) were selected from the negative herd, matched for age, lactation, BCS and DIM with the positive and exposed animals. For the miRNA-Seq experiment, in addition to these animals, 7 negative cows from the same negative herd were added to enlarge the control group. Sample preparation, RNA extraction and quality control: Samples used were taken during obligatory routine animal sanitary controls by an authorized veterinarian. No ethical approval was required, in compliance with the European Directive 2010/63/UE and the Italian regulation D. Lgs n. 26/2014. Whole blood was collected into PAXgene® Blood RNA Tubes (PreAnalytiX GmbH) and kept at room temperature for at least two hours before freezing at -20°C for 24 hours before being transferred at -80°C until processing, as recommended by the manufacturer. Total RNA was extracted using the PAXgene® Blood miRNA Kit (PreAnalytiX GmbH) according to the manufacturer’s protocol. Total RNA was eluted in a final volume of 80 μL. RNA concentration was measured by NanoDrop™ 1000 spectrophotometer (Thermo Scientific) and RNA integrity was assessed with an Agilent 2200 TapeStation system (Agilent Technologies). RNA-Seq library preparation and sequencing: Libraries were prepared with the Illumina Truseq RNA sample prep kit (Illumina Inc., USA) following manufacture’s protocol and the size and yield evaluated using an Agilent TapeStation 2200. Libraries were then quantified with an ABI9700 qPCR instrument using the KAPA Library Quantification Kit in triplicate, according to the manufacture’s protocol (Kapa Biosystems, Woburn, MA, USA) and then normalized to 10 nM as recommended by Illumina for cluster generation on the Hiseq2000. Fifteen libraries were prepared and equimolar amounts of 5 samples were mixed before NaOH denaturation. Each of the pools was run in a lane of a Hiseq Flowcell. The Illumina Truseq PE cluster kit v3 was used to generate clusters on the grafted Illumina Flowcell and the hybridized libraries were sequenced on a Hiseq2000 with a 100 cycles of paired-end sequencing module using the Truseq SBS kit v3 (Illumina Inc., USA). miRNA-Seq library preparation and sequencing: For each sample, 5 μl of RNA were used to prepare a library with the TruSeq SmallRNA kit (Illumina Inc., USA) following the manufacturer’s instructions. In order to minimize primer dimers formation, half of the TruSeq Small RNA sample reagents were used, followed by 11 PCR cycles of PCR to amplify the library. 10 μl of unique indexed libraries were pooled and DNA fragments from 140 to 160 bp (the length of miRNA inserts plus the 3′ and 5′ adaptors) were selected from a Pippin 3% Agarose Dye free Gel cassette (BluPippin, Sage Science, MA, USA), which were then recovered in 40 μL of Pippin elution buffer. Fragments were purified by Qiagen MinElute PCR Purification kit (Qiagen, CA, USA). The indexed libraries were quantified as described above. Two pools of 11 sample libraries were prepared and 10 μL of each pool at a final concentration of 2 nM were used in a lane for 50bp Single-Read sequencing for a total of 2 lanes of an Illumina HiSeq2000. Preliminary quality control of raw reads was carried out with FastQC software v0.11.2 [89]. Illumina raw sequences were trimmed using Trimmomatic [90] and PCR primers and Illumina adapter sequences were removed. Minimum base quality 15 over a 4 base sliding window was required, then only sequences longer than 36 nucleotides were included in the downstream analysis. Reads which successfully passed trimming were mapped against the Bos taurus UMD3.1.68 reference genome sequence, using STAR [91] aligner to obtain BAM alignment files. The BAM files were sorted and indexed using Samtools [92]. In order to quantify counts for each sample a list of genes and relative abundance of mapping reads were extracted using htseq-count [93]. These count files were used for downstream statistical analysis. Trimming and quality control were performed as described for RNA-Seq analysis, with the difference that a 15 nucleotide minimum sequence length was required. Reads which passed the quality control were used for novel small-RNA discovery using Mirdeep2 [43]. Both “cow” and “human” known small-RNAs (mature and precursors) were downloaded from MirBase [94] and used as support datasets to help the discovery process. Absolute positions of novel miRNAs on the bovine reference genome were retrieved by BLAST [95]. Counts of all the known and the novel miRNAs were used to quantify expression levels for each sample using Mirdeep2. This pipeline produced a list of small-RNA IDs and the relative abundance of mapping reads (counts) for each sample which was used in the downstream statistical analysis. Expression Analysis of mRNAs and miRNAs: Statistical analyses to compare mRNA and miRNA expression profiles were performed using the "R" statistical environment edgeR [96], vegan [97] and gplots [98] packages. Genes that had at least one count per million in at least 3 samples were included in the gene expression analysis. A general linear model was used in the edge R Package to generate lists of mRNAs and miRNAs with statistically significant different expression among the three comparisons: PP vs NN, NP vs NN and PP vs NP. Differentially expressed mRNAs and miRNAs were defined as having a False Discovery Rate (FDR) below 0.05. Real-time RT-qPCR for RNA-Seq validation: To validate the RNA-Seq gene expression data RT-qPCR was performed on 4 genes selected from those which were differentially expressed in both PP vs NN and NP vs NN comparisons and had a FDR<0.01. Three hundred ng of total RNA were reverse transcribed into cDNA using the QuantiTect Reverse Transcription Kit (Qiagen), RT-qPCR was performed using 5 ng cDNA, and 250 nM each primer in GoTaq qPCR Master Mix (Promega). The reaction mixtures were incubated in 384-wells plates at 95°C for 10 min, followed by 45 cycles of 95°C for 15 seconds and 62°C for 1 min using a CFX 384 Real-Time PCR Detection System (Bio-rad Laboratories, USA). All reactions were performed in triplicate and “no-template” controls were included. Following amplification, a melting curve analysis was performed to verify the specificity of the reactions. GAPDH and ACTB were chosen as the reference genes, as they were stably expressed in the RNAseq data, and assayed in the same samples to normalize the data. In order to determine the efficiency and the dynamic range of the reaction, for each primer pair a standard curve was constructed from triplicate assays for 4 dilutions from 30 ng to 0.03 ng of pooled cDNA samples. Primers for Real-time RT-qPCR are listed in Table 7 and were designed using Primer-BLAST (NCBI) or deduced by the literature. The relative expression level of each selected gene was calculated according to the Pfaffl method [99] and was reported as normalized fold expression relative to the control (NN cows). Log2 fold-change data of RT-qPCR and RNA-Seq analyses were compared to validate results. The Pearson correlation coefficient between the two analyses was calculated using IBM Statistical Package for Social Sciences software (SPSS, ver. 21; IMB SPSS Inc., Chicago, IL, USA) and differences were considered as statistically significant if the p-Value was <0.05. Table data removed from full text. Table identifier and caption: 10.1371/journal.pone.0164461.t007 Primers used for RT-qPCR. Functional Analysis of gene expression data: The RNA-Seq differentially expressed (DE) gene lists obtained for each comparison (PP vs NN, NP vs NN, PP vs NP) were submitted to the Qiagen Ingenuity Pathway Analysis (IPA; Ingenuity Systems Inc., USA). The DE genes were filtered for an FDR≤0.05 and a fold change ≥1 or ≤-1 and were used as input in three separate data sets. Each gene was mapped to its corresponding gene object in the Ingenuity Knowledge Base. In the Core Analysis, IPA assigns functional information and biological relevance by analyzing RNA expression data in the context of known biological response and regulatory networks. In addition, canonical pathways, biological functions, and networks were investigated to identify major biological pathways associated with MAP infection. The upstream regulators were investigated based on an activation z-score, which was considered as significant when below -2 (inhibited) or above 2 (activated). Prediction of miRNA target genes and correlation with mRNA data: Due to the uncertainty associated with the prediction of miRNA target genes [40], two approaches were used. For known miRNAs only the gene targets identified by both methods were considered in the subsequent analysis. miRDB software [101] was used with default threshold settings for the initial prediction, based on miRNA-target matching the miRNA seed region and 3’ UTR mRNAs. The method was applied for both known and novel differentially expressed miRNAs and a threshold of 90 was set to select a panel of potential target genes. Predicted targets of known miRNAs were then analyzed with the TargetScan software v7.0 [102]. RNA-Seq expression levels were investigated for the genes that were shared between the two lists. Predicted gene targets of the novel miRNAs were obtained only from the miRDB prediction software, because target identification was done directly from the novel miRNA sequences. Targets which in RNA-Seq results had FDR <0.05 and opposite pattern of expression with respect to related miRNAs were considered as putatively affected by miRNA expression.
rdf:type