- Open Access
Analysis of long non-coding RNA expression profile of bovine monocyte-macrophage infected by Mycobacterium avium subsp. paratuberculosis
BMC Genomics volume 23, Article number: 768 (2022)
Mycobacterium avium subsp. paratuberculosis (MAP) is the causative agent of paratuberculosis. As a potential zoonotic pathogen, MAP also seriously threatens human health and social security. At present, long non-coding RNA (lncRNA) has attracted wide attention as an useful biomarker in various diseases. Therefore, our study analyzed the lncRNA expression profiles and lncRNA-mRNA regulatory network of MAP infected bovine monocytes-macrophages and uninfected bovine cells by high-throughput sequencing. A total of 4641 differentially expressed lncRNAs genes were identified, including 3111 up-regulated genes and 1530 down-regulated genes. In addition, lncRNA-mRNA interaction analysis was performed to predict the target genes of lncRNA. Among them, after MAP infection, 86 lncRNAs targeted to mRNA, of which only 6 genes were significantly different. The results of Gene Ontology (GO) enrichment analysis showed that the differentially expressed genes significantly enriched in functional groups were related to immune regulation. Multiple signal pathways including NF-κB, NOD-like receptor, Cytokine-cytokine receptor, Toll-like receptor signaling pathway, Chemokine signaling pathway, and other important biochemical, metabolic and signal transduction pathways were enriched in Kyoto Encyclopedia of Genes and Genomes (KEGG). In this study, analysis of macrophage transcriptomes in response to MAP infection is expected to provide key information to deeply understand role of the pathogen in initiating an inappropriate and persistent infection in susceptible hosts and molecular mechanisms that might underlie the early phases of paratuberculosis.
Paratuberculosis, also known as Johne’s disease (JD), can cause chronic enteritis in cattle and other ruminants. This disease is endemic in many ruminant herds worldwide  and has caused significant economic losses to the dairy industry. In a survey of 48 countries, paratuberculosis was found to be very common in livestock. In about half of the countries, more than 20% herds are infected with paratuberculosis . Paratuberculosis has been found of more than 68% of cattle herds in the United States . In Canada, the economic losses associated with JD in dairy herds were estimated to be more than US$35–57/year/cow . The disease may cause profuse chronic diarrhea, submandibular edema, weight loss despite normal appetite, malnutrition, lethargy, and even death . Mycobacterium avium subsp. paratuberculosis (MAP), pathogen of paratuberculosis, is a slowly growing fastidious acid-fast bacterium belonging to the Mycobacterium avium complex (MAC). MAP can infect a variety of wild ruminants, as well as common mammals including cattle and sheep, where cattle are the most common host for MAP infection [6,7,8]. The most common symptoms of this disease in cattle are loss of milk production, weight loss and diarrhoea, whereas in sheep and goats, the symptoms are emaciation, anorexia and severe disability . In addition to cattle, sheep and goats, MAP has been found in ruminants such as white-tailed deer , red deer , elk  and non-ruminants such as rabbits and foxes , and even primates such as mandrills and macaques , indicating the wide adaptability of this mycobacterium to host species. This extensive host species may lead to the wide spread of MAP, which makes it difficult to control JD on a global scale. Due to the widespread and impact of JD in global livestock and the huge economic losses it cause, it has been recognized by the World Organization for Animal Health as one of the most important infectious diseases. In addition, MAP may be associated with the development of human Crohn’s disease, type I diabetes, multiple sclerosis, or rheumatoid arthritis [15,16,17,18].
The main reasons why MAP is prevalent and difficult to eradicate are the amazing speed of transmission and the widespread and susceptibility of the host. The pathogen can be transmitted by fecal–oral route, feces and milk of subclinical or clinically infected animals, allowing MAP to spread to susceptible calves, the environment, and dairy products for human consumption [19, 20]. Milk containing MAP is particularly worrying because there is increasing evidence that MAP is a zoonotic pathogen . MAP is an obligate intracellular pathogen that can multiply in host’s macrophages. A mouse model showed that the junction between MAP and intestinal mucosa and the translocation of MAP through intestinal mucosa were mediated by M cells and intestinal cells after MAP was ingested . In the submucosa, MAP was uptaken by macrophages that play a critical role in the host–pathogen interaction of JD. They can not only mediate the destructive effect of MAP, but when macrophages are attacked and subverted by MAP, they can be transformed into safe havens for pathogens to survive, proliferate and spread . In different periods of infection, MAP promotes or inhibits apoptosis respectively .
Long non-coding RNA (lncRNA) is a type of non-coding RNA with a length greater than 200 nt, which is mainly transcribed from the antisense strand and spacer region of protein-coding genes. LncRNA has been shown to participate in many important regulatory processes such as X chromosome silencing, genome imprinting, chromatin modification, transcription activation, transcription interference, and nuclear transport . Besides, as a competitive endogenous RNA (ceRNA), lncRNA can affect the pathophysiological environment of different organisms and ultimately affect gene expression . The lncRNA-mediated ceRNA network affects the ischemic stroke (IS) process, especially the immune response after stroke . LncRNA has been reported to be involved in a variety of diseases and is considered as molecular biomarkers for disease diagnosis and treatment [28,29,30,31]. LncRNAs are also regarded as critical regulators of gene expression at both the transcriptional and post-transcriptional levels. The lncRNA located in nucleus participates in gene regulation at epigenetic and transcription level, for example, DNA methylation . At the same time, lncRNA located in the cytoplasm participates in gene regulation at post-transcriptional and translational levels, such as regulation of mRNA metabolism . The expression level of lncRNA is usually lower than that of mRNA, but lncRNA has stronger tissue specificity , indicating an indispensable role in the process of cell type specificity . So far, there are few studies on the lncRNA expression profile of MAP-infected bovine peripheral blood mononuclear macrophages. Therefore, compared with studies in human diseases and mouse disease models, there are still many unsolved issues of lncRNA’s immune regulation and genetic changes in response to bacteria and other external infections in cattle.
The potential zoonotic threat of MAP and the huge economic losses caused by MAP prompt the formulation of effective management strategies to prevent the spread of the disease. Therefore, understanding the pathogenic mechanism of MAP and its transmission mode is crucial for vaccine development and better control of the disease progression. It is well known that macrophages are the target cells where MAP is able to survive and multiply and play a crucial role in the host–pathogen interaction. Currently, macrophages have been used as pathogenic models to study the pro- or anti-inflammatory transition of immune cell effector function following MAP exposure due to its crucial role in the development and control of MAP disease. To better study the pathogenic mechanism and immune methods of MAP in cattle, high-throughput sequencing was performed to analyze the lncRNA expression profile and lncRNA-mRNA interaction network of MAP infected bovine monocyte-derived macrophages (MDM).
Cell total RNA extraction results
The total RNA concentration and quality test results are shown in Supplementary Table 1 and Supplementary Fig. 1. The RNA quality is good and meets the requirements for continuing to construct a library and sequencing.
Analysis of the overall expression level of lncRNA
The expression level of lncRNA mainly uses FPKM (Fragments Per Kilobase of exon model per Million mapped reads) to measure the abundance value of gene expression, and the expression abundance of known genes in different samples is counted by FPKM. FPKM represents the number of sequenced fragments contained in every thousand transcript sequenced bases per million sequenced bases. Simply put, the FPKM value can be understood as the expression level of lncRNA. The five types of class code definition information are shown in Supplementary Table 2. All lncRNA expression profile information is shown in Supplementary Table 3. We made statistics on the percentage of lncRNA different class code in each sample, and the pie chart is shown in Fig. 1. The score statistical box diagram of lncRNA CPC (Coding Potential Calculator) and CNCI (Coding-Non-Coding Index) in each sample is shown in Fig. 2A and B. Among the 4928 novel non-coding transcripts predicted by CPC and CNCI, most novel lncRNAs were defined on chromosomes 3 and 10, with 241 and 285 lncRNAs, respectively.
Visualization results of lncRNA genome in different samples
In order to show the distribution of lncRNA candidates in the chromosomes more intentionally, we used circos software to perform genome mapping on the lncRNAs obtained from the screening. It is mainly divided into two parts for display, one is for genome mapping according to different classifications of lncRNA, and the other is for genome mapping according to lncRNA in different samples. When mapping, every chromosome is based on every 25 mb. When mapping lncRNA genome visualization in different samples, the expression quantity of lncRNA in each segment is counted and plotted, as shown in Fig. 3A. The number of lncRNAs in each segment should be counted and plotted when analyzing the genome visualization of different types of lncRNAs, as shown in Fig. 3B.
Analysis of differential lncRNA expression results
The experimental results showed that there were a total of 4641 differentially expressed lncRNA genes (Supplementary Table 4), including 3111 up-regulated genes and 1530 down-regulated genes. There were 102 genes with significant differences among the up-regulated genes, and 91 genes with significant differences among the down-regulated genes. In order to show the overall distribution of differentially expressed genes more intuitively, volcanic maps (Fig. 4) and heat maps (Fig. 5) were drawn for all genes in differential expression analysis.
Comparison of lncRNA and mRNA structural characteristics
LncRNA is also called long-chain non-coding RNA, with a length greater than 200 nt, which is mainly transcribed from the antisense strand and spacer region of protein-coding genes. They regulate the expression level of genes in the form of RNA at multiple levels (epigenetic regulation, transcription regulation, post-transcriptional regulation, etc.). The biological content is quite abundant, accounting for about 4–8% of RNA (mRNA accounts for about 1–2%). Therefore, the follow-up functional annotation of lncRNA can be made by analyzing the length and expression correlation of lncRNA and mRNA (Supplementary Table 5 and Fig. 6). It has been reported that lncRNA has no obvious ORF and has very low protein coding potential [36, 37]. We compare the length of the ORF of lncRNA and mRNA translation and discuss the differences between the two (Supplementary Tables 6 and 7). The principle of ORF analysis is mainly based on the six-frame translation principle of nucleic acid (Fig. 7A and B). In addition, the number of lncRNA and mRNA exons (Supplementary Table 8 and Fig. 8) and the expression level of lncRNA and mRNA (Supplementary Table 9 and Fig. 9) were also counted.
lncRNA and mRNA interaction analysis
lncRNA target gene prediction
There are two main types of lncRNA regulation. One is cis-regulation, that is, lncRNA regulates the expression of its neighboring genes. The cis-regulation of target genes by lncRNA is mainly based on the prediction of the positional relationship, which defines the differentially expressed lncRNA and differentially expressed mRNA in the upper and lower reaches of the chromosome in the range of 100 kbp, and this type of lncRNA constitutes cis-regulation [38, 39]. Specifically, cis-regulation mainly relies on cis-acting elements. Cis-acting elements refer to sequences that can affect gene expression in the flanking sequences of genes, including promoters, enhancers, regulatory sequences and inducible elements, etc. The role is to participate in the regulation of gene expression in the nucleus, and cis-acting elements are usually transcribed into non-coding RNA; the second type of regulation is trans regulation, that is, lncRNA regulates the expression of genes across chromosomes. The trans regulation of target genes is mainly based on the amount of free energy required to form a secondary structure between lncRNA and mRNA sequences. If the binding of two sequences requires very low free energy, there may be an interaction between them. In this article, we mainly conduct predictive analysis of differential mRNA and differential lncRNA regulated by cis. The schematic diagram of lncRNA regulation mechanism is shown in Fig. 10.
In this paper, only the upstream and downstream 100 kbp mRNA and lncRNA are predicted in cis-regulation. A total of 86 lncRNAs were targeted to mRNA in MAP-infected and control MDM cells (Supplementary Table 10). Among them, only 6 targets have significant differences, see Table 1. Among them, the differential lncRNAs are located on the 2nd, 5th, 6th, 11th, 14th and 19th chromosomes, respectively, and the corresponding target mRNAs are CYP1B1, IMPA1, FXR2, ECE1, BHLHE41 and IL-8.
Differential gene GO and KEGG enrichment analysis
After infecting bovine monocyte-macrophage cells with MAP, in order to understand the specific functional roles of the differentially expressed genes, GO and KEGG enrichment analysis were performed on the differentially expressed genes.
GO has a total of three ontology, which describe the molecular function, cellular component, and biological process of the gene. The statistical description of these three ontologies in the text is shown in Supplementary Table 11. Among all the selected lncRNA differential genes, we divided the GO annotations corresponding to these differential genes into three categories according to Biological Process, Cellar Component and Molecular Function, and classified the GO functions in each category according to the target genes annotated. The numbers are sorted from high to low, and the top 25, 15, and 10 GO categories are selected for mapping (Fig. 11). In the classification of Biological Process, the most significant group was the regulation of single stranded viral RNA replication via double stranded DNA intermediate (GO: 0,045,091) GO functional group, followed by the trabecular meshwork development (GO: 0,002,930) GO functional group, with respectively 1 annotated gene IL-8 and 1 annotated gene CYP1B1. In the Molecular Function classification, the most important GO functional groups are the lithium ion binding functional group (GO: 0,031,403) with 1 annotated gene IMPA1 and the interleukin-8 receptor binding (GO: 0,005,153) functional group with 1 annotated gene IL-8. In the classification of Cellar Component, the most prominent functional group was the polysome (GO: 0,005,844) group with 1 annotated gene FXR2, followed by the ruffle (GO: 0,001,726) function group with 1 annotated gene TLN1.
The basic unit of GO is term, and each term corresponds to an attribute. The GO functional significance enrichment analysis first maps all the significantly differentially expressed genes to each term of the Gene Ontology database, calculates the number of genes in each term, and then applies the hypergeometric test to find out the significant difference compared with the whole genome background. Significantly enriched GO entries among differentially expressed genes. In this experiment, ggplot2 was used for GO enrichment analysis and Top20 GO_term was plotted, and the results were displayed in a scatter plot (Fig. 12).
In organisms, different genes coordinately perform their biological functions. Pathway-based analysis helps to further understand the biological functions of genes. KEGG is the main public pathway database. The results of KEGG enrichment analysis in this study are shown in Supplementary Table 12. Use ggplot2 to perform enrichment analysis on KEGG, and select the top 20 pathways according to the enriched P value to display the results in a scatter plot (Fig. 13). The KEGG pathway map shows that the enriched target genes are involved in multiple important biochemical, metabolic and signal transduction pathways, such as NOD-like receptor signaling pathway (ko04621), NF-κB signaling pathway (ko04064), Toll-like receptor signaling pathway (ko04620), Chemokine signaling pathway (ko04062), and Cytokine-cytokine receptor interaction (ko04060).
MAP can cause paratuberculosis in a variety of domestic and wild animals, which has resulted in serious economic losses and threats to human health and public health security. The disease is chronic, consumptive, and incurable. MAP is mainly transmitted through the fecal–oral route. Besides, it can also be spread via saliva, uterine fluid or semen [43, 44]. Another study has shown that MAP can also be transmitted through aerosols . There is also a hazard of lateral transmission between herds.
Compared with the research progress of lncRNA in other species such as humans and mice, there are few studies on the mechanism of lncRNA regulating bovine immune response. Recently, with the application of high-throughput sequencing technology, studies on the identification of lncRNA expression in cattle immune-related tissues, cells and differences in expression during the period of bacterial/virus infection are gradually increasing. High-throughput sequencing technology has greatly enhanced the ability to understand interaction between host macrophages and mycobacterium pathogens. Therefore, high-throughput sequencing was used in this study to analyze the transcriptome of MAP-infected bovine mononuclear-macrophages. As mentioned earlier, most of the transcriptional changes in MDM induced by MAP infection occurred within first 6 h . Therefore, in our study, the lncRNA expression profile and lncRNA-mRNA expression network of MDMs 6 h after infection were analyzed. In 2019, RNA-Seq technology was applied to study the changes of mRNA and lncRNA genome expression profile of MAP infected bovine macrophages in vitro by Pooja Gupta et al. . The results showed that differentially expressed lncRNA was highly positively correlated with adjacent differential protein coding genes, and some differentially expressed lncRNA adjacent protein coding genes were related to immune response related pathways. This study is the first time to put forward a new view on the role of lncRNAs and its interaction with mRNAs in the process of cattle MAP infection. In 2021, Andrew Marete et al.  studied the potential role of lncRNA in mononuclear macrophages of dairy cows naturally infected with MAP. The results showed that the lncRNA target genes were significantly enriched in biological process involved in immunity and nucleic acid regulation. The study demonstrates the potential role of lncRNA in host immunity and the potential candidate genes and pathways that may play in response to MAP infection. Coincidentally, our results showed that differentially expressed target genes of lncRNA are also related to inflammation and immune response. Therefore, lncRNA plays an important role in MAP infected macrophages. By analyzing the reasons, the complex numbers of infection of MAP-infected cells in these two papers were 100:1 and 10:1, respectively, while the complex numbers of MAP infection in this paper were 5:1. Therefore, it is speculated that this may be caused by the difference in the complex numbers of infection. Additionally, in this study, the proportion of new lncRNA found on chromosomes 1, 2, 7 and X exceeds 60% of the total lncRNA. X chromosome accounts for the highest proportion, and the newly found lncRNA on chromosome 18 is 54.8%, which are both consistent with the research results of Andrew Marete et al. .
Presently, more and more evidences showed that lncRNA widely attended host’s response to pathogen invasion. Early studies have shown that the presence of several lncRNAs in cattle has an immunomodulatory effect on host’s immune response to Staphylococcus aureus infection. In a cell model of bovine S. aureus mastitis, lncRNA-TUB and H19 mediated Escherichia coli-induced inflammatory factor secretion and S. aureus adhesion to epithelial cells [49, 50]. Gupta et al.  found that in bovine paratuberculosis induced by Mycobacterium paratuberculosis, lncRNA (XLOC-033995) can regulate the expression level of its adjacent inflammatory signal factor TNFAIP3 and participate in immune response related to NF-κB, organelle fission pathways, which affects the inflammatory response process of macrophages to infection, and ultimately regulates the pathogenesis of bovine paratuberculosis. During breeding process, bovine viral diarrhea virus (BVDV) infection can severely affect digestive system and lead to persistent diarrhea and enteritis of a cow. Due to BVDV replicates in the infected Madin-darby bovine kidney cells (MDBK), more and more genes are activated and participate in the immune response, where the number of lncRNAs involved in regulation also increases significantly . Moreover, lncRNA is also related to other diseases, such as liver metastasis, Alzheimer’s disease and colorectal cancer [53,54,55].
In this study, MAP was used to infect bovine monocyte-macrophages for 6 h. The most obvious expression increase in monocyte-macrophages is IL-8, a chemokine that in charge of maintaining balance between physiological reactions and pathological processes in many body systems. Monocytes and macrophages represent principal cellular sources of IL-8 . IL-8 is a strong neutrophil chemotactic agent, which can enhance the killing effect of neutrophils and macrophages on mycobacteria . Simultaneously, these immune cells secrete IL-8 when stimulated by Mycobacterium tuberculosis. In addition, the expression levels of CYP1B1, IMPA1, FXR2, ECE1, LYZ, BHLHE41 and MARCH3 were also significantly changed. CYP1B1 is a member of the CYP superfamily and crucial for endogenous metabolic pathways regulation and is expected to be a therapeutic target for the treatment of metabolic diseases . Owing to involvement in key aspects of cell metabolism, CYP1B1 is used in the diagnosis and predictive diagnosis of various diseases. Based on importance and function of mRNA more than conventional translation, CYP1B1 may also serve as a universal cancer marker. The viral target Endothelin-converting enzyme 1 (ECE1), as an IFNα2a-stimulated protein was exclusively upregulated at the surface of CD4+ T cells, is a new type of cell type specific IFN stimulating factor . BHLHE41 belongs to the basic helix‐loop‐helix family member e41. It is a negative regulatory factor that regulates the inhibition of circadian rhythm, cell differentiation and apoptosis [60, 61]. Membrane-associated RING-CH 3 (MARCH3) is a member of the MARCH ubiquitin ligase family. Several members of this family have been shown to ubiquitinate and down-regulate transmembrane proteins, such as IL-1RAcP . Furthermore, research by Heng Lin et al. showed that MARCH3 is a key negative regulator of IL-1β triggering signal. Overexpressed MARCH3 inhibits IL-1β from activating the expression of NF-κB and inflammatory genes, while inadequate MARCH3 has the opposite effect .
In GO analysis, lncRNAs are mainly enriched in protein homodimerization activity, positive regulation of angiogenesis, interleukin-8 receptor binding, lithium ion binding, trabecular meshwork development, regulation of single stranded viral RNA replication via double stranded DNA intermediate, and endothelial cell–cell adhesion, etc. The results of differential gene enrichment analysis in KEGG pathway showed that most of the genes enriched involved pathways related to the activation of host immune responses to mycobacteria and their regulatory mechanisms. The recognition of Mycobacterium paratuberculosis by macrophages is mediated by host pathogen recognition receptors, including Toll-like receptors (TLRs) and NOD-like receptors (NLRs). TLRs are essential in the recognition of pathogen associated molecular patterns (PAMPs). TLR is a member of mammalian pattern recognition receptor, which is a key component of pathogen recognition mechanism of inflammatory response induced by a microorganism or microbial cell component. TLRs can not only mediate the up-regulation of interleukins, chemokines, costimulatory molecules, adhesion and pro-inflammatory/anti-inflammatory cytokines, but also guide phagocytes to process and present antigens, and induce their own expression . The intracellular NOD-like receptor (NLR) family contains more than 20 members in mammals and is important in the recognition of intracellular ligands. It can be expressed in many cell types, including immune cells and epithelial cells, and thus playing an important role in host innate immune responses and immune homeostasis. Most NLRs are capable of acting as pattern-recognition receptors (PRRs), recognizing various ligands from microbial pathogens, host cells, and environmental sources and activating inflammatory responses. A study showed that NOD2 played an important role in regulating the inflammatory response and intracellular growth of Mycobacterium tuberculosis in human primary macrophages . Additionally, the NF-κB pathway was significantly enriched. The NF-κB pathway helps control cell survival, differentiation and proliferation. Abnormal activation of NF-κB is associated with various diseases such as cancer, autoimmune disease, neurodegenerative disease and cardiovascular disease. Of these diseases, the role of NF-κB signaling in chronic inflammatory and autoimmune diseases, such as Crohn’s disease, ulcerative colitis, and rheumatoid arthritis, is well-defined and known. According to our results, cytokine-cytokine receptor interactions were also significantly enriched. Cytokines, primarily derived from mononuclear phagocytes and other antigen-presenting cells (APCs), are particularly effective in promoting cellular infiltration and resident tissue damage characteristic of inflammation. Cytokines are virtually involved in every aspect of immunity and inflammation, including innate immunity, antigen presentation, myeloid differentiation, cell recruitment and activation, and adhesion molecule expression. Depending on the source of cytokines, their functions are also different that mainly refer to mediating cytotoxic, humoral, cellular or allergic immunity. Moreover, IL-17 signaling pathway was also significantly enriched. The IL-17 cytokine family was named CTLA8 when it was discovered in 1993, and is secreted by T helper cell 17 (Th17). IL-17 plays a key role in protecting the host against extracellular pathogens. The IL-17 family signals via their correspondent receptors and activates downstream pathways that include NF-κB, MAPKs and C/EBPs to induce the expression of antimicrobial peptides, cytokines and chemokines.
In a word, by transcriptome sequencing, a large number of biological information and multiple cell signaling pathways in response to infection was obtained. Based on results of the study, it can provide a basis for comparing the differences of macrophages’ responses to different mycobacteria, and help to clarify the key cellular pathways involved in the early stage of infection. Moreover, screening gene expression information related to infection can also provide a theoretical basis for paratuberculosis diagnosis improvement.
Research on the regulatory relationship network and target gene function in our study is mainly based on lncRNA, analysis and prediction of lncRNA-mRNA transcriptome data. Among the possible modes of action, lncRNA may exert regulatory effect by inhibiting or stimulating gene expression. The potential role of lncRNA as a biomarker has been confirmed in Mycobacterium tuberculosis infection, and this may serve as a promising method in Mycobacterium tuberculosis infection. Furthermore, some lncRNAs are highly regulated by infection, but their functional role still needs to be studied and verified from the aspect of regulation mechanism. Whether the sequenced lncRNAs are specific and sensitive to MAP infection, and whether they can be potential biomarkers for the specific detection of infected animals remain to be verified in the future.
Briefly, in future studies, we will compare and analyze the differences in protein expression between MAP-infected and control groups, and perform association analysis of the three in combination with mRNA and lncRNA to improve the regulatory network of lncRNA. This will provide a new perspective for the pathogenic mechanism study of MAP or novel defense strategies of host cells.
Material and methods
Purification of bovine monocyte derived macrophages
Six healthy Holstein cows (3–6 years old) with no history of Johne’s disease were selected. Collect 200 mL of peripheral blood from the bovine jugular vein and place it in a collection container containing the anticoagulant sodium citrate. The collected blood was then subjected to density gradient centrifugation using lymphocyte separation medium. Peripheral blood mononuclear cells (PBMCs) in the buff layer were collected and resuspended in PBS containing 0.5% BSA and 2 mM EDTA (pH 7.2). Isolation of CD14+ cells was performed using CD14 MACS magnetic beads according to the manufacturers instructions. After magnetic bead sorting, the collected cells were cultured in RPMI 1640 medium containing 10% FBS (Invitrogen, Shanghai, China) and added with 1% non-essential amino acids and 50 ug/mL gentamicin (Sangon Biotech, Shanghai, China). The culture conditions were 37 ℃ and 5% CO2, and the culture density was 1 × 106/well in 24-well cell culture plate. The culture medium was changed on the third day when the monocytes were cultured. At this time, the cell morphology was uniform as shown in Supplementary Fig. 2A. Through flow cytometry, the purity of monocytes after sorting can reach more than 95%, and the results are shown in Supplementary Fig. 2B. The cells were incubated overnight before being used for uptake and replication experiments and allowed to differentiate into the early adherent macrophage phenotype confirmed under the microscope (Supplementary Fig. 2C). When the cells differentiated to the 5th day, the monocyte-derived macrophages were infected by MAP in a 5: 1 infection plural number.
The MAP-10 strain is preserved by the China Institute for Veterinary Drug Control and cultured in Middlebrook 7H9 supplemented with 2 mg/ml Mycobactin J and a 10% Middlebrook OADC enrichment (BD Biosciences, Shanghai, China). The absorbance of the bacterial suspension (OD 600 nm) was measured to evaluate the growth of the bacteria until the collected bacteria could be used for infection. After 8 weeks, bacteria were collected, dissolved in RPMI 1640 cell culture medium, washed and suspended again to prevent bacterial clumps. Finally, colony count was performed  and used to infect cells. For infection, one group was infected with MAP, at a multiplicity of infection of 5:1, the other group was not infected, and each group has three replicates. At 6 h post-infection, the cells were harvested, and the infected and non-infected control MDMs were lysed and stored at -80 ℃ until required for RNA extraction.
RNA extraction, library preparation, sequencing and data analysis
Total RNA was extracted using Trizol reagent (Invitrogen, CA, USA) following the manufacturer’s procedure. The total RNA quantity and purity were analysis of Bioanalyzer 2100 and RNA 6000 Nano LabChip Kit (Agilent, CA, USA) with RIN number > 9.2. Approximately 1 ug of total RNA were used to prepare small RNA library according to protocol of TruSeq Small RNA Sample Prep Kits (Illumina, San Diego, USA).
Ribosomal RNA was removed and approximately 10 ug of total RNA representing the specific adipose type was used to deplete ribosomal RNA. After RNA was purified, the cleaved RNA fragments were reverse transcribed according to the scheme of mRNA-Seq sample preparation kit (Illumina, San Diego, USA) to create the final cDNA library. Then we performed double-terminal sequencing on an Illumina Hiseq 4000 (lc-bio, China) according to the scheme recommended by the supplier.
Transcript assembly and lncRNA screening
After assembling reads using the latest transcript assembly software StringTie, known mRNAs and transcripts smaller than 200 bp were removed, and lncRNA prediction was performed on the remaining transcripts. The prediction software is CPC and CNCI. If these remaining transcripts had the potential to encode proteins, we classified them as novel mRNAs, which were then filtered after defining them as mRNAs. Through some filtering conditions, the lncRNA sequence is finally obtained.
After the final transcriptome was generated, estimate the expression levels of all transcripts, and then lncRNA was screened. Transcripts overlapping with known mRNAs and transcripts shorter than 200 bp were discarded first. Then, CPC and CNCI were used to predict the transcripts with coding potential. All transcripts with CPC score < -1 and CNCI score < 0 were removed. The other types of transcripts with code names I, J, O, U, X are lncRNAs.
Raw RNA-Seq fastq data for the mRNA and lncRNA discussed in this article have been deposited at NCBI and can be found under GEO series accession number GSE193595.
Differential lncRNA expression analysis
For biological replicate results, the Ballgown package in R language was used to perform differential analysis of lncRNA genes assembled and quantified by StringTie. For the screening of differentially expressed genes, the differential lncRNAs were screened with log2 (Foldchange) > 1 (two-fold difference fold) and P value ≤ 0.05 as the threshold.
lncRNA-mRNA network construction, target prediction and functional analysis
According to relevant studies, lncRNA and mRNA are different in many aspects, such as structural characteristics (length distribution, exon number, ORF length) and expression levels. Therefore, this section mainly conducts a systematic comparative analysis on the structural characteristics and expression levels of lncRNAs and mRNAs. To explore the function of lncRNAs, we predicted the cis-target genes of lncRNAs. LncRNAs may play a cis role acting on neighboring target genes. In this study, coding genes in 100,000 upstream and downstream were selected by perl script. Then, we showed functional analysis of the target genes for lncRNAs by using the scripts in house. Significance was expressed as a P value ≤ 0.05.
In this experiment, only cis-regulation predictions were made for the upstream and downstream 100 kbp mRNA and lncRNA. GO and KEGG enrichment analysis was performed on the targeted differential mRNA.
Availability of data and materials
Raw sequencing files and associated metadata have been deposited in NCBI’s Sequence Read Archive (GSE193595).
Groenendaal H, Zagmutt FJ. Scenario analysis of changes in consumption of dairy products caused by a hypothetical causal link between Mycobacterium avium subspecies paratuberculosis and Crohn’s disease. J Dairy Sci. 2008;91(8):3245–58.
Whittington R, Donat K, Weber MF, Kelton D, Nielsen SS, Eisenberg S, Arrigoni N, Juste R, Saez JL, Dhand N, et al. Control of paratuberculosis: who, why and how. A review of 48 countries. BMC Vet Res. 2019;15(1):198.
Kabara E, Coussens PM. Infection of Primary Bovine Macrophages with Mycobacterium avium Subspecies paratuberculosis Suppresses Host Cell Apoptosis. Front Microbiol. 2012;3:215.
Rasmussen P, Barkema HW, Beaulieu E, Mason S, Hall DC. Estimation of the value of Johne’s disease (paratuberculosis) control to Canadian dairy producers. Prev Vet Med. 2021;189:105297.
Caldeira J, Faria A, Diaz-Miranda EA, Zilch TJ, Da CCS, Okano DS, Guimaraes JD, Pena JL, Barbosa WF, Junior AS, et al. Interaction of Mycobacterium avium subsp. paratuberculosis with bovine sperm. Theriogenology. 2021;161:228–36.
Quintas H, Minguez GO, Vila AG, Perez V, Coelho AC. A serosurvey of Mycobacterium avium subsp. paratuberculosis infection of goats in the North of Portugal. Acta Vet Hung. 2021;70:1–8.
Ly A, Kirkeby C, Sergeant E, Plain KM, Smith M, Dhand NK. Comparison of the current abattoir surveillance system for detection of paratuberculosis in Australian sheep with quantitative PCR tissue strategies using simulation modelling. Prev Vet Med. 2021;196:105495.
Wiszniewska-Laszczych A, Liedtke KG, Szteyn JM, Lachowicz T. The Effect of Mycobacterium avium subsp. Paratuberculosis Infection on the Productivity of Cows in Two Dairy Herds with a Low Seroprevalence of Paratuberculosis. Animals (Basel). 2020;10(3):490.
Pourmahdi BM, Haji HM, Ghorbanpoor M, Elhaei SH, Bagheri S, Roveyshedzadeh S. Comparison of Mycobacterium avium subsp. paratuberculosis infection in cattle, sheep and goats in the Khuzestan Province of Iran: Results of a preliminary survey. Vet Med Sci. 2021;7(5):1970–9.
Chiodini RJ, Van Kruiningen HJ. Eastern white-tailed deer as a reservoir of ruminant paratuberculosis. J Am Vet Med Assoc. 1983;182(2):168–9.
Reyes-Garcia R, Perez-de-la-Lastra JM, Vicente J, Ruiz-Fons F, Garrido JM, Gortazar C. Large-scale ELISA testing of Spanish red deer for paratuberculosis. Vet Immunol Immunopathol. 2008;124(1–2):75–81.
Jessup DA, Abbas B, Behymer D. Paratuberculosis in tule elk in California. J Am Vet Med Assoc. 1981;179(11):1252–4.
Salgado M, Herthnek D, Bolske G, Leiva S, Kruze J. First isolation of mycobacterium avium subsp. Paratuberculosis from wild guanacos (Lama guanicoe) on Tierra del Fuego Island. J Wildl Dis. 2009;45(2):295–301.
Mcclure HM, Chiodini RJ, Anderson DC, Swenson RB, Coutu T. Mycobacterium paratuberculosis Infection in a Colony of Stumptail Macaques (Macaca arctoides). J Infect Dis. 1987;155(5):1011–9.
Sechi LA, Dow CT. Mycobacterium avium ss. paratuberculosis Zoonosis - The Hundred Year War - Beyond Crohn’s Disease. Front Immunol. 2015;6:96.
Niegowska M, Rapini N, Piccinini S, Mameli G, Caggiu E, Manca BM, Sechi LA. Type 1 Diabetes at-risk children highly recognize Mycobacterium avium subspecies paratuberculosis epitopes homologous to human Znt8 and Proinsulin. Sci Rep. 2016;6:22266.
Bo M, Erre GL, Niegowska M, Piras M, Taras L, Longu MG, Passiu G, Sechi LA. Interferon regulatory factor 5 is a potential target of autoimmune response triggered by Epstein-barr virus and Mycobacterium avium subsp. paratuberculosis in rheumatoid arthritis: investigating a mechanism of molecular mimicry. Clin Exp Rheumatol. 2018;36(3):376–81.
Hayashi F, Isobe N, Cossu D, Yokoyama K, Sakoda A, Matsushita T, Hattori N, Kira JI. Elevated mycobacterium avium subsp. paratuberculosis (MAP) antibody titer in Japanese multiple sclerosis. J Neuroimmunol. 2021;360:577701.
Steuer P, Tejeda C, Moroni M, Verdugo C, Collins MT, Salgado M. Attempted Control of Paratuberculosis in Dairy Calves by Only Changing the Quality of Milk Fed to Calves. Animals (Basel). 2021;11(9):2569.
Hassan AA, Khan I, Ganz S, Wehrend A, Failing K, Eisenberg T, Abdulmawjood A, Bulte M. Assessing efficacy of N-Acetyl-l-Cysteine-Sodium Hydroxide on bacterial viability and enhanced recovery of Mycobacterium avium subsp. paratuberculosis from bovine colostrum. J Microbiol Methods. 2020;175:105968.
Kuenstner L, Kuenstner JT. Mycobacterium avium ssp. paratuberculosis in the Food Supply: A Public Health Issue. Front Public Health. 2021;9:647448.
Bermudez LE, Petrofsky M, Sommer S, Barletta RG. Peyer’s Patch-Deficient Mice Demonstrate That Mycobacterium avium subsp. paratuberculosis Translocates across the Mucosal Barrier via both M Cells and Enterocytes but Has Inefficient Dissemination. Infect Immun. 2010;78(8):3570–7.
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.
Behar SM, Divangahi M, Remold HG. Evasion of innate immunity by Mycobacterium tuberculosis: is death an exit strategy? Nat Rev Microbiol. 2010;8(9):668–74.
Lin Y, Shen Y, Chen J, Hu C, Zhou Z, Yuan C. The Function of LncRNA FTX in Several Common Cancers. Curr Pharm Des. 2021;27(20):2381–6.
Yoon JH, Abdelmohsen K, Gorospe M. Functional interactions among microRNAs and long noncoding RNAs. Semin Cell Dev Biol. 2014;34:9–14.
Li S, Cao Y, Zhang H, Lu X, Wang T, Xu S, Kong T, Bo C, Li L, Ning S, et al. Construction of lncRNA-Mediated ceRNA Network for Investigating Immune Pathogenesis of Ischemic Stroke. Mol Neurobiol. 2021;58(9):4758–69.
Cao M, Li H, Zhao J, Cui J, Hu G. Identification of age- and gender-associated long noncoding RNAs in the human brain with Alzheimer’s disease. Neurobiol Aging. 2019;81:116–26.
Ghafouri-Fard S, Shirvani-Farsani Z, Hussen BM, Taheri M. The critical roles of lncRNAs in the development of osteosarcoma. Biomed Pharmacother. 2021;135:111217.
Radhakrishnan R, Kowluru RA. Long Noncoding RNA MALAT1 and Regulation of the Antioxidant Defense System in Diabetic Retinopathy. Diabetes. 2021;70(1):227–39.
Taniue K, Akimitsu N. The Functions and Unique Features of LncRNAs in Cancer Development and Tumorigenesis. Int J Mol Sci. 2021;22(2):632.
He Y, Dan Y, Gao X, Huang L, Lv H, Chen J. DNMT1-mediated lncRNA MEG3 methylation accelerates endothelial-mesenchymal transition in diabetic retinopathy through the PI3K/Akt/mTOR signaling pathway. Am J Physiol Endocrinol Metab. 2021;320(3):E598–608.
Cai R, Sun Y, Qimuge N, Wang G, Wang Y, Chu G, Yu T, Yang G, Pang W. Adiponectin AS lncRNA inhibits adipogenesis by transferring from nucleus to cytoplasm and attenuating Adiponectin mRNA translation. Biochim Biophys Acta Mol Cell Biol Lipids. 2018;1863(4):420–32.
Mukherjee N, Calviello L, Hirsekorn A, de Pretis S, Pelizzola M, Ohler U. Integrative classification of human coding and noncoding genes through RNA metabolism profiles. Nat Struct Mol Biol. 2017;24(1):86–96.
Cabili MN, Dunagin MC, McClanahan PD, Biaesch A, Padovan-Merhar O, Regev A, Rinn JL, Raj A. Localization and abundance analysis of human lncRNAs at single-cell and single-molecule resolution. Genome Biol. 2015;16:20.
Sikora M, Marycz K, Smieszek A. Small and Long Non-coding RNAs as Functional Regulators of Bone Homeostasis, Acting Alone or Cooperatively. Mol Ther Nucleic Acids. 2020;21:792–803.
Gaiti F, Degnan BM, Tanurdzic M. Long non-coding regulatory RNAs in sponges and insights into the origin of animal multicellularity. RNA Biol. 2018;15(6):696–702.
Chao Y, Jin J, Wang L, Jin X, Yang L, Zhang B. Transcriptome Analysis of lncRNA-mRNA Interactions in Chronic Atrophic Gastritis. Front Genet. 2020;11:612951.
Font-Cunill B, Arnes L, Ferrer J, Sussel L, Beucher A. Long Non-coding RNAs as Local Regulators of Pancreatic Islet Transcription Factor Genes. Front Genet. 2018;9:524.
Kanehisa MAGS. KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Res. 2000;28:27–30.
Kanehisa M. Toward understanding the origin and evolution of cellular organisms. Protein Sci. 2019;28(11):1947–51.
Kanehisa M, Furumichi M, Sato Y, Kawashima M, Ishiguro-Watanabe M. KEGG for taxonomy-based analysis of pathways and genomes. Nucleic Acids Res. 2022;gkac963.
Borujeni MP, Hajikolaei MRH, Ghorbanpoor M, Sahar HE, Roveyshedzadeh S. Comparison of Mycobacterium avium subsp. paratuberculosis infection in cattle, sheep and goats in the Khuzestan Province of Iran: Results of a preliminary survey. Vet Med Sci. 2021;7(5):1970–9.
Abbas M, Munir M, Khaliq SA, Haq M, Khan MT, Qureshi Z. Detection of Paratuberculosis in Breeding Bulls at Pakistani Semen Production Units: A Continuous Source of Threat. Isrn Veterinary Science. 2011;2011(1):501235.
Eisenberg S, Nielen M, Santema W, Houwers DJ, Heederik D, Koets AP. Detection of spatial and temporal spread of subsp. in the environment of a cattle farm through bio-aerosols. Vet Microbiol. 2010;143(2–4):284–92.
Casey ME, Meade KG, Nalpas NC, Taraktsoglou M, Browne JA, Killick KE, Park SD, Gormley E, Hokamp K, Magee DA, et al. Analysis of the Bovine Monocyte-Derived Macrophage Response to Mycobacterium avium Subspecies Paratuberculosis Infection Using RNA-seq. Front Immunol. 2015;6:23.
Gupta P, Peter S, Jung M, Lewin A, Hemmrich-Stanisak G, Franke A, von Kleist M, Schutte C, Einspanier R, Sharbati S, 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(1):1571.
Marete A, Ariel O, Ibeagha-Awemu E, Bissonnette N. 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. 2021;8:209.
Wang H, Wang X, Li X, Wang Q, Qing S, Zhang Y, Gao MQ. A novel long non-coding RNA regulates the immune response in MAC-T cells and contributes to bovine mastitis. FEBS J. 2019;286(9):1780–95.
Li X, Wang H, Zhang Y, Zhang J, Qi S, Zhang Y, Gao MQ. Overexpression of lncRNA H19 changes basic characteristics and affects immune response of bovine mammary epithelial cells. PeerJ. 2019;7:e6715.
Gupta P, Peter S, Jung M, Lewin A, Bruegge JZ. Analysis of long non-coding RNA and mRNA expression in bovine macrophages brings up novel aspects of Mycobacterium avium subspecies paratuberculosis infections. Sci Rep UK. 2019;9(1):1571.
Qm A, Llb F, Yan TE, Qiang FD, Sheng LB, Sh C, Jq B, Cc B, Wei NC. Analyses of long non-coding RNAs and mRNA profiling through RNA sequencing of MDBK cells at different stages of bovine viral diarrhea virus infection. 2017.
Huang L, Lin H, Kang L, Huang P, Huang J, Cai J, Xian Z, Zhu P, Huang M, Wang L, et al. Aberrant expression of long noncoding RNA SNHG15 correlates with liver metastasis and poor survival in colorectal cancer. J Cell Physiol. 2019;234(5):7032–9.
Zhang Z. Long non-coding RNAs in Alzheimer’s disease. Curr Top Med Chem. 2016;16(5):511–9.
Gupta RA, Shah N, Wang KC, Kim J, Horlings HM, Wong DJ, Tsai MC, Hung T, Argani P, Rinn JL, et al. Long non-coding RNA HOTAIR reprograms chromatin state to promote cancer metastasis. Nature. 2010;464(7291):1071–6.
Sikora J, Smycz-Kubanska M, Mielczarek-Palacz A, Kondera-Anasz Z. Abnormal peritoneal regulation of chemokine activation-The role of IL-8 in pathogenesis of endometriosis. Am J Reprod Immunol. 2017;77(4):e12622.
Krupa A, Fol M, Dziadek BR, Kepka E, Wojciechowska D, Brzostek A, Torzewska A, Dziadek J, Baughman RP, Griffith D, et al. Binding of CXCL8/IL-8 to Mycobacterium tuberculosis Modulates the Innate Immune Response. Mediators Inflamm. 2015;2015:124762.
Carrera AN, Grant M, Zordoky BN. CYP1B1 as a therapeutic target in cardio-oncology. Clin Sci (Lond). 2020;134(21):2897–927.
Soday L, Potts M, Hunter LM, Ravenhill BJ, Weekes MP. Comparative Cell Surface Proteomic Analysis of the Primary Human T Cell and Monocyte Responses to Type I Interferon. Front Immunol. 2021;12:600056.
Nagata T, Minami K, Yamamoto M, Hiraki T, Idogawa M, Fujimoto K, Kageyama S, Tabata K, Kawahara K, Ueda K, et al. BHLHE41/DEC2 Expression Induces Autophagic Cell Death in Lung Cancer Cells and Is Associated with Favorable Prognosis for Patients with Lung Adenocarcinoma. Int J Mol Sci. 2021;22(21):11509.
Sato F, Bhawal UK, Yoshimura T, Muragaki Y. DEC1 and DEC2 Crosstalk between Circadian Rhythm and Tumor Progression. J Cancer. 2016;7(2):153–9.
Chen R, Li M, Zhang Y, Zhou Q, Shu HB. The E3 ubiquitin ligase MARCH8 negatively regulates IL-1beta-induced NF-kappaB activation by targeting the IL1RAP coreceptor for ubiquitination and degradation. Proc Natl Acad Sci U S A. 2012;109(35):14128–33.
Lin H, Gao D, Hu MM, Zhang M, Wu XX, Feng L, Xu WH, Yang Q, Zhong X, Wei J, et al. MARCH3 attenuates IL-1beta-triggered inflammation by mediating K48-linked polyubiquitination and degradation of IL-1RI. Proc Natl Acad Sci U S A. 2018;115(49):12483–8.
Bhide MR, Mucha R, Mikula IJ, Kisova L, Skrabana R, Novak M, Mikula IS. Novel mutations in TLR genes cause hyporesponsiveness to Mycobacterium avium subsp. paratuberculosis infection. BMC Genet. 2009;10:21.
Brooks MN, Rajaram MV, Azad AK, Amer AO, Valdivia-Arenas MA, Park JH, Núñez G, Núñez G. NOD2 controls the nature of the inflammatory response and subsequent fate of Mycobacterium tuberculosis and M. bovis BCG in human macrophages. Cell Microbiol. 2011;13:402–18.
Lamont EA, Talaat AM, Coussens PM, Bannantine JP, Grohn YT, Katani R, Li LL, Kapur V, Sreevatsan S. Screening of Mycobacterium avium subsp. paratuberculosis mutants for attenuation in a bovine monocyte-derived macrophage model. Front Cell Infect Microbiol. 2014;4:87.
We are grateful for the owner of the farm who provided cattle blood. We also thank the editor and reviewers for their helpful comments and insights.
This work was supported by funding from the Natural Science Foundation of Jilin Province (No.20220101349JC) and the National Natural Science Foundation of China (Grants: 30471285, 30972190 and 31572555).
Ethics approval and consent to participate
All experimental and sample collection procedures were approved by Animal Care and Use Committee of Jilin Agricultural University, with approval number GB/T 35892–2018. All methods were carried out under relevant guidelines and regulations (declaration of helsinki). All methods are reported in accordance with ARRIVE guidelines for the reporting of animal experiments.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
About this article
Cite this article
Bao, Y., Wu, S., Yang, T. et al. Analysis of long non-coding RNA expression profile of bovine monocyte-macrophage infected by Mycobacterium avium subsp. paratuberculosis. BMC Genomics 23, 768 (2022). https://doi.org/10.1186/s12864-022-08997-5
- M. avium subsp. paratuberculosis
- High-throughput sequencing