Mycoplasma hyopneumoniae strain 98, a Danish field isolate strain (provided by Dr. N. Friis, National Veterinary Laboratory, Copenhagen) was used. The strain was grown in 100 ml FRIIS medium (Friis, 1975) in 250 ml closed glass bottles at 37°C with agitation (100 RPM). Cultures were sampled during exponential growth indicated by an increase in titer measured using flow cytometry (FACSMicroCount, BD). Mycoplasma cells were pelleted using centrifugation (3 min at 9000 g’s) and the cell pellet was directly frozen (<-15°C) for DNA extraction or submerged in RNA later (Ambion) and stored at 2–8°C for RNA extraction. DNA extraction for genome sequencing from the pelleted frozen cells was done using the Gentra Puregene bacterial kit (Qiagen). Genome sequencing was done using Illumina HiSeq2500 (paired-end, 500 MB, 250 bp read length). Genome assembly was done using the Ray algorithm (Boisvert et al., 2010). Annotation of the genome was done using SAPP (Koehorst et al., 2016) in which Prodigal v. 2.6.2 (Hyatt et al., 2010) was used for gene calling and InterProScan v.5.17-56.0 (Jones et al., 2014) for protein domain annotation. A reciprocal best-blast hit analysis (Wolf and Koonin, 2012) was done using BLAST + (Camacho et al., 2009) to find orthologous proteins in M. hyopneumoniae strain 232 (Minion et al., 2004), filtering for matches with a minimum of 70% identity and 50% query coverage. For interpretation of the differentially expressed bacterial genes we mention the locus tags of genes present in strain 232 found with the BLAST + analysis. NCBI reference genomes used to compare the protein domain repertoire were: NC_006360 (strain 232), NC_007295 (strain J), NC_007332 (strain 7448), NC_017509.1 (strain 168), NC_021283 (strain 168L), and NC_021831 (strain 7422). For strain 11 accession MWWN00000000 was used. The genome sequence of strain 98 is available from the NCBI under accession number WTQC00000000.

Animal studies were performed after approval by an ethical commission and according to national regulations in The Netherlands. In total thirteen pigs were used for this study from three animal trials performed at MSD Animal Health. Nine healthy pigs from a M. hyopneumoniae free herd were challenged intra-tracheally on two consecutive days with 10 ml of M. hyopneumoniae strain 98 culture containing ± 10⁷ CCU/ml. Pigs were anesthetized, euthanized and exsanguinated 3 weeks after challenge. Within a short time-period (<5 min) after the necropsy the lungs were removed and tissue samples were obtained and directly submerged in RNA later (Ambion) or infected lobes were flushed with RNA later using a method previously described (Mengeling et al., 1995) but adapted for this specific purpose. Briefly, an infected lobe was selected for sampling based on the presence of lesions in the lobe. A plastic pipette (7 ml total volume) was used to dispense 5 ml of RNA later (Ambion) into the bronchus toward the selected lobe. Directly after dispensing, the fluid was retrieved from the lobe and stored at 2–8°C for further processing. Four healthy pigs were necropsied unchallenged at 7 weeks ± 3 days of age.

RNA extraction from tissue samples was done using the RNeasy mini kit (Qiagen) following the manufacturer’s protocol. Lesions in infected lungs were identified by discoloration of tissue and from the affected tissue, total RNA was extracted for sequencing. RNA extraction from flush samples (0.75 ml) was done with 7.5 ml of Trizol LS reagent (Ambion) following the manufacturer’s protocol. For bacterial pellets, RNA later was removed and 3.75 ml Trizol LS reagent was added. RNA was precipitated using glycogen as co-precipitant as described in the Trizol extraction protocol. After RNA re-suspension, the quantity was determined using NanoDrop. RNA samples were treated with Turbo DNase (Ambion) to remove potential DNA contamination following the manufacturer’s protocol. Enrichment of bacterial RNA in the flush samples and culture samples was done using the MICROBEnrich kit (Ambion) and the quality of all samples was determined using Experion RNA StdSense (Bio-Rad).

Kallisto (Bray et al., 2016) was used to determine estimated counts on gene transcripts using pseudoalignment to a pig reference transcript file from Ensembl (all cDNA version 11.1, FTP location)¹. To analyze the read distribution of prokaryotic reads on the M. hyopneumoniae strain 98 genome, rRNA was removed using the Ribo-Zero Gold rRNA Removal Kit (Epidemiology). rRNA-depleted RNA was fragmented to an average length of 100–200 base pairs and converted to double-stranded complementary DNA (cDNA). Strand specific library preparation was done using a protocol based on the “dUTP (deoxyuridine triphosphate) method.” The Illumina stranded TruSeq RNA-seq library preparation kit was used. Sequencing of the library was done using the Illumina HiSeq 2500: single-end reads, 50 cycles. Quality assessment of reads was done as previously described (Kamminga et al., 2017b). Reads were aligned to the strain 98 genome using TopHat v2.1.0 (Kim et al., 2013) with the standard settings for strand-specific alignment enabled (fr-firststrand). Accepted hits found by TopHat were sorted using samtools v1.3.1 (Li et al., 2009), from the sorted hits a BED-file was created using bamToBed and coverage per gene was calculated using coverageBed, both programs from the bedtools suite v2.17.0 (Quinlan and Hall, 2010). Reads mapping to non-coding RNAs, were determined using the same methods after creation of a ncRNA annotation file for the strain 98 genome based on the best blast hit (Camacho et al., 2009) with the annotated ncRNAs in strain 232 (Lloréns-Rico et al., 2016). Raw data files used for this study are available from the NCBI bioproject PRJNA593525.

The R package edgeR (Robinson et al., 2009; McCarthy et al., 2012) was used for analysis of differentially expressed genes and differentially expressed ncRNAs (for simplicity we refer to both as genes) between the in vitro and in vivo condition. Only genes with a read count of >100 counts per million (CPM) in two or more datasets were kept as described by Rienksma et al. (2015). Data was normalized for RNA composition by calculating scale factors (based on the trimmed mean of M-values; Robinson and Oshlack, 2010) which were used to correct the library size. Common dispersion was calculated based on conditional maximum likelihood, correcting for library size differences based on pseudocounts estimated per quantile, first by assuming a poisson distribution to calculate the common estimated dispersion. This initially calculated common dispersion was used to get a final estimate for the pseudocounts and a final estimate of the common dispersion (Robinson and Smyth, 2008). Two types of variation contribute to the total variance in the probability distribution: technical and biological variation. Biological variation (BCV) represents the variation between biological replicates that would remain if sequencing depth is infinite. Dispersion per gene (biological variation) was calculated using an empirical Bayes estimation with weighted likelihood (Robinson and Smyth, 2007). Differential expression was calculated for each gene using a pairwise exact testing method taking into account the tag-wise dispersion estimates. False discovery rates were controlled using the algorithm by Benjamini and Hochberg (1995). Genes with an FDR < 0.01 and a LOG2 fold change larger than 2 were considered to be significantly different. For eukaryotic genes an FDR threshold of 0.01 and fold change threshold of 1.5 was used to check for significance. A more conservative cut-off value for the fold change of differentially expressed bacterial genes was used because low RNA quantities were obtained from infected lung flushes.

Genome assembly using the filtered Illumina FastQ sequence reads (average Phred quality score: 36.17) resulted in a draft genome sequence of 19 scaffolds. The total genome size was 880.620 bp and the GC% was 28.5%. We compared the protein domain content of strain 98 to the domain content of other Mycoplasma hyopneumoniae strains (232, J, 11, 7448, 7422, 168, and 168L, Supplementary Table S1) obtained after genome annotation using SAPP (Koehorst et al., 2016; Kamminga et al., 2017a). We found a total core domainome size of 843 unique protein domains and an accessory domainome of 23 protein domains (Supplementary Table S2). With a reciprocal-best blast hit analysis 601 orthologous proteins were identified in strain 98 when compared to strain 232. Based on the composition of the accessory domainome, strain 98 is mostly related to M. hyopneumoniae strain 7448 (Supplementary Figure S1). Overall, the differences were small and no novel protein domains were identified in the strain 98 genome.

Comparison of RNA sequencing datasets of three infected tissue samples with four non-infected tissue samples showed 424 genes that were significantly differentially expressed in the infected animals versus the controls (Supplementary Table S3). Reactome pathways containing these 424 genes (filtered for significance based on p < 0.01) were analyzed for overrepresentation of differentially expressed genes. In total 17 pathways were enriched (Table 1), related to interleukin signaling, detection of chemokines, tissue damage (regulation TLR, O-fucosylation = Peters syndrome, TSR), STAT3 activation, B cell receptor (BCR) signaling and complement activation (Figure 1). Specifically, interleukin-10, interleukin-4, interleukin-13, and interleukin-18 related pathways were found up-regulated. Increased expression of pro-inflammatory cytokines was expected in infected lung tissue and increased levels of IL-4, IL-10, and IL-18 were previously determined using immunohistochemical and ELISA methods (Rodríguez et al., 2004; Thanawongnuwech et al., 2004; Choi et al., 2006; Muneta et al., 2006). Increased expression of IL-13 was not previously described during M. hyopneumoniae infection and should be verified by measuring protein levels. We did not identify upregulation of all interleukins previously found upregulated during infection (IL-1, IL-2, IL-6, IL-8, IL-12) or of TNF-alpha production. This could be because we used a different analysis method, had animals with a different genetic background or because we specifically sequenced a lesion where infection is already causing tissue damage.

The pathway for regulation of Toll-like receptors by endogenous ligands was upregulated in infected lung tissue. Endogenous ligands are formed when tissue damage occurs. Specific genes identified were TLR1, TLR2, TLR6, LBP, and SFTPA1. A role for TLR2 and TLR6 was described previously (Muneta et al., 2003; Sun et al., 2013) and upregulation of these genes is caused by the presence of lipoproteins on the cell membrane of M. hyopneumoniae. A possible upregulation of lipopolysaccharide binding protein (LBP) was not previously described and needs to be further investigated. LBP might be upregulated because glycolipids are present on the M. hyopneumoniae membrane and/or because oxidized phospholipids are present that interact with LBP (Erridge and Spickett, 2013). SFTPA1 produces pulmonary surfactant-associated protein A1 which could play a role in killing of Mycoplasma by alveolar macrophages (Hickman-Davis et al., 1999). Note that because only tissue from lesions was sequenced by RNAseq, it was not possible to distinguish the immune response to exogenous compounds from the response to endogenous compounds.

Whether the complement system plays a role in protection against M. hyopneumoniae is not known. However, pulmonary alveolar type II epithelial cells have been shown to produce complement proteins in the lung (Pandya and Wilkes, 2014). We identified an increased expression of the C3 and C5 components of the complement system which could be caused by an increased concentration of pro-inflammatory cytokines in the lung. The role of the complement system during M. hyopneumoniae infection should be further investigated.

We found upregulation of a pathway which plays a role in B cell receptor signaling and B cell development (RUNX1 regulates transcription of genes involved in BCR signaling). Increased expression in pathways related to B cell receptor signaling is expected when there is a serological response to infection. In previous studies, elevated TNF levels have been measured in infected lungs which could stimulate activation of B cells (Thacker et al., 2000; Thanawongnuwech and Young, 2001). At the moment of necropsy, seroconversion was found for all M. hyopneumoniae challenged pigs using the IDEXX M. hyo Ab test.

Lung tissue samples were obtained from infected animals and sequenced with high capacity. Here, 0.1% of the reads were found to match with the M. hyopneumoniae genome (Supplementary Table S4). As a rule of thumb, we assumed that sufficient coverage of a bacterial genome for analysis of gene expression levels is obtained with ≥ 1 million non-rRNA bacterial reads (Westermann et al., 2012). Since insufficient numbers of bacterial reads were obtained, dual RNASeq was considered not economically feasible for this sample type.

To overcome this hurdle, we developed alternative methods to isolate bacterial RNA from lungs infected with M. hyopneumoniae (Figure 2A). Application of differential lysis techniques to isolate bacterial RNA (Rienksma et al., 2015; Robbe-Saule et al., 2017), was not an option for Mycoplasma hyopneumoniae since this bacterium does not contain a cell wall. As previously reported, we identified with immunohistochemistry that M. hyopneumoniae is found specifically on the surface of the respiratory epithelium (Figure 2A). To isolate sufficient amounts of bacteria, we decided to flush an infected lung lobe with RNA later and isolate total RNA from the lung flush. Again, we tested if this would be a feasible method to perform dual RNASeq but from non-infected lung flush samples insufficient host RNA concentrations were isolated. We hypothesized that low cell numbers were present at the epithelial surface of non-infected lungs where infected lungs contained a higher cell number as a result of infection. Since we no longer pursued dual RNAseq, an enrichment step for bacterial RNA was implemented in which pig mRNA was removed based on the presence of a poly-A tail (Figure 2A).

The average RNA quantity isolated from infected lung flush samples was low: 4.2 μg (range: 0.17–18.51 μg). Due to the low quantity of RNA it was not possible to assess the RNA quality based on the ratio between 23S and 16S rRNA (Supplementary Figures S2–S6). Therefore, we analyzed the quality of the in vivo data with multivariate data analysis. A ten times higher read mapping percentage was achieved by sequencing RNA from a flush sample and even higher read mapping percentages were reached by applying the enrichment step (Figure 2B and Supplementary Table S4). To assess if the enrichment step influenced the expression levels for bacterial genes, we sequenced a flush sample from the same pig with and without applying the enrichment step, all other treatment steps were done in parallel. There was a strong positive correlation between the expression levels in the non-enriched sample compared to the enriched sample (Pearson coefficient = 0.97, Supplementary Figure S7) which showed that the enrichment step and the repeat of the other purification steps did not influence the bacterial read distribution.

We further analyzed the variation in the read distribution per bacterial gene in the various sample types with principal component analysis (PCA) based on the read counts per gene. This analysis showed a separation of in vitro samples from in vivo samples and again showed that enrichment for bacterial RNA did not influence the general expression landscape in biological replicates (Supplementary Figure S8). This analysis also showed that 69% of the variance in the dataset followed from the sample origin and only 19% variance was found between in vivo samples, showing high consistency for both sample types. The method we developed enabled sequencing of bacterial mRNA from infected lung tissue with sufficient read numbers.

A disadvantage of the lung flush method is that the exact location where the bacterial cells resided cannot be determined. We consider it likely that by applying shear force when flushing the lung, cells were removed that were both unbound and bound to the cilia, however an exact ratio could not be determined. It is also not clear if multilayered cell structures (biofilms) were removed by flushing or mainly single cells not bound to the cilia. However, since biofilms were present on the lung surface where cilia were degraded (Raymond et al., 2018), we consider it likely that by flushing both loose cells as well as biofilm structures can be removed for sequencing. The PCA method did not show a major variation between flush samples indicating that, although we do not know exactly the location of the bacterial cells, the lung flush method can be used to consistently obtain bacterial culture from multiple lobes in the pig lung.

We analyzed the distribution of sense and antisense mapping reads per gene and found on average 97.5 ± 1.2% sense mapping reads (Supplementary Table S5). Most genes needed for cilium adhesion were highly expressed during infection (Table 2). This was expected because the bacterium specifically binds the ciliated epithelium (Figure 2A; Zielinski et al., 1990; Blanchard et al., 1992; Young et al., 2000). The same genes were also highly expressed in vitro which may suggest constitutive expression. Surprisingly, the gene encoding the cilium adhesin P102 (17_1, MHP_RS00920) was repressed during infection (see below). Genes related to cell division and transport, specifically a subunit of the ascorbate transporter and MFS transporter, were also found highly expressed in both conditions. Key metabolic genes found highly expressed under both conditions were related to the final stages of glycolysis (pyruvate dehydrogenase and lactate dehydrogenase) and myo-inositol metabolism. This indicates an important role for these metabolic pathways under both growth conditions. Our results agreed well with data reported by Siqueira et al. (2014) of highly expressed genes in M. hyopneumoniae strain 7448 determined with RNA sequencing.

The overlap observed in highly expressed genes could be partly explained by the choice of the cultivation medium. In our study Friis medium was used for in vitro cultivation of M. hyopneumoniae. This medium contains animal serum, animal derived peptones and yeast extract. The composition of the growth medium will influence the in vitro transcriptome of M. hyopneumoniae. Since there are components in serum that will also be present in the lung when tissue damage occurs, bacteria could encounter similar components in the lung as in the medium which could result in similar expression levels for certain genes. If a medium was chosen with less animal-derived components, more differences could be observed between in vitro and in vivo growth, but such a medium is not readily available.

Sense mapping reads were used to analyze differentially expressed genes after filtering for genes with a read count of >100 counts per million (CPM) in two or more datasets (Rienksma et al., 2015). By filtering we removed 203 genes (29% of genes in the genome). Based on FDR < 0.01 and LOG2 fold change larger than two, we found 23 genes up-regulated in vivo, 30 genes down-regulated in vivo and 445 genes not differentially expressed (Figure 3 and Supplementary Table S6). The average biological coefficient of variance (BCV) calculated for all genes in the dataset using EdgeR was 0.30, meaning that on average 30% of the variation in gene abundance is the result of biological variation between biological replicates. This value for BCV is higher than expected for experiments with genetically identical organisms (McCarthy et al., 2012). The higher variation is possibly the result of varying conditions in the lung, the PCA analysis showed that the majority of the variation within samples types is present in the in vivo datasets (Supplementary Figure S8).

Previously microarrays have been used to study the transcriptome of M. hyopneumoniae during lung infection (Madsen et al., 2008). We found no correlation (Spearman’s coefficient = -0.14, Supplementary Figure S10) between fold changes for 162 differentially expressed genes reported in this microarray study and fold changes observed in our study. This could be due to: (1) different M. hyopneumoniae strains were used, (2) pigs with a different genetic background and age were used, (3) different methods were used to isolate bacteria and bacterial RNA and (4) presence of host RNA could have caused non-specific binding in the microarray study. It was recently shown that differentially expressed genes at low transcript levels could be identified with RNAseq and validated with qPCR while differentially expressed genes identified with microarrays could not be validated with qPCR when transcript levels were low (Wang et al., 2014). Since the lung flush samples contained a high quantity of background RNA and low bacterial transcript levels, RNAseq seems a more suitable method to determine bacterial transcript levels in this sample type.