Andrew Marete1†
Olivier Ariel
Eveline Ibeagha-Awemu
Nathalie Bissonnette- 1Agriculture and Agri-Food Canada, Sherbrooke Research and Development Centre, Sherbrooke, QC, Canada
- 2Faculty of Science, Sherbrooke University, Sherbrooke, QC, Canada
Mycobacterium avium ssp. paratuberculosis (MAP) causes chronic enteritis in most ruminants. The pathogen MAP causes Johne's disease (JD), a chronic, incurable, wasting disease. Weight loss, diarrhea, and a gradual drop in milk production characterize the disease's clinical phase, culminating in death. Several studies have characterized long non-coding RNA (lncRNA) in bovine tissues, and a previous study characterizes (lncRNA) in macrophages infected with MAP in vitro. In this study, we aim to characterize the lncRNA in macrophages from cows naturally infected with MAP. From 15 herds, feces and blood samples were collected for each cow older than 24 months, twice yearly over 3–5 years. Paired samples were analyzed by fecal PCR and blood ELISA. We used RNA-seq data to study lncRNA in macrophages from 33 JD(+) and 33 JD(–) dairy cows. We performed RNA-seq analysis using the “new Tuxedo” suite. We characterized lncRNA using logistic regression and multilayered neural networks and used DESeq2 for differential expression analysis and Panther and Reactome classification systems for gene ontology (GO) analysis. The study identified 13,301 lncRNA, 605 of which were novel lncRNA. We found seven genes close to differentially expressed lncRNA, including CCDC174, ERI1, FZD1, TWSG1, ZBTB38, ZNF814, and ZSCAN4. None of the genes associated with susceptibility to JD have been cited in the literature. LncRNA target genes were significantly enriched for biological process GO terms involved in immunity and nucleic acid regulation. These include the MyD88 pathway (TLR5), GO:0043312 (neutrophil degranulation), GO:0002446 (neutrophil-mediated immunity), and GO:0042119 (neutrophil activation). These results identified lncRNA with potential roles in host immunity and potential candidate genes and pathways through which lncRNA might function in response to MAP infection.
Introduction
One of the most economically significant diseases in livestock is paratuberculosis (1). The etiological agent of paratuberculosis or Johne's disease (JD) is Mycobacterium avium ssp. paratuberculosis (MAP). Whitlock and Buergelt described JD as chronic, wasting, incurable, and infectious ruminant enteritis (2). For dairy producers, MAP infection translates to significant financial losses related to reduced milk production, decreased pregnancy rates, increased replacement costs, and decreased slaughtered carcass weight (3, 4), not to mention diminished animal welfare. Around 24–66% of dairy herds in Canada are MAP infected (5). Horizontal transmission of infection via the fecal–oral route is the most important mode of spread of infection due to the high amounts of MAP excreted in the feces. After ingestion of contaminated water or food, MAP reaches the gastrointestinal tract; MAP shows an evident tropism for this site (1, 6). The first route of MAP entry is through the ileum and jejunum's organized lymphoid tissue, the Peyer's patch in the intestinal mucosa and submucosa (7). These early events of MAP infection occur in two functional stages: (1) invasion through the intestinal barrier via MAP discharge from epithelial M cells and (2) phagocytosis and survival in macrophages of the lamina propria (7–9). It is known that MAP uses tissue-resident macrophages as its primary reservoir for survival and multiplication (10–12). Interestingly, genetic variations in numerous candidate genes expressed in macrophages are associated with resistance/susceptibility to MAP infection, notably the NOD2 (13, 14), IL10 (15–18), SLC11A1, and Toll-like receptor genes (19, 20).
With a slow progression of the disease, the pathogenesis of JD makes diagnosis difficult, more so in the subclinical stage of infection before the clinical signs appear (21). The first clinical signs include gradual weight loss despite normal appetite or, sometimes, increased appetite. During this clinical period, there is a decrease in milk production, accompanied by a concomitant weight loss with a sometimes more fluid consistency of manure. Diarrhea may be intermittent at first, with periods of normal manure consistency leading to chronic diarrhea (2). During the prolonged incubation subclinical period of 4–7 years, scarce clinical signs are observed (22). Difficulty in JD diagnosis is further hindered by host genetics (23), herd management, MAP strain, and infectious dose (24).
MAP employs complex mechanisms to control macrophages, which turn into a duel that lasts for years with unpredictable disease progression. The T helper type 1 (Th1)-mediated response that usually effectively controls non-pathogenic intracellular mycobacterial infections fails for MAP infection (25). The pathogenesis of JD is still under investigation because macrophage-MAP cross talk in the subclinical stage of the disease is partially resolved. It is paramount to consider the study of alternative molecular avenues for identifying the product resulting from bacterial MAP infection that might become potential biomarkers of JD and evolve therapeutic tools.
Previous reports indicate that bacteria interfere with mammalian regulatory RNA expression that is not translated to protein (such as long non-coding RNA–lncRNA) to modify immune signaling, autophagy, or apoptosis machinery (26, 27). lncRNAs are now emerging as important regulators of innate and adaptive immune responses (28–30). In humans, while they are widely investigated in aging, cancers, and epigenetics (31–33), there is growing evidence that lncRNAs interfere in the pathogenetic mechanisms of multifactorial disease like Crohn's disease and inflammatory bowel disease (34, 35).
In bovine, few studies have examined the occurrence of lncRNA in tissues, including muscle, skin, various tissues, and the mammary gland (36–40). Previous studies also report that lncRNA is involved with host cell response toward bacterial infections, including paratuberculosis (27, 41). LncRNA is unique from other RNA based on size (>199 nucleotides) and limited evidence of protein-coding potential (36, 37, 42–44). However, lncRNA's novel nature means there is no consensus on the best way to classify the protein-coding and non-coding potential of the lncRNA, so we use both logistic regression (45) and multilayered neural networks (46). The former implements human-designed features, such as open reading frame (ORF) length and integrity, GC content, and hexamer usage bias, whereas the latter identifies multilayered deep patterns solely on sequence information.
This study aims to use available deep learning and logistic regression approaches to study lncRNA associated with MAP in Canadian Holstein and provide novel insight into lncRNA's regulatory function in macrophages of dairy cattle during MAP infection. To this end, we investigated the presence of potentially novel lncRNA candidates and their role in MAP infection using RNA sequencing.
Materials and Methods
Cow Selection and Johne's Disease Diagnosis
According to the Canadian Council on Animal Care guidelines for institutional animal use, we carried out all animal procedures and obtained ethical approval for the study from the Agriculture and Agri-Food Canada Animal Ethics Committee (protocol 362). To select cows for RNA-Seq, we analyzed fecal and blood samples from 15 commercial dairy herds positive for JD located in the province of Quebec, Canada, as described in the companion project (47).
Briefly, from each herd, we sampled cows twice yearly. The cows were older than 24 months to be enrolled in the study. They had calved twice or more at culling. We collected one fecal sample of a volume equivalent of 100 mL using a single-use veterinary glove. Consecutively, we also drew two blood samples per cow in dry tubes for serum collection (SST Serum Separation Tubes 8.5 mL; BD Biosciences, Ontario, Canada). Within 1 h of sampling, the tubes were centrifuged at 1,300 × g at 20°C for 10 min and kept at 4°C during the transport to the laboratory. Sera were collected and then stored at −80°C until ELISA analysis. According to the manufacturer's instructions, we processed sera using the IDEXX MAP Ab test kit (IDEXX Laboratories, USA). As described by Collins (48), we transformed optical density values into sample-to-positive (S/P) ratios and selected samples with an S/P ratio of at least 55% as positive. The presence of MAP in feces was tested by qPCR, and JD cows were confirmed infectious using the BD MGIT ParaTB culture medium and the BACTEC 960 detection 960 system described in Fock-Chow-Tho et al. (47). Cows that presented concordant serological and fecal culture or qPCR statuses, either positive or negative, were retained. A cow was designated JD (+) when a fecal culture and ELISA were positive at least two sampling periods. During macrophage analysis, the mean age of JD(–) cows was 6.4 ± 1.5 years, and 5.1 ± 1.8 years for the JD(+) cows. While JD(+) cows were promptly culled, the JD(–) cows were kept on-farm for >7 years to confirm their status definitively. In total, we selected 66 cows for RNA-Seq analysis, of which 33 were JD(+), and 33 were JD(–).
RNA Isolation, Library Preparation, and Sequencing
The monocyte-derived macrophages (MDM) were prepared in the absence of FBS and granulocyte-macrophage colony-stimulating factor (GM-CSF), or M-CSF, to avoid activation or bias in the differentiation toward M1 or M2 polarization, as described in Ariel et al. (49). Freshly isolated monocytes were confirmed exempt of MAP from both JD(+) and JD(–) cows, confirmed using qPCR and fluorescent microscopy. DNA was extracted from adherent monocytes, MDM, and PBMC using ZR Fecal DNA MiniPrep kit (Zymo Research Corp., Irvine, CA, USA), and qPCR was performed using the VETMAX Gold MAP Detection Kit (Life Technology Inc., Burlington, Ontario, Canada) as described previously (47). The absence of MAP in JD(+) and JD(-) MDM was also confirmed in vitro using fluorescence microscopy as described (49). To profile the lncRNA in resting macrophages (CTL, i.e., non-infected) and in response to MAP infection, the MDM were also infected with MAP. Our previous experimental design was used: 1 h, 4 h, 8 h, and 24 h post-infection with MAP at the multiplicity of infection of 10 (49).
In summary, we extracted total RNA from MDM from each experimental time point (CTL and MAP-infected at 1, 4, 8 h, and 24 hpi) in 66 cows [33 JD(–) and 33 JD(+)] using the RNeasy kit (Qiagen) total RNA extraction protocol. We quantified the RNA yield using a NanoDrop spectrophotometer (Thermo Fisher) and assessed RNA quality using the Bioanalyzer RNA 6000 kit (Agilent Technologies). We used the Ribo-Zero Gold kit to remove ribosomal RNA and Illumina TruSeq Stranded Total mRNA Sample Preparation kit (Illumina) to generate cDNA libraries. After quality control (size and absence of primer dimers) and qPCR library quantification [Kappa Library Quantification kit (Roche)], we performed paired-end sequencing using the Illumina HiSeq 2500 platform running HiSeq Control Software (v2.2.68). A subset of the sequencing data from the 66 cows is available from the Gene Expression Omnibus repository (accession number GSE98363). All processes followed manufacturer recommendations.
Transcriptome Assembly, Novel lncRNA Prediction, and Differential Expression (DE) Analysis
Figure 1A illustrates the steps of transcriptome assembly. Illumina adapter sequences were trimmed from each RNA-Seq read using Trimmomatic (V0.39), keeping reads longer than 36 bp and with a Phred score ≥30. Reads were mapped to UMD 3.1.1 bovine genome assembly using HISAT2 (v2.2.0) and transcripts assembled using StringTie (v2.1.0). We merged assembled transcripts from all cows using the—merge option of StringTie, resulting in non-redundant assembled transcripts. Using Gffcompare (v0.10.1), transcripts were then compared with Ensembl bovine gene annotation (release 94) to identify transcripts overlapping with known protein-coding and non-coding regions. To identify lncRNA, we used the transcript classification codes of Gffcompare to select transcripts categorized as “u” and with a length of ≥200 nt.
Figure 1. (A) The study workflow from raw fastq to lncRNA identification and DE of the nearest genes. (B) Two-graph receiver operating characteristic curve (ROC) to determine the optimal coding-probability cutoff value. (C) Combinatorial effects of Fickett score, hexamer score, and ORF size on coding transcripts (brown dots) and non-coding genes (blue dots). (D) Mapping statistics (vertical axis) showing unique (red), multi-mapped (yellow-green), and overall (blue) alignment rate for all cows (horizontal axis). (E) Snapshot of lncRNA statistics identified by two tools with known lncRNAs in parentheses.
Transcripts were analyzed using two approaches: (1) estimation of transcript coding probability and (2) differential expression (DE) analysis. Coding probability was estimated using two tools: RNAsamba (46) and CPAT (v2.0.0) (45). Both tools were tested using the Bos taurus dataset of known and unknown protein-coding sequences to train the models. RNAsamba computes RNA sequences' coding potential using a neural network classification model resulting in sequences classified as coding or non-coding based on an estimated coding score. The approach of CPAT uses a logistic regression model. Also, CPAT evaluates each base's unequal content frequency and asymmetrical distribution in the positions of codons in one sequence, i.e., the Fickett score and usage bias of adjacent amino acids in proteins, namely, the hexamer score. The models' respective outputs were evaluated using 20-fold cross-validation to determine the coding probability cutoff (Figures 1B,C). Using usearch (50), transcripts with a coding probability of ≤ 0.4 and ORF ≤ 300 bp were selected. Sequences were further filtered out if they blasted against the Swiss-Prot database (e < 1 × 10−05). The final dataset was compared to the NONCODEV5 database (51). Those annotated were qualified as known bovine lncRNA, and the remaining isoform transcripts (code = “j”) were classified as novel lncRNA.
We used DEseq2 (v1.26.0) for DE analysis using raw read counts of each sample from the final retained dataset. DESeq2 calculates each sample size factor to correct for library size and RNA composition bias (52). We considered lncRNA as truly expressed if normalized counts ≥5 in at least 10% of our libraries and an FDR < 0.5. To explore the functions of significant DE lncRNA, we used bedtools to obtain the genes 100 kb of each lncRNA. We performed gene ontology (GO) term enrichment analysis using the Panther classification system (53) and Reactome pathways (54). Significance was expressed as P-value, with a lower P-value indicating higher significance.
Results and Discussion
In most dairy cows, JD progression is slow, primarily due to the ability of MAP to lodge in intestinal tissue-resident phagocytic cells and escape the immune system's surveillance. Myeloid cells, including monocytes, macrophages, and neutrophils, work in concert with lymphoid cells to initiate and amplify innate and adaptive immunity. Previous reports indicate that lncRNA plays a significant role in regulating the immune response toward several bacterial pathogens known to induce Mycobacterium tuberculosis infection in human macrophages (55). We hypothesized that lncRNA could be part of the mechanisms employed by MAP to control macrophages. The current study provides information (1) on potentially lncRNA-targeted genes affected by MAP infection to weaken the host and (2) on the biological pathways that might lead to susceptibility to MAP infection. In our study, RNA-Seq data from macrophages of 33 cows diagnosed JD positive (+), and 33 JD negative (–) cows were used, among which the study of differential gene expression of 12 cows was previously described (49).
Prediction of lncRNA From Expressed in Bovine Macrophages
To investigate the potential role of lncRNA in macrophages for JD susceptibility, we used RNA-Seq data of macrophages from JD(+) and JD(–) cows. Previous studies have demonstrated the importance of a strong correlation between read alignments and accurate transcript assembly and quantification since misaligned reads usually decrease the number of reconstructed genes (56). In our previous study (49), we analyzed the differentially expressed (DE) genes in macrophages from JD(+) and JD(–) cows using a TopHap-Cufflink pipeline. In the current study, we used HISAT2 because this program, while aligning RNA-Seq reads to the genome, discovers transcript splice sites and provides an accurate representation of all transcript isoforms, creating a rigorous representation of lncRNA (56).
Identifying and inferring lncRNA's biological role is challenging, more so for dairy cattle, where the functional annotation of lncRNA is limiting (57). This study implemented a computational pipeline based on the “new Tuxedo” package (56). As illustrated in Figure 1A, we used HISAT2, StringTie, and Gffcompare to align transcripts to the reference genome, assemble the transcripts, and produce transcript statistics, thus allowing us to obtain results comparable to previous lncRNA studies in cattle (38, 40, 58). Furthermore, we obtained promising alignment statistics (>92% concordant alignment for all samples, Figure 1D) with HISAT2, which was more accurate compared to its preceding software (59, 60), and exhibited a faster search algorithm (61). After mapping a minimum of 13 million paired reads per cow to the UMD 3.1.1 bovine reference genome, on average, 71% were uniquely aligned, 14% were multi-mapped, and the average overall alignment rate was 92% (Figure 1D). Of the lncRNAs identified in this study, 45.05% had an average length of 600 bp long with a range of 200–1,000 bp, most of which were mostly intergenic. This result concurred with previous reports that lncRNAs are mainly located between genes (i.e., intergenic) with a smaller overlap within genic regions (42, 44). Moreover, though previous studies indicate intergenic lncRNA may act in cis or trans to regulate gene activities (62, 63), the ever-continuing curation of the bovine functional lncRNA annotation posed a challenge to study how the identified lncRNA may act in trans to regulate distant genes.
With the low multi-mapped read rate and high overall alignment rate, we were confident that false-positive alignments would not disrupt StringTie's flow algorithm and skew the expression estimates of assembled transcripts. Parsed mapped files identified 47,683 potential transcripts, of which 13,301 were putative lncRNA (i.e., had a code =“u”). CPAT and RNAsamba predicted 3,894 and 6,593 as non-coding transcripts, respectively, with an intersection of 2,814 non-coding transcripts (Figure 1E). Most of the lncRNA were <1,000 base pairs, with the shortest lncRNA being on chromosomes 12, 26, 27, and X, with ENSBTAG00000046640, LBX1, ERI1, and AMELX being the closest genes, respectively (Table 1). Chromosomes X and 18 report the highest number of lncRNAs identified in bovine macrophages, with 1,385 and 1,054 lncRNA, respectively (Figure 2).
Table 1. Location of the highest DE lncRNAs in macrophages from JD(+) vs. JD(–) cows identified on Bos taurus autosomes 1–29 and chromosome X.
Figure 2. Distribution of putative, known, and novel lncRNAs in chromosomes 1–29 and X.
Novel lncRNA in Bovine Macrophages
Among the 16,970 lncRNA predicted by the two tools, 3,669 were novel non-coding transcripts, and 385 were known lncRNA. Interestingly, most novel lncRNAs identified in bovine macrophages were mapped on chromosomes X and 18 (Figure 2 and Supplementary File 1). Chromosomes 26 and 28 had the least number of novel non-coding transcripts, with 148 and 161 lncRNAs, respectively (Supplementary File 1). CPAT predicted a smaller number of overall non-coding transcripts, but 35% of these were novel, whereas RNAsamba predicted a larger number of non-coding transcripts, of which 33% were novel lncRNA (Figure 1E).
Differential Expression (DE) of lncRNA in Macrophages From JD(+) and JD(–) Cows
The differential expression (DE) analysis method usually has the most substantial impact on results (64, 65). For this study, we used DESeq2 to study DE of lncRNA detected in macrophages from JD(+) and JD(–) cows. DESeq2 performs robustly in comparison to other existing DE tools (66). As presented in Table 1 and Figure 3, DE analysis using the two predictor tools identified the lncRNA having the greatest significance (P < 0.05) for all chromosomes. The two longest DE lncRNAs were on chromosomes 4 and 10, within a 0.5-Mb region (P = 3.45 × 10−03), at 49 Kb from the FZD1 gene, and within a 0.1-Mb region (P = 4.45 × 10−04) close to the FLRT2 gene, respectively (Table 2 and Figure 3). Interestingly, the Frizzled Class Receptor 1 gene (FZD1) was DE in our previous study, where fold change was estimated using FPKM: FZD1 was 15.56 times more expressed in macrophages from JD(+) cows than JD(–) macrophages (49). The gene encodes a transmembrane domain protein acting as a receptor for Wnt signaling proteins, which are essential for regulating pro-inflammatory cytokines in response to bacteria and mycobacterial infection (67).
Figure 3. Normalized count distributions in JD(+) macrophages (purple) and JD(–) macrophages (blue) for 10 highly expressed genes close to lncRNA transcripts.
Table 2. The top 10 most significant DE lncRNAs in macrophages ranked by the FDR corrected P-value.
Ten of the predicted lncRNAs were highly DE at P < 0.05 (Table 2). The most highly DE lncRNA was downregulated by ~3-fold change (FC) (Log2 −1.58 ±0.4) in JD(+) compared to JD(–) macrophages. This lncRNA is located within a 0.4-Kb region on chromosome 24 and is close to the predicted twisted gastrulation protein homolog one gene (TWSG1) (Table 2). TWSGI was found expressed in bovine macrophages and, interestingly, was previously significantly downregulated by 0.70 Log2FC in response to MAP infection at 4 h post-infection (4 hpi) (49). The second highest DE was an lncRNA of 0.3 Kb located on chromosome 20 with a Log2FC of −1.23 (±0.4). The negative Log2FC indicates that macrophages from JD(+) cows have 2.35 times fewer lncRNA transcripts than macrophages from JD(–) cows. This lncRNA is 13 Kb away from Bos taurus biorientation of chromosomes in cell division 1 gene (BOD1). Although BOD1 is a provisional NCBI gene (accession no. NM_001076200), BOD1 expression was detected in bovine macrophages (49) but was however not found DE between JD(–/+) groups. Two lncRNAs with the greatest Log2FC are located in a 31-Kb region of chromosome 9 (FC = 2.35, P = 6.65 × 10−04, vicinity of AKAP12) and in a 108-Kb region of chromosome 10 (Log2FC = 2.81, P = 4.45 × 10−04, vicinity of FLRT2). These lncRNAs' impact on AKAP12 and FLRT2 is unlikely because these computationally predicted genes were not found expressed in macrophages (49). However, these lncRNAs' role in JD should not be excluded considering their potential activity on trans-target genes. Overall, there were 255 DE genes within the neighborhood of the lncRNA transcripts (Supplementary Figure 1).
Functions of lncRNA on the Different Biological Systems in Macrophages
The 255 DE genes were used to explore the functions of significant DE lncRNAs. Part of the Reactome pathways and Panther classification system analysis are presented in Table 3 with some of the Gene Ontology (GO) terms. A detailed report is found in Supplementary File 1. The analysis revealed 14 significant enrichment in the Reactome pathways, among the 1,287 enriched pathways. The highest significant enriched Reactome pathway was the Nuclear Receptor transcription pathway (R-HSA-383280, P = 1.69 × 10−5). Interestingly, part of this pathway is the Nuclear Receptor Subfamily 3 Group C Member 1 gene (NR3C1). NR3C1 was previously found expressed in bovine macrophages and being upregulated nearly 2-fold in response to MAP infection (Log2FC 0.8; P < 0.05) (49). This gene encodes a glucocorticoid receptor and chiefly binds small diffusible signaling molecules in the cytoplasm. Upon ligand binding, it migrates to the nucleus to bind glucocorticoid response elements in the promoters of glucocorticoid-responsive genes (68). While its role in JD is not reported, upon glucocorticoid-receptor bindings, various physiological functions, notably those involved in metabolism, inflammatory processes, and stress, are affected (69). This NR3C1 gene was also identified in the Panther classification system (Table 4 and Supplementary File 1). The NR3C1 molecular function falls in “glucocorticoid receptor activity” (GO:0004883; P = 8.45 × 10−03) of the biological activity “glucocorticoid mediated signaling pathway” (GO:0043402) with the highest fold (118.25) enrichment pathway (Table 3). As presented in Table 4, this nuclear receptor NR3C1 gene was found in “RNA polymerase II cis-regulatory region sequence-specific DNA binding” (GO:0000978; P = 1.16 × 10−04), in “cis-regulatory region sequence-specific DNA binding” (GO:0000987; P = 1.56 × 10−04), and in “sequence-specific DNA binding” (GO:0043565; P = 1.57 × 10−04) pathways (70–72).
Table 3. Enriched Reactome pathways using 255 DE lncRNA.
Table 4. Gene ontology results using 255 DE lncRNAs.
A second significant enriched pathway, which includes genes found expressed in bovine macrophages, is the RNA Polymerase I Promoter Escape (Table 3). Interestingly, the TWISTNB gene, which encodes the RNA polymerase I subunit F, was upregulated to 2.42 FC (P = 5 × 10−04) in macrophages from JD(+) cows (49), suggesting increased recruitment of Pol I to rDNA promoters in macrophages from JD(+) cows.
The complete list of GO terms from categories of the Panther classification system, notably molecular function, cellular components, and biological processes, is presented in Supplementary File 1. The biological process category reported 3,144 pathways, of which 2,117 had 2-fold enrichment or better; molecular function reported 148 pathways with 132 having at least a 2-fold enrichment; and 153 pathways had a cellular component with 26 pathways having a 2-fold enrichment or better (Supplementary File 1).
It is well-known for M. tuberculosis (73) and suggested for MAP (12) that preventing acidification or fusion of the phagosome to the lysosome is a surviving strategy. Interestingly, several lncRNAs were identified in cellular components associated with lytic vacuole membrane (GO:0098852), lysosomal membrane (GO:0005765), and endocytic vesicle membrane (GO:0030666). Two genes located in the lysosomal membrane were previously found downregulated in macrophages in response to MAP infection, namely, the Solute Carrier Family 2 Member 8 gene (SLC2A8) and the Solute Carrier Family 48 Member 1 (SLC48A1) (data shown). It is of particular interest that SLC48A1 encodes to a heme transporter in the context that MAP (70) as for other mycobacteria (71, 72) relies on the host for the acquisition of iron (Fe) which is critical for their growth. Other cellular components were the extracellular exosome (GO:0070062, P = 3.18 × 10−25) and vesicle (GO:1903561, P = 5.63 × 10−24) which were the most significant cellular components with a similar high fold enrichment (~30). Interestingly, the histamine N-methyltransferase gene (HNMT) was significantly downregulated by 2.7-fold in macrophages infected by MAP at 4 hpi (49).
One of the longest significant lncRNAs, as mentioned above, is in the vicinity of the FZD1-encoded transmembrane protein, identified associated with the cellular component focal adhesion (GO:0005925). The molecule FZD1 functions as protein binding (GO:0005515), a critical Wnt/β-catenin negative feedback loop for repression of Toll-like receptor (TLR)-triggered inflammatory responses (74). This cell receptor in macrophages performs Wnt signaling to regulate pro-inflammatory cytokines in response to bacteria and mycobacterial infection (67). Interestingly, FZD1 and its ligand Wnt3a are involved in reprogramming Mycobacterium tuberculosis-infected macrophages (75). This is particularly interesting for JD while supporting our hypothesis developed based on our previous study (49) and other studies (76) that JD(+) macrophages may be responsive because of tolerance, i.e., epigenetic reprogramming. FZD1 was found more expressed in macrophages from JD(+) cows than JD(–) macrophages and might explain the phenotypes observed for JD(+) macrophages (49).
lncRNA's Putative Role of JD-Associated lncRNA in the Immune Response
Though not among the most significant GO terms, immunoreceptor activity (GO:0140375) had a 5.18-fold enrichment. This high fold enrichment is of interest because this function is responsible for receiving a signal and transmitting it in a cell to initiate an immune response during an invasion by a pathogen. The top 10 GO terms involved with biological processes were associated with neutrophils, for instance, neutrophil degranulation (GO:0043312, P = 7.08 × 10−12), which is involved in regulated exocytosis of secretory granules, neutrophil-mediated immunity (GO:0002446, P = 3.2 × 10−09), and neutrophil activation (GO:0042119, P = 2.61 × 10−08). For Mtb infection, secreted products from neutrophils regulate the macrophage activity (77). Since neutrophils are part of the first line of innate immunity of healthy cows (78) and are the second type of cell migrating in lesions of experimentally MAP-infected calves (79), the involvement of these pathways to support neutrophils in their role for pathogen clearance is highly relevant. Interestingly, the “negative regulation of T cell-mediated cytotoxicity” (GO: GO:0001915, P = 3.1 × 10−04) was among the top enriched biological processes. Both genes associated with this pathway, notably CEACAM1 and interleukin (IL) 7 receptor gene (IL7R), were expressed in bovine macrophages, and most interestingly is IL7R that was up to 7-fold increased in macrophages in response to MAP infection at 4–8 hpi (49). Aberrant plasma IL7 and soluble IL7 receptor levels indicate impaired T-cell response in human tuberculosis (80), and gene expression of IL7R had significant discriminatory power between tuberculosis-positive and -negative patients (81). The relevance of studying the lncRNA associated with IL7R in the pathogenesis of MAP also merits that IL7R is an immune biomarker validated for detecting clinical tuberculosis (82).
Gene ontology (GO) indicated that lncRNA influences, among others, the inflammation, the evidence of which is well-documented for macrophages (83–85). The myeloid differentiation primary response (MyD88) pathway was confirmed to be affected by MAP infection in our previous study (49). Interestingly, previous studies have shown that mice with knocked-out Myd88 are highly susceptible to infection by M. tuberculosis (86), implying that the MyD88-dependent toll-like receptor signaling pathway could be a mycobacterial target for pathogen evasion of host responses. As expected, lncRNA target genes were significantly enriched for biological process GO terms involved in this immune regulation (GO:0034146, GO:0002755). Previous studies report that lncRNA may bind to their target genes; hence, it is unsurprising that the nucleic acid regulation GO terms were enriched (87). GO terms (e.g., GO:0043312, GO:0002446, and GO:0042119) were characterized by pathways associated with neutrophils reported to provide the first line of cellular defense against bacterial colonization in cattle (88).
We observed little overlap between predicted lncRNA candidates and the previously published cattle non-coding RNA. This little overlap may partially be because the lncRNA library is not fully curated in dairy cattle or that previously published non-coding RNA was identified in different tissues since lncRNAs show tissue- and cell-specific expressions. This study does have limitations; one is the incomplete bovine annotation where many unannotated genes exist, both protein- and non-protein-coding. In the current state, this would render the unknown transcripts in our data moot. Collaborative projects may eventually lead to creating a comprehensive map of functional elements in cattle's genome, allowing better identification of potential bovine lncRNA.
A good example is the Functional Annotation of Animal Genomes (FAANG) (89). Previous studies have identified lncRNA using in vitro infected macrophages (41). However, our study's novelty aimed to identify genome-wide lncRNAs using primary monocyte-derived macrophages for naturally MAP infected cows.
Our study's rigor also comes from the culture protocol of macrophages that do not bias macrophages' polarization. In Gupta et al. (41), macrophages were cultured in the presence of fetal bovine serum (FBS) and were differentiated in the presence of macrophage colony-stimulating factor (M-CSF). The effect of a culture medium on polarization when supporting monocytes' differentiation has been rigorously documented (90). It turns out that the presence of fetal bovine serum (FBS) and stimulating factors, such as GM-CSF or M-CSF allow a differentiation of the monocytes that activate and polarize them (90). The presence of FBS affects the functionality of monocyte-derived macrophages (91). FBS contains a considerable amount of immunoregulatory cytokines and bioactive molecules, notably the transforming growth factor-β (TGF-β) in notable concentrations (92) that yields a dominant immunosuppressive phenotype in the presence of M-CSF (93). This is all the more important since polarization impacts the lncRNA profile in macrophages (30). An additional novelty of our study is that macrophages were differentiated in the absence of FBS and that no stimulating factors were used to avoid polarization bias. It might explain that we identified 3,669 novel lncRNAs in bovine macrophages compared to 397 novel lncRNAs in the previous study (41).
Dairy and beef studies continuously use underlying biological information of point mutations (e.g., SNP) and genomic features to discover which genomic variants are theoretically enriched (94). The bovine lncRNAs reported here could further extend the underlying biological information by including a new class for genomic variants exclusive in lncRNA regions to enrich bovine health traits.
Although human data studies show that single-exon lncRNAs are more likely to be conserved (95), we did not discriminate against lncRNAs from abundant lowly expressed single-exonic fragments. Furthermore, following stringent filtering criteria based on other genomic features like ORF length, protein-coding potential, and expression levels, we identified 16,970 lncRNAs, 3,669 of which were novel lncRNAs with 255 potential cis target genes.
Conclusion
In summary, we have provided the lncRNA expression profile of macrophages from JD cows, together with many potential co-regulated candidate protein-coding genes. We have identified 13,589 lncRNAs, 3,669 of which are novel. Among those lncRNAs significantly upregulated in JD(+) macrophages, three were linked (<1 kb) to genes expressed in bovine macrophages, notably TWSG1, DPH1, and BOD1. Bacteria interfere with mammalian regulatory RNA expression, including lncRNA, to modify molecular and cellular signaling. We identified lncRNAs as having the potential to play a significant role in regulating several cellular activities, including the immune response. Although the mechanisms of MAP are currently unknown, we speculate that many of the identified lncRNAs are essential participants in the bovine innate immune response. At the same time, some could support MAP in macrophages in its duel to escape the surveillance of the immune system that lasts for years.
Data Availability Statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.
Ethics Statement
The animal study was reviewed and approved by Agriculture and Agri-Food Canada Animal Ethics Committee. Written informed consent was obtained from the owners for the participation of their animals in this study.
Author Contributions
NB: project administration, supervision, and funding acquisition. AM, OA, and NB: methodology. AM: data curation and formal analysis. AM and NB: writing. All authors provided the input in interpreting the results and have read and agreed to the published version of the manuscript.
Funding
This study was funded by Agri-Food and Agriculture Canada (Project AAFC J-000075 and AAFC J-000079). This research was supported by a contribution from the Dairy Research Cluster 3 (Dairy Farmers of Canada, Canadian Dairy Network and Agriculture, and Agri-Food Canada) under the Canadian Agricultural Partnership AgriScience Program (Project AAFC J-002095).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Acknowledgments
We appreciate the invaluable collaboration of the veterinarians and dairy producers in the Ontario and Quebec Provinces.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fvets.2021.639053/full#supplementary-material
Supplementary File 1. (A) Distribution of lncRNA per chromosome; (B) General Transfer Format (gtf) of identified lncRNA; (C) Reactome pathways; (D) General biological programs accomplished by the multiple molecular activities (Biological processes; GO terms); (E) Molecular-level activities performed by gene products (molecular functions; GO terms), and (F) The locations relative to cellular structures in which a gene product performs a function (Cellular component; GO terms).
Supplementary Figure 1. Illustration of gene network including genes within 100 kb of highly expressed lncRNA.
References
1. Whittington R, Donat K, Weber MF, Kelton D, Nielsen SS, Eisenberg S, et al. Control of paratuberculosis: who, why and how. A review of 48 countries. BMC Vet Res. (2019) 15:198. doi: 10.1186/s12917-019-1943-4
2. Whitlock RH, Buergelt C. Preclinical and clinical manifestations of paratuberculosis (including pathology). Vet Clin North Am Food Anim Pract. (1996) 12:345–56. doi: 10.1016/S0749-0720(15)30410-2
3. Smith RL, Strawderman RL, Schukken YH, Wells SJ, Pradhan AK, Espejo LA, et al. Effect of Johne's disease status on reproduction and culling in dairy cattle. J Dairy Sci. (2010) 93:3513–24. doi: 10.3168/jds.2009-2742
4. Garcia AB, Shalloo L. Invited review: the economic impact and control of paratuberculosis in cattle. J Dairy Sci. (2015) 98:5019–39. doi: 10.3168/jds.2014-9241
5. Corbett CS, Naqvi SA, Bauman CA, De Buck J, Orsel K, Uehlinger F, et al. Prevalence of Mycobacterium avium ssp. paratuberculosis infections in Canadian dairy herds. J Dairy Sci. (2018) 101:11218–28. doi: 10.3168/jds.2018-14854
6. Merkal RS, Larsen AB, Kopecky KE, Kluge JP, Monlux WS, Lehmann RP, et al. Experimental paratuberculosis in sheep after oral, intratracheal, or intravenous inoculation: serologic and intradermal tests. Am J Vet Res. (1968) 29:963–9.
7. Momotani E, Whipple DL, Thiermann AB, Cheville NF. Role of M cells and macrophages in the entrance of Mycobacterium paratuberculosis into Domes of Ileal Peyer's patches in calves. Vet Pathol. (1988) 25:131–7. doi: 10.1177/030098588802500205
8. Arsenault RJ, Maattanen P, Daigle J, Potter A, Griebel P, Napper S. From mouth to macrophage: mechanisms of innate immune subversion by Mycobacterium avium subsp. paratuberculosis. Vet Res. (2014) 45:54. doi: 10.1186/1297-9716-45-54
9. Koets AP, Eda S, Sreevatsan S. The within host dynamics of Mycobacterium avium ssp. paratuberculosis infection in cattle: where time and place matter. Vet Res. (2015) 46:61. doi: 10.1186/s13567-015-0185-0
10. Bannantine JP, Stabel JR. Killing of Mycobacterium avium subspecies paratuberculosis within macrophages. BMC Microbiol. (2002) 2:2. doi: 10.1186/1471-2180-2-2
11. Weiss DJ, Souza CD, Evanson OA, Sanders M, Rutherford M. Bovine monocyte TLR2 receptors differentially regulate the intracellular fate of Mycobacterium avium subsp. paratuberculosis and Mycobacterium avium subsp avium. J Leuk. Biol. (2008) 83:48–55. doi: 10.1189/jlb.0707490
12. Hussain T, Shah SZ, Zhao D, Sreevatsan S, Zhou X. The role of IL-10 in Mycobacterium avium subsp. paratuberculosis infection. Cell Commun Signal. (2016) 14:29. doi: 10.1186/s12964-016-0152-z
13. Pinedo PJ, Buergelt CD, Donovan GA, Melendez P, Morel L, Wu R, et al. Association between CARD15/NOD2 gene polymorphisms and paratuberculosis infection in cattle. Vet Microbiol. (2009) 134:346–52. doi: 10.1016/j.vetmic.2008.09.052
14. Ruiz-Larranaga O, Garrido JM, Iriondo M, Manzano C, Molina E, Koets AP, et al. Genetic association between bovine NOD2 polymorphisms and infection by Mycobacterium avium subsp. paratuberculosis in Holstein-Friesian cattle. Anim Genet. (2010) 41:652–5. doi: 10.1111/j.1365-2052.2010.02055.x
15. Pinedo PJ, Buergelt CD, Donovan GA, Melendez P, Morel L, Wu R, et al. Candidate gene polymorphisms (BoIFNG, TLR4, SLC11A1) as risk factors for paratuberculosis infection in cattle. Prev Vet Med. (2009) 91:189–96. doi: 10.1016/j.prevetmed.2009.05.020
16. Korou LM, Liandris E, Gazouli M, Ikonomopoulos J. Investigation of the association of the SLC11A1 gene with resistance/sensitivity of goats (Capra hircus) to paratuberculosis. Vet Microbiol. (2010) 144:353–8. doi: 10.1016/j.vetmic.2010.01.009
17. Ruiz-Larranaga O, Garrido JM, Manzano C, Iriondo M, Molina E, Gil A, et al. Identification of single nucleotide polymorphisms in the bovine solute carrier family 11 member 1 (SLC11A1) gene and their association with infection by Mycobacterium avium subspecies paratuberculosis. J Dairy Sci. (2010) 93:1713–21. doi: 10.3168/jds.2009-2438
18. Verschoor CP, Pant SD, You Q, Schenkel FS, Kelton DF, Karrow NA. Polymorphisms in the gene encoding bovine interleukin-10 receptor alpha are associated with Mycobacterium avium ssp. paratuberculosis infection status. BMC Genet. (2010) 11:23. doi: 10.1186/1471-2156-11-23
19. Cinar MU, Hizlisoy H, Akyüz B, Arslan K, Aksel EG, Gümüşsoy KS. Polymorphisms in toll-like receptor (TLR) 1, 4, 9 and SLC11A1 genes and their association with paratuberculosis susceptibility in Holstein and indigenous crossbred cattle in Turkey. J Genet. (2018) 97:1147–54. doi: 10.1007/s12041-018-1008-7
20. Juste RA, Vazquez P, Ruiz-Larranaga O, Iriondo M, Manzano C, Agirre M, et al. Association between combinations of genetic polymorphisms and epidemiopathogenic forms of bovine paratuberculosis. Heliyon. (2018) 4:e00535. doi: 10.1016/j.heliyon.2018.e00535
21. Clarke CJ. The pathology and pathogenesis of paratuberculosis in ruminants and other species. J Compar Pathol. (1997) 116:217–61. doi: 10.1016/S0021-9975(97)80001-1
22. Sweeney RW. Pathogenesis of paratuberculosis. Vet Clin North Am Food Anim Pract. (2011) 27:537–46. doi: 10.1016/j.cvfa.2011.07.001
23. Kirkpatrick BW. Genetics of Host Susceptibility to Paratuberculosis. Cambridge, MA: CABI (2010).
24. Stevenson K. Genetic diversity of Mycobacterium avium subspecies paratuberculosis and the influence of strain type on infection and pathogenesis: a review. Vet Res. (2015) 46:64. doi: 10.1186/s13567-015-0203-2
25. Stabel JR. Host responses to Mycobacterium avium subsp. paratuberculosis: a complex arsenal. Anim Health Res Rev. (2006) 7:61–70. doi: 10.1017/S1466252307001168
26. Duval M, Cossart P, Lebreton A. Mammalian microRNAs and long noncoding RNAs in the host-bacterial pathogen crosstalk. Semin Cell Dev Biol. (2017) 65:11–9. doi: 10.1016/j.semcdb.2016.06.016
27. Zur Bruegge J, Einspanier R, Sharbati S. A long journey ahead: long non-coding RNAs in bacterial infections. Front Cell Infect Microbiol. (2017) 7:95. doi: 10.3389/fcimb.2017.00095
28. Fitzgerald KA, Caffrey DR. Long noncoding RNAs in innate and adaptive immunity. Curr Opin Immunol. (2014) 26:140–6. doi: 10.1016/j.coi.2013.12.001
29. Elling R, Chan J, Fitzgerald KA. Emerging role of long noncoding RNAs as regulators of innate immune cell development and inflammatory gene expression. Eur J Immunol. (2016) 46:504–12. doi: 10.1002/eji.201444558
30. Ahmad I, Valverde A, Ahmad F, Naqvi AR. Long noncoding RNA in myeloid and lymphoid cell differentiation, polarization and function. Cells. (2020) 9:269. doi: 10.3390/cells9020269
31. Kim J, Kim KM, Noh JH, Yoon JH, Abdelmohsen K, Gorospe M. Long noncoding RNAs in diseases of aging. Biochim Biophys Acta. (2016) 1859:209–21. doi: 10.1016/j.bbagrm.2015.06.013
32. Li L, Song X. The working modules of long noncoding RNAs in cancer cells. Adv Exp Med Biol. (2016) 927:49–67. doi: 10.1007/978-981-10-1498-7_2
33. Talebian S, Daghagh H, Yousefi B, Ozkul Y, Ilkhani K, Seif F, et al. The role of epigenetics and non-coding RNAs in autophagy: a new perspective for thorough understanding. Mech Ageing Dev. (2020) 190:111309. doi: 10.1016/j.mad.2020.111309
34. Zacharopoulou E, Gazouli M, Tzouvala M, Vezakis A, Karamanolis G. The contribution of long non-coding RNAs in inflammatory bowel diseases. Dig Liver Dis. (2017) 49:1067–72. doi: 10.1016/j.dld.2017.08.003
35. Ge Q, Dong Y, Lin G, Cao Y. Long noncoding RNA antisense noncoding RNA in the INK4 locus correlates with risk, severity, inflammation and infliximab efficacy in Crohn's disease. Am J Med Sci. (2019) 357:134–42. doi: 10.1016/j.amjms.2018.10.016
36. Weikard R, Hadlich F, Kuehn C. Identification of novel transcripts and noncoding RNAs in bovine skin by deep next generation sequencing. BMC Genomics. (2013) 14:789. doi: 10.1186/1471-2164-14-789
37. Billerey C, Boussaha M, Esquerre D, Rebours E, Djari A, Meersseman C, et al. Identification of large intergenic non-coding RNAs in bovine muscle using next-generation transcriptomic sequencing. BMC Genomics. (2014) 15:499. doi: 10.1186/1471-2164-15-499
38. Koufariotis LT, Chen YP, Chamberlain A, Vander Jagt C, Hayes BJ. A catalogue of novel bovine long noncoding RNA across 18 tissues. PLoS ONE. (2015) 10:e0141225. doi: 10.1371/journal.pone.0141225
39. Tong C, Chen Q, Zhao L, Ma J, Ibeagha-Awemu EM, Zhao X. Identification and characterization of long intergenic noncoding RNAs in bovine mammary glands. BMC Genomics. (2017) 18:468. doi: 10.1186/s12864-017-3858-4
40. Ibeagha-Awemu EM, Li R, Dudemaine PL, Do DN, Bissonnette N. Transcriptome analysis of long non-coding RNA in the bovine mammary gland following dietary supplementation with linseed oil and safflower oil. Int J Mol Sci. (2018) 19:3610. doi: 10.3390/ijms19113610
41. Gupta P, Peter S, Jung M, Lewin A, Hemmrich-Stanisak G, Franke A, et al. Analysis of long non-coding RNA and mRNA expression in bovine macrophages brings up novel aspects of Mycobacterium avium subspecies paratuberculosis infections. Sci Rep. (2019) 9:1571. doi: 10.1038/s41598-018-38141-x
42. Cabili MN, Trapnell C, Goff L, Koziol M, Tazon-Vega B, Regev A, et al. Integrative annotation of human large intergenic noncoding RNAs reveals global properties and specific subclasses. Genes Dev. (2011) 25:1915–27. doi: 10.1101/gad.17446611
43. Chen G, Yin K, Shi L, Fang Y, Qi Y, Li P, et al. Comparative analysis of human protein-coding and noncoding RNAs between brain and 10 mixed cell lines by RNA-Seq. PLoS ONE. (2011) 6:e28318. doi: 10.1371/journal.pone.0028318
44. Derrien T, Johnson R, Bussotti G, Tanzer A, Djebali S, Tilgner H, et al. The GENCODE v7 catalog of human long noncoding RNAs: analysis of their gene structure, evolution, and expression. Genome Res. (2012) 22:1775–89. doi: 10.1101/gr.132159.111
45. Wang L, Park HJ, Dasari S, Wang S, Kocher JP, Li W. CPAT: coding-potential assessment tool using an alignment-free logistic regression model. Nucleic Acids Res. (2013) 41:e74. doi: 10.1093/nar/gkt006
46. Camargo AP, Sourkov V, Pereira GAG, Carazzolle MF. RNAsamba: neural network-based assessment of the protein-coding potential of RNA sequences. NAR Genom Bioinform. (2020) 2:lqz024. doi: 10.1093/nargab/lqz024
47. Fock-Chow-Tho D, Topp E, Ibeagha-Awemu EA, Bissonnette N. Comparison of commercial DNA extraction kits and quantitative PCR systems for better sensitivity in detecting the causative agent of paratuberculosis in dairy cow fecal samples. J Dairy Sci. (2017) 100:572–81. doi: 10.3168/jds.2016-11384
48. Collins MT. Interpretation of a commercial bovine paratuberculosis enzyme-linked immunosorbent assay by using likelihood ratios. Clin Diagn Lab Immunol. (2002) 9:1367–71. doi: 10.1128/CDLI.9.6.1367-1371.2002
49. Ariel O, Gendron D, Dudemaine PL, Gevry N, Ibeagha-Awemu EM, Bissonnette N. Transcriptome profiling of bovine macrophages infected by Mycobacterium avium spp. paratuberculosis depicts foam cell and innate immune tolerance phenotypes. Front Immunol. (2019) 10:2874. doi: 10.3389/fimmu.2019.02874
50. Edgar RC. Search and clustering orders of magnitude faster than BLAST. Bioinformatics. (2010) 26:2460–1. doi: 10.1093/bioinformatics/btq461
51. Fang S, Zhang L, Guo J, Niu Y, Wu Y, Li H, et al. NONCODEV5: a comprehensive annotation database for long non-coding RNAs. Nucleic Acids Res. (2018) 46:D308–14. doi: 10.1093/nar/gkx1107
52. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. (2014) 15:550. doi: 10.1186/s13059-014-0550-8
53. Mi H, Muruganujan A, Casagrande JT, Thomas PD. Large-scale gene function analysis with the PANTHER classification system. Nat Protoc. (2013) 8:1551–66. doi: 10.1038/nprot.2013.092
54. Jassal B, Matthews L, Viteri G, Gong C, Lorente P, Fabregat A, et al. The reactome pathway knowledgebase. Nucleic Acids Res. (2020) 48:D498–503. doi: 10.1093/nar/gkz1031
55. Yang X, Yang J, Wang J, Wen Q, Wang H, He J, et al. Microarray analysis of long noncoding RNA and mRNA expression profiles in human macrophages infected with Mycobacterium tuberculosis. Sci Rep. (2016) 6:38963. doi: 10.1038/srep38963
56. Pertea M, Kim D, Pertea GM, Leek JT, Salzberg SL. Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat Protoc. (2016) 11:1650–67. doi: 10.1038/nprot.2016.095
57. Weikard R, Demasius W, Kuehn C. Mining long noncoding RNA in livestock. Anim Genet. (2017) 48:3–18. doi: 10.1111/age.12493
58. Kern C, Wang Y, Chitwood J, Korf I, Delany M, Cheng H, et al. Genome-wide identification of tissue-specific long non-coding RNA in three farm animal species. BMC Genomics. (2018) 19:684. doi: 10.1186/s12864-018-5037-7
59. Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. (2013) 14:R36. doi: 10.1186/gb-2013-14-4-r36
60. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. (2015) 12:357–60. doi: 10.1038/nmeth.3317
61. Baruzzo G, Hayer KE, Kim EJ, Di Camillo B, Fitzgerald GA, Grant GR. Simulation-based comprehensive benchmarking of RNA-seq aligners. Nat Methods. (2017) 14:135–9. doi: 10.1038/nmeth.4106
62. Ponjavic J, Oliver PL, Lunter G, Ponting CP. Genomic and transcriptional co-localization of protein-coding and long non-coding RNA pairs in the developing brain. PLoS Genet. (2009) 5:e1000617. doi: 10.1371/journal.pgen.1000617
63. Kotake Y, Nakagawa T, Kitagawa K, Suzuki S, Liu N, Kitagawa M, et al. Long non-coding RNA ANRIL is required for the PRC2 recruitment to and silencing of p15(INK4B) tumor suppressor gene. Oncogene. (2011) 30:1956–62. doi: 10.1038/onc.2010.568
64. Teng M, Love MI, Davis CA, Djebali S, Dobin A, Graveley BR, et al. A benchmark for RNA-seq quantification pipelines. Genome Biol. (2016) 17:74. doi: 10.1186/s13059-016-1060-7
65. Williams CR, Baccarella A, Parrish JZ, Kim CC. Empirical assessment of analysis workflows for differential expression analysis of human samples using RNA-Seq. BMC Bioinformatics. (2017) 18:38. doi: 10.1186/s12859-016-1457-z
66. Merino GA, Conesa A, Fernandez EA. A benchmarking of workflows for detecting differential splicing and differential expression at isoform level in human RNA-seq studies. Brief Bioinform. (2019) 20:471–81. doi: 10.1093/bib/bbx122
67. Schaale K, Neumann J, Schneider D, Ehlers S, Reiling N. Wnt signaling in macrophages: augmenting and inhibiting mycobacteria-induced inflammatory responses. Eur J Cell Biol. (2011) 90:553–9. doi: 10.1016/j.ejcb.2010.11.004
68. Bray PJ, Cotton RG. Variations of the human glucocorticoid receptor gene (NR3C1): pathological and in vitro mutations and polymorphisms. Hum Mutat. (2003) 21:557–68. doi: 10.1002/humu.10213
69. Vitellius G, Trabado S, Bouligand J, Delemer B, Lombes M. Pathophysiology of glucocorticoid signaling. Ann Endocrinol. (2018) 79:98–106. doi: 10.1016/j.ando.2018.03.001
71. Olakanmi O, Kesavalu B, Abdalla MY, Britigan BE. Iron acquisition by Mycobacterium tuberculosis residing within myeloid dendritic cells. Microb Pathog. (2013) 65:21–8. doi: 10.1016/j.micpath.2013.09.002
72. Chao A, Sieminski PJ, Owens CP, Goulding CW. Iron acquisition in Mycobacterium tuberculosis. Chem Rev. (2019) 119:1193–220. doi: 10.1021/acs.chemrev.8b00285
73. Harriff MJ, Purdy GE, Lewinsohn DM. Escape from the phagosome: the explanation for MHC-i processing of mycobacterial antigens? Front Immunol. (2012) 3:40. doi: 10.3389/fimmu.2012.00040
74. Li Y, Shi J, Yang J, Ma Y, Cheng L, Zeng J, et al. A Wnt/beta-catenin negative feedback loop represses TLR-triggered inflammatory responses in alveolar epithelial cells. Mol Immunol. (2014) 59:128–35. doi: 10.1016/j.molimm.2014.02.002
75. Neumann J, Schaale K, Farhat K, Endermann T, Ulmer AJ, Ehlers S, et al. Frizzled1 is a marker of inflammatory macrophages, and its ligand Wnt3a is involved in reprogramming Mycobacterium tuberculosis-infected macrophages. FASEB J. (2010) 24:4599–612. doi: 10.1096/fj.10-160994
76. Khare S, Lawhon SD, Drake KL, Nunes JE, Figueiredo JF, Rossetti CA, et al. Systems biology analysis of gene expression during in vivo Mycobacterium avium paratuberculosis enteric colonization reveals role for immune tolerance. PLoS ONE. (2012) 7:e42127. doi: 10.1371/journal.pone.0042127
77. Scordo JM, Arcos J, Kelley HV, Diangelo L, Sasindran SJ, Youngmin E, et al. Mycobacterium tuberculosis cell wall fragments released upon bacterial contact with the human lung mucosa alter the neutrophil response to infection. Front Immunol. (2017) 8:307. doi: 10.3389/fimmu.2017.00307
78. Alhussien MN, Dang AK. Potential roles of neutrophils in maintaining the health and productivity of dairy cows during various physiological and physiopathological conditions: a review. Immunol Res. (2019) 67:21–38. doi: 10.1007/s12026-019-9064-5
79. Khare S, Nunes JS, Figueiredo JF, Lawhon SD, Rossetti CA, Gull T, et al. Early phase morphological lesions and transcriptional responses of bovine ileum infected with Mycobacterium avium subsp. paratuberculosis. Vet Pathol. (2009) 46:717–28. doi: 10.1354/vp.08-VP-0187-G-FL
80. Lundtoft C, Afum-Adjei Awuah A, Rimpler J, Harling K, Nausch N, Kohns M, et al. Aberrant plasma IL-7 and soluble IL-7 receptor levels indicate impaired T-cell response to IL-7 in human tuberculosis. PLoS Pathog. (2017) 13:e1006425. doi: 10.1371/journal.ppat.1006425
81. Mihret A, Loxton AG, Bekele Y, Kaufmann SH, Kidd M, Haks MC, et al. Combination of gene expression patterns in whole blood discriminate between tuberculosis infection states. BMC Infect Dis. (2014) 14:257. doi: 10.1186/1471-2334-14-257
82. Jenum S, Dhanasekaran S, Lodha R, Mukherjee A, Kumar Saini D, Singh S, et al. Approaching a diagnostic point-of-care test for pediatric tuberculosis through evaluation of immune biomarkers across the clinical disease spectrum. Sci Rep. (2016) 6:18520. doi: 10.1038/srep18520
83. Carpenter S, Aiello D, Atianand MK, Ricci EP, Gandhi P, Hall LL, et al. A long noncoding RNA mediates both activation and repression of immune response genes. Science. (2013) 341:789–92. doi: 10.1126/science.1240925
84. Hu G, Gong AY, Wang Y, Ma S, Chen X, Chen J, et al. LincRNA-Cox2 promotes late inflammatory gene transcription in macrophages through modulating SWI/SNF-mediated chromatin remodeling. J Immunol. (2016) 196:2799–808. doi: 10.4049/jimmunol.1502146
85. Covarrubias S, Robinson EK, Shapleigh B, Vollmers A, Katzman S, Hanley N, et al. CRISPR/Cas-based screening of long non-coding RNAs (lncRNAs) in macrophages with an NF-kappaB reporter. J Biol Chem. (2017) 292:20911–20. doi: 10.1074/jbc.M117.799155
86. Arsenault RJ, Li Y, Maattanen P, Scruten E, Doig K, Potter A, et al. Altered Toll-like receptor 9 signaling in Mycobacterium avium subsp. paratuberculosis-infected bovine monocytes reveals potential therapeutic targets. Infect Immun. (2013) 81:226–37. doi: 10.1128/IAI.00785-12
87. Flintoft L. Non-coding RNA: structure and function for lncRNAs. Nat Rev Genet. (2013) 14:598. doi: 10.1038/nrg3561
88. Hammon DS, Evjen IM, Dhiman TR, Goff JP, Walters JL. Neutrophil function and energy status in Holstein cows with uterine health disorders. Vet Immunol Immunopathol. (2006) 113:21–9. doi: 10.1016/j.vetimm.2006.03.022
89. Andersson L, Archibald AL, Bottema CD, Brauning R, Burgess SC, Burt DW, et al. Coordinated international action to accelerate genome-to-phenome with FAANG, the functional annotation of animal genomes project. Genome Biol. (2015) 16:57. doi: 10.1186/s13059-015-0622-4
90. Emam M, Tabatabaei S, Sargolzaei M, Sharif S, Schenkel F, Mallard B. The effect of host genetics on in vitro performance of bovine monocyte-derived macrophages. J Dairy Sci. (2019) 102:9107–16. doi: 10.3168/jds.2018-15960
91. Rey-Giraud F, Hafner M, Ries CH. In vitro generation of monocyte-derived macrophages under serum-free conditions improves their tumor promoting functions. PLoS ONE. (2012) 7:e42656. doi: 10.1371/journal.pone.0042656
92. Oida T, Weiner HL. Depletion of TGF-beta from fetal bovine serum. J Immunol Methods. (2010) 362:195–8. doi: 10.1016/j.jim.2010.09.008
93. Mia S, Warnecke A, Zhang XM, Malmstrom V, Harris RA. An optimized protocol for human M2 macrophages using M-CSF and IL-4/IL-10/TGF-beta yields a dominant immunosuppressive phenotype. Scand J Immunol. (2014) 79:305–14. doi: 10.1111/sji.12162
94. Koufariotis L, Chen YP, Bolormaa S, Hayes BJ. Regulatory and coding genome regions are enriched for trait associated variants in dairy and beef cattle. BMC Genomics. (2014) 15:436. doi: 10.1186/1471-2164-15-436
Keywords: bovine, genomics, long non-coding RNA, macrophages, paratuberculosis (Mycobacterium avium ssp. paratuberculosis), Johne's disease, MAP disease
Citation: Marete A, Ariel O, Ibeagha-Awemu E and Bissonnette N (2021) Identification of Long Non-coding RNA Isolated From Naturally Infected Macrophages and Associated With Bovine Johne's Disease in Canadian Holstein Using a Combination of Neural Networks and Logistic Regression. Front. Vet. Sci. 8:639053. doi: 10.3389/fvets.2021.639053
Received: 08 December 2020; Accepted: 15 February 2021;
Published: 22 April 2021.
Edited by:
Miguel Salgado, Austral University of Chile, ChileReviewed by:
Marta Alonso-Hearn, Animalien Osasuna, NEIKER-Instituto Vasco de Investigación y Desarrollo Agrario, SpainMohanned Naif Alhussien, National Dairy Research Institute (ICAR), India
Copyright © 2021 Marete, Ariel, Ibeagha-Awemu and Bissonnette. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Nathalie Bissonnette, nathalie.bissonnette@canada.ca
†ORCID: Andrew Marete orcid.org/0000-0003-2301-4168