New insights into host adaptation to swine respiratory disease revealed by genetic differentiation and RNA sequencing analyses

Abstract Swine respiratory disease (SRD) causes massive economic losses in the swine industry and is difficult to control and eradicate on pig farms. Here, we employed population genetics and transcriptomics approaches to decipher the molecular mechanism of host adaptation to swine respiratory disease. We recorded two SRD‐related traits, the enzootic pneumonia‐like (EPL) score and lung lesion (LL) levels, and performed four body weight measurements, at ages of 150, 180, 240, and 300 days, in a Chinese Bamaxiang pig herd (n = 314) raised under consistent indoor rearing conditions. We divided these animals into disease‐resistant and disease‐susceptible groups based on the most likely effects of both SRD‐related traits on their weight gain, and performed genetic differentiation analyses in these two groups. Significant loci showing the top 1% of genetic differentiation values, exceeding the threshold of p = 0.005 set based on 1,000‐times permutation tests, were defined as candidate regions related to host resistance or susceptibility to SRD. We identified 107 candidate genes within these regions, which are mainly involved in the biological processes of immune response, fatty acid metabolism, lipid metabolism, and growth factor signaling pathways. Among these candidate genes, TRAF6, CD44, CD22, TGFB1, CYP2B6, and SNRPA were highlighted due to their central regulatory roles in host immune response or fat metabolism and their differential expression between healthy lung tissues and lung lesions. These findings advance our understanding of the molecular mechanisms of host resistance or susceptibility to respiratory disease in pigs and are of significance for the breeding pigs resistant to respiratory disease in the swine industry.


| INTRODUC TI ON
Swine respiratory disease (SRD), also referred as the swine/porcine respiratory disease complex (SRDC/PRDC), is prevalent in modern intensive pig farms worldwide and can cause poor porcine growth performance, in turn leading to serious economic losses to the swine industry (Opriessnig, Gimenez-Lirola, & Halbur, 2011). This disease is arguably the most important health concern for swine producers at present (Brockmeier, Halbur, & Thacker, 2002). Multiple infectious agents are usually involved in swine respiratory disease (Brockmeier et al., 2002;Opriessnig et al., 2011). These agents can be divided into primary bacterial pathogens, such as Mycoplasma hyopneumoniae (Thacker, 2004) and Bordetella bronchiseptica (Brockmeier, 2004); primary viral pathogens, such as porcine circovirus type 2 (Kim, Chung, & Chae, 2003) and porcine reproductive and respiratory syndrome virus (Thacker, Halbur, Ross, Thanawongnuwech, & Thacker, 1999); and secondary or opportunistic pathogens, such as Pasteurella multocida (Ciprian et al., 1988). In addition, environmental factors and management conditions, such as a high animal stocking density and negligent housing maintenance, could increase risk to affect swine respiratory disease (Maes et al., 2008;Stark, 2000).
On pig farms, swine respiratory disease caused by viral pathogens can often be effectively controlled by timely vaccination. However, swine respiratory disease caused by microbial pathogens, especially M. hyopneumoniae, is difficult to prevent and is prone to developing into a chronic, endemic disease that recessively exists in pig populations over the long term, resulting in reduction of the porcine growth rate and decreasing feed efficiency (Brockmeier et al., 2002). Under natural farming conditions, pigs are always subject to an adverse disease environment for long-term recessive swine respiratory disease caused by bacterial agents, most likely by M. hyopneumoniae. It was shown that 93.4% of pig lungs exhibited pneumonic lesions among 4,508 pigs at slaughter, and 98.1% of enzootic pneumonia (EP)-affected lungs showed lesions similar to pneumonic lesions caused by M. hyopneumoniae (Wallgren, Beskow, Fellström, & Renström, 1994), which demonstrated that the pigs had been subjected to an adverse SRD environment.
Great efforts have been devoted to controlling the prevalence of swine respiratory disease, including improving the feeding environment and management, implementing vaccination schedules, and creating pathogen eradication schemes (Bargen, 2004;Maes, Verdonck, Deluyker, & Kruif, 1996;Tzivara, Kritas, Bourriel, Alexopoulos, & Kyriakis, 2007). On the other hand, increasing host immunity or resistance to the pathogens of swine respiratory disease would be an alternative important approach to control disease development. Recently, several pilot studies have been conducted to explore the genetic mechanisms of susceptibility/resistance to swine respiratory disease (Fang et al., 2013(Fang et al., , 2015Huang et al., 2017;Okamura et al., 2012). Using quantitative genetic methods, five significant and 18 suggestive quantitative trait loci (QTL) for swine respiratory disease traits and immune capacity traits were identified in a Landrace pig line selected for SRD resistance over five generations (Okamura et al., 2012). Genome-wide association studies for an SRD trait have been performed in Chinese Erhualian pigs, and five novel loci were identified on swine chromosomes 2, 8, 12, and 14 . The CXCL6, CXCL8, KIT, and CTBP2 genes were identified as possibly corresponding to host resistance or susceptibility to the pathogens of swine respiratory disease . Some other candidate genes, such as TLR4 (Fang et al., 2013) and CYP1A1 (Fang et al., 2016), have also been reported to be relevant to host resistance or susceptibility to swine respiratory disease. In addition, the primary immunodeficiency, Toll-like receptor signaling, and steroid metabolism pathways might play important roles in the regulation of the inflammatory response to M. hyopneumoniae infection in pigs (Fang et al., 2015). Previous genetic studies have provided important clues to elucidate the biological mechanism of host resistance to the pathogens of respiratory disease in pigs, which have suggested that the genetic basis of host resistance to swine respiratory disease is complex.
Bamaxiang pigs are a well-known Chinese indigenous breed distributed around Bama Yao Autonomous County in Guangxi Province, China. These pigs are characterized by two-end black coat color, small body size, and good meat quality (Wang et al., 2011). Owing to their small body size, Bamaxiang pigs can serve as model animals for biomedical research (Liu, Zeng, Shang, Cen, & Wei, 2008). In the present study, Bamaxiang pigs were used to investigate host resistance or susceptibility to swine respiratory disease. We recorded body weights with different ages and two traits representing the extent to which the pigs suffered from swine respiratory disease in a herd of Bamaxiang pigs under consistent indoor farming conditions.
According to the effects of both SRD traits on weight gain, which were estimated via the maximum likelihood method, the pigs were divided into disease-susceptible and disease-resistant groups. Then, the polygenic effects on animal adaptation to the respiratory disease were explored through genetic differentiation analyses with a whole-genome high-density SNP dataset. Finally, a transcriptomics study was conducted to compare the expression of candidate genes between healthy lung tissues and lung lesion tissues.
The present study attempted to employ population genetics and transcriptomics methods to elucidate the genetic mechanism of host adaption to a chronic disease in mammals. Our findings provide novel insights into the genetic architecture of host adaptation to swine respiratory disease, will be beneficial for the breeding of pigs resistant to swine respiratory disease in the swine industry, and may provide a useful reference for revealing the molecular genetic basis of host resistance to chronic refractory diseases in other mammals.

| Ethics statement
All procedures performed in this study involving animals were in compliance with guidelines for the care and use of experimental animals established by the Ministry of Agriculture of China. The ethics committee of Jiangxi Agricultural University specifically approved this study.

| Experimental animals
The experimental animals investigated in the present study were bought from the core-breeding farm of Bamaxiang pigs in Bama Yao Autonomous County, Guangxi Province. These animals consisted of 155 males and 159 females from 43 sire and 61 dam families, which were selected so that the genetic diversity of the population was as high as possible in this study. These piglets were immunized with several common vaccines according to their instructions, including those for classical swine fever, foot-and-mouth disease, pseudorabies, porcine reproductive and respiratory syndrome, and porcine circovirus type 2. All male pigs were castrated postweaning. At the age of approximately 2 months, the pigs were transported to Nanchang, Jiangxi Province, and were raised under consistent indoor rearing conditions on a farm in two rows of head-to-head pens, where each pen (approximately 10 m 2 ) hosted 5-8 pigs. All animals were fed twice a day (10:00 a.m. and 16:00 p.m.) with a diet containing 16% crude protein, 3,100 kJ of digestible energy, and 0.78% lysine. Water was available ad libitum from automatic faucets. During the study period, the affected pigs with severe respiratory symptoms were treated with the medicine florfenicol following the provided instructions. All pigs were uniformly slaughtered at 300 ± 3 days at an abattoir for trait measurements.

| Phenotypic measurements
To estimate the extent to which the pigs suffered from swine respiratory disease, we recorded two related traits. The first trait was the enzootic pneumonia-like (EPL) score, and the other was lung lesion (LL) levels.
The EPL score was described previously .
In brief, the respiratory status of all pigs was recorded twice a day within one hour of the feeding time. When we found a pig standing with a dry cough or lying down with obvious abdominal breathing, accompanied by fast breath and a lack of appetite, an EPL score of one was given to the individual. The sum of the EPL scores was then calculated for each individual during the rearing period from 100 to 300 days of age, which was treated as the value of the EPL score trait for that pig. In the Bamaxiang pig population, the smallest EPL score was zero, indicating that the animal did not suffer from asthma; the largest EPL score was 112, indicating that the pig suffered from the most severe asthma.
Immediately after slaughter, each pig lung was placed in a brightly lit position, and photographs of both the anterior and posterior sides were taken. Lung lesion levels were estimated based on lung photographs referring to a scoring criterion similar to that described in a previous report (Qu, Liu, Yao, & Jin, 2012). First, we evaluated the proportion of the lesion area in different parts of pig lung, including the anterior and posterior sides. Then, different weights were assigned to different lung parts; that is, weights of 10%, 10%, 25%, 10%, 10%, 25%, and 10% were assigned to the left apical, left cardiac, left diaphragmatic, right apical, right cardiac, right diaphragmatic, and accessory lobes, respectively. As all six parts of the lung except for the accessory lobe can be observed from both the anterior and posterior sides, a weight of 50% was equally assigned to the anterior and posterior sides, and a weight of 100% was assigned to the posterior side of the accessory lobe. Finally, the lung lesion level in each pig was obtained by summing the lesion area proportions of these seven lobes multiplied by their own weights (Supporting information Figure S1).
All the Bamaxiang pigs were weighed using an electronic scale, and their weights were recorded at the ages of 150, 180, 240, and 300 days.

| SNP genotyping and quality control
Genomic DNA was extracted from the ear tissue of each animal using a standard phenol/chloroform approach. DNA quality and concentrations were determined via agarose gel electrophoresis and using a NanoDrop 1,000 spectrophotometer (Thermo Fisher, USA). Before genotyping, DNA samples were diluted to a final concentration of 50 ng/ml. All experimental pigs were genotyped using Affymetrix Axiom customized genotyping arrays containing 1,348,804 porcine SNPs (Zhang et al., 2017). Quality control was performed with Plink v1.07 (Purcell et al., 2007) with parameters of an SNP call rate>90% and a minor allele frequency (MAF) > 0.01. After the quality control step, a total of 782,364 SNPs were used for genetic differentiation analyses, among which 25,427 SNPs were located on the X chromosome.

| Animal grouping
To perform genetic differentiation analyses, the Bamaxiang pigs were divided into disease-susceptible and disease-resistant subgroups based on their phenotypic data related to respiratory disease and the effects of respiratory disease on their weight gain. Before grouping the animals, we conducted independent maximum likelihood tests to estimate the reasonable effects of the EPL score and LL levels on porcine weight gain. The maximum likelihood test was described previously . Taking the EPL score as an example, we first determined a threshold of EPL score (x): If the EPL score of an individual was less than or equal to x, the individual was assigned to the disease-resistant group; otherwise, it was assigned to the disease-susceptible group. t test values between the two groups were calculated for body weight on all test days. Then, the likelihood statistic L(x) was calculated according to the following model: where p(x) denotes the p-value of the t test between the two groups at a given presumed x, and i ranged from 1 to 4, indicating the time point (150,180,240, and 300 days, respectively). The most likely classification threshold of the EPL score was determined based on the maximum likelihood statistic L(x). Then, the pigs were assigned to the disease-susceptible group if their EPL scores were greater than the most likely classification threshold but were assigned to the disease-resistant group if their EPL scores were less than or equal to the threshold.
The same strategy was used for animal grouping based on the effects of LL levels on porcine weight gain. Additionally, in the third grouping method, we clustered the animals in both disease-resistant groups according to EPL scores and LL levels in the disease-resistant group, while the animals in both disease-susceptible groups according to EPL scores and LL levels were assigned to the disease-susceptible group.

| Genetic differentiation analyses
Two parameters, population genetic differentiation and the delta allele frequency, were estimated between the disease-susceptible and disease-resistant groups in the Bamaxiang pig population. Unbiased genetic differentiation estimates of the fixation index (F ST ) were calculated as described by Akey, Zhang, Zhang, Jin, and Shriver (2002) using the SNP dataset. Briefly, F ST was estimated as follows: where MSG and MSP denote the observed mean square errors for loci within and between populations, respectively, and n c is the average sample size across samples, which incorporates and corrects for the variance in the sample size over population.
In the above formulae, n i denotes the sample size in the ith population; p Ai is the frequency of SNP allele A in the ith population; and − p A is a weighted average of p A across populations. Because the range of F ST was originally defined as 0 to 1 (Wright, 1951), negative F ST values with no biological interpretation were set to 0. The mean F ST value in each 50 kb window, with a sliding window size of 25 kb, was calculated to represent the genetic differentiation extent of a locus.
Finally, the mean F ST values of 93,425 loci in total were obtained.
The top 1% of loci according to genetic differentiation values were defined as candidate regions related to host resistance or susceptibility to swine respiratory disease.
The delta allele frequency is the difference in the frequencies of an allele between the two populations. To verify the results of the population differentiation analysis, we calculated the absolute delta minor allele frequency (deltaMAF) of each SNP within the candidate loci between the disease-susceptible and disease-resistant subgroups, treating the disease-susceptible subgroup as the reference population. The mean deltaMAF values were calculated in each 50 kb window, with a sliding window size of 25 kb, in the same manner as for the F ST computation.

| Permutation test
The permutation test was performed similar to a previous report (Churchill & Doerge, 1994). Individuals in the experiment were indexed from 1 to n. The data were shuffled by computing a random permutation of the indices 1,..., n and assigning the ith phenotypic value to the individual whose index was given by the ith element of the permutation. The shuffled data were then analyzed for F ST calculation according to the initial grouping state. The resulting test statistics at each analysis point are stored, and the entire procedure of shuffling and F ST calculation was repeated 1,000-times. At the end of this process, we stored the results of genetic differentiation analyses on 1,000 shuffled datasets. For each time, the 100(1-α) percentile was found. The mean of these 100(1-α) percentile values was set as the threshold of significance level α. For example, the mean of the 99.5 percentile values of 1,000-times permutation tests was set as the threshold of p = 0.005.

| Characterization of candidate genes
Candidate genes within the above candidate regions were searched in Build 10.2 of the pig reference genome via the Ensembl Genome Browser (http://www.ensembl.org/Sus_scrofa/Info/Index). The genes that were orthologous to human genes were annotated with Ensembl Biomart (https://www.ensembl.org/biomart) with filter parameters of "orgholog confidence==1" and "ortholog_one2one." GO terms and KEGG pathways were enriched using the default options of the ClueGO plug-in (Bindea et al., 2009) of the Cytoscape 3.5.0 platform (Shannon et al., 2003), setting human orthologous EnsemblGeneIDs as the gene list and the human Gene Ontology database as the query database.
Additionally, we input these human orthologous genes to query the STRING 10.5 database (Szklarczyk et al., 2015), searching protein-protein interaction (PPI) networks related to pig adaption to swine respiratory disease. The topological properties of PPI networks, such as betweenness, node degree, and eigenvector, were analyzed using the CentiScape v2.2 plug-in (Scardoni, Petterlini, & Laudanna, 2009) of the Cytoscape 3.5.0 platform. The intersecting genes passing the thresholds of the three parameters were considered to be critical nodes. These nodes with high centrality are also referred to as key node genes and act as key connector proteins of the major network processes. (2)

| Gene expression analysis
Five infected tissue samples from the lung lesions and nine noninfected tissue samples from healthy lung parts were collected from eleven pigs in a commercial black pig population (Supporting information Table S1) under common indoor farming conditions. RNA sequencing data were obtained from these fourteen lung tissue samples; RNA-seq reads from each sample were aligned to Build 10.2 of the pig reference genome with STAR v.2.5.3 (Dobin et al., 2013), using genomic sequence and transcript annotations derived from Release 89 of the Ensembl database. Transcripts were assembled with Cufflinks v.2.1.1 (Trapnell et al., 2010). Gene counts were calculated with FeatureCounts implemented in Subread v.1.5.2 (Liao, Smyth, & Shi, 2014) and were then normalized and used for differential gene expression analysis between infected and noninfected lung samples with the DESeq2 package (Love, Huber, & Anders, 2014). The P-values from the differential gene expression analyses were adjusted via Benjamini-Hochberg multiple corrections.
An adjusted p-value <0.05 was set as the level of significance.

| Phenotypic statistics of EPL scores and LL levels
Although the Bamaxiang pigs were immunized with several common vaccines, including the classical swine fever, foot-and-mouth disease, pseudorabies, porcine reproductive and respiratory syndrome, and porcine circovirus type 2 vaccines, this herd suffered from swine respiratory disease. Similar to the previous study of Erhualian herd , no typical clinical symptoms of these viral diseases were observed during the entire rearing time on the farm, and the disease-tolerant period was long, lasting for almost the entire farming stage in Nanchang, and the hot, wet environment could intensify the incidence of respiratory disease (Supporting information Figure S2). We again defined this respiratory disease infecting the Bamaxiang herd as "enzootic pneumonia-like (EPL) disease".
Among the 314 Bamaxiang individuals, only 4.8% (15 individuals) did not appear to exhibit any respiratory symptoms, while 50.3% suffered from severe dyspnea and coughs at least 16 times (Figure 1a).
The average EPL score of the Bamaxiang herd was 22.3 (Supporting information Table S2) Table S2). Surprisingly, the phenotypic correlation between EPL scores and LL levels was small (r 2 = 0.152, p = 6.79 × 10 -3 ) in the Bamaxiang pig herd (Figure 1d). For example, some pigs with zero EPL scores suffered from severe lung lesions, presenting an LL level of more than 0.50, while some pigs with no lung lesions suffered from severe dyspnea. These findings suggested that there was not a strong relationship between the two traits of swine respiratory disease.
To clarify whether the relationship between EPL scores and LL levels was constant throughout the observed growth period, we divided the growth period into six intervals: age ≤120 days; 120 days <age ≤ 150 days; 150 days <age ≤ 180 days; 180 days <age ≤ 210 days; 210 days <age ≤ 240 days; and age ≥240 days. We performed correlation analysis between EPL scores and LL levels in each interval. We found that there were no correlations between EPL scores and LL levels in the two age intervals of less than or equal to 150 days. The highest correlation (r 2 = 0.279, p = 5.13 × 10 -7 ) existed within the interval of 180 days <age ≤ 210 days, while those for the other three intervals were approximately 0.117 and were all less than the overall correlation value of 0.152 (Supporting information Figure S3). These findings indicated that the clinical respiratory symptoms presented a stronger relationship with the extent of lung lesions during the growth period from 180 to 210 days than in other growth stages in the Bamaxiang pigs.

| Effects of EPL scores and LL levels on weight gain
It has been repeatedly reported that swine respiratory disease affects weight gain in pig farming (Hoy, 1994;Huang et al., 2017). In the present study, we also investigated the effect of respiratory disease on body weights in the Bamaxiang population. We found that both EPL scores and LL levels showed a negative correlation with body weight at the ages of 150, 180, 240, and 300 days (Table 1), and the EPL score generally had a larger negative effect on body weight than LL levels at the above four ages.
To perform genetic differentiation analyses in this Bamaxiang population, we applied three strategies for dividing the herd into two groups, in which the pigs were divided into disease-susceptible and disease-resistant groups based on EPL scores or LL level, or both EPL scores and LL levels (Table 1). In the grouping method based on the EPL score, the likelihood statistic was maximized when the threshold of the EPL score (x) was set to 4 ( Figure 1a). We compared the body weights of the disease-susceptible pigs (EPL score>4) and disease-resistant pigs (EPL score ≤4). Significant differences were found between the two groups at the ages of 150, 180, and 240 days (p < 0.01), but not at the age of 300 days (p = 0.08 likelihood statistic was maximized when the threshold of the LL level (y) was set to 0.204 (Figure 1c). Significant differences were found between the two groups at days 180, 240, and 300 (p < 0.01), but not at day 150 (p = 0.11). In the grouping method based on both EPL scores and LL levels (Figure 1d), significant differences were found between the two groups at all ages (150, 180, 240, and 300 days).
The body weights of the disease-susceptible group decreased by 16.7% (150 days), 18.7% (180 days), 11.9% (240 days), and 6.4% (300 days). These findings indicated that this respiratory disease has a persistent influence on the weight gain of disease-susceptible pigs during the growth period, and the combined effect of EPL scores and LL levels on weight gain was greater than each single effect.

| Candidate loci detected through genetic differentiation analyses
We estimated the genetic differentiation of the fixation index (F ST ) within these three-paired animal groups for 93,425 loci throughout the genome (Supporting information Figure S4). Under each grouping method, the F ST distribution was highly skewed toward zero and exhibited a long tail (Figure 2a). In the grouping method based on the EPL score, the maximum F ST value was 0.1109, and the mean was 0.0051; based on LL levels, the maximum F ST value was 0.0354, and the mean was 0.0015; and based on both EPL scores and LL levels, To validate the reliability of the candidate loci we identified, we also examined the mean deltaMAF in the candidate regions. We

| Overview of candidate genes for host adaptation to SRD
A total of 459, 456, and 527 candidate genes (Supporting information  Figure S4). Of the above 459, 456, and 527 candidate genes, 314 were commonly detected in the EPL and EPL_LL grouping methods; 116 were shared between the EPL and LL grouping methods; 170 were shared between the LL and EPL_LL grouping methods; and 107 were shared among all three grouping methods ( Figure 3a). Using these 107 common candidate genes as an input gene list, 28 significant GO terms and two significant KEGG pathways were enriched (Supporting information Table S4), which were clustered into ten GO groups (Figure 3b). To investigate the interaction networks of candidate genes related to host adaptation to SRD, we used the protein-protein interaction (PPI) database from STRING (Szklarczyk et al., 2015) to construct undirected graphs of gene coding proteins obtained from the three different grouping methods. We detected 295 nodes and 528 edges of interactions in the protein network for the F I G U R E 3 The 107 common candidate genes and their enriched gene ontology terms and KEGG pathways. (a) Venn diagram showing that a total of 107 candidate genes were common to the top 1% of loci obtained under the three different grouping methods. A list of these 107 genes is provided in Supporting information Table S7. (b) Gene ontology terms and KEGG pathways were enriched using these 107 common genes. The enriched terms are color-coded to reflect relatedness in the ontology or functional proximity. All GO and KEGG terms are listed in Supporting information Table S4 EPL score grouping methods using 459 candidate gene coding pro-  Figure S5C). Taking all of these findings together, we observed the intersection of the results of the three different grouping methods and ultimately obtained seven common genes (TGFB1, SNRPA, TRAF6, CYP2A6, CD44, CYP2B6, and PLAUR) with high centrality values (Supporting information Figure S5D).
We further investigated the SNP effects of these common genes on both SRD-related traits of EPL scores and LL levels. First, we examined the genotypes of the SNPs with largest F ST values located in the candidate loci harboring these common genes, which were defined here as the genotypes of the corresponding genes. Five candidate genes (TGFB1, SNRPA, TRAF6, CD44, and CYP2B6) exhibited only two genotypes (Supporting information Table S5). Association analysis revealed significant differences in EPL scores between the two genotypes of TGFB1, SNRPA, TRAF6, CYP2B6, and CD44, with p-values of 1.01 × 10 -2 , 1.13 × 10 -3 , 4.35 × 10 -4 , 7.37 × 10 -4 , and 1.67 × 10 -4 , respectively. Significant differences in LL levels also existed between the two genotypes of TGFB1, TRAF6, CYP2B6, and CD44, with p-values of 7.48 × 10 -3 , 2.57 × 10 -3 , 2.41 × 10 -2 , and 1.37 × 10 -2 , respectively. In addition, we found another gene coding for a family member of the sialic acid binding lectins among 107 common candidate genes, the CD22 gene, showing significant differences between the EPL scores of its three genotypes (p = 1.28 × 10 -3 ). When the combined effects of CD22 and CD44 were considered, a more significant difference (p = 4.78 × 10 -5 ) in EPL scores was found among the individuals with different genotypes (Supporting information Table S6) than was observed for the independent effects of CD22 and CD44.

| Differential gene expression in the lungs
To gain insight into the mechanism whereby pig lungs respond to SRD infection, we generated and analyzed transcriptomic data from nine healthy lung tissues and five affected lung lesions collected from a commercial black pig population under common indoor farming conditions. Hierarchical cluster analyses showed that the transcriptomic profiles of the healthy and affected lung samples differed from each other, and the healthy and affected lung samples were grouped on their own branches (Supporting information Figure S6A). A total of 23,382 Sus scrofa reference genes downloaded from the Ensembl database (ftp://ftp.ensembl.org/pub/release-75/gtf/sus_scrofa/) and 20,564 predicted genes generated with Cufflinks software were used for the differential expression analysis. Among these genes, 5,894 reference genes (2,716 upregulated, 3,178 downregulated) and 2,137 predicted genes (1,416 upregulated, 721 downregulated) were detected as being differentially expressed between the infected and noninfected lung tissues. In summary, 4,132 and 3,899 genes were up-and downregulated in healthy tissues compared with infected tissues, respectively (Supporting information Figure S6B). Notably, among the above 107 common candidate genes shared by the three grouping methods and 158 key node genes highlighted through PPI analysis, 72 genes were identified as differentially expressed genes (Supporting information Table S7) Table S7).

| The Bamaxiang pigs were under selection pressure in the SRD environment
In the present study, we defined an SRD environment as the biological and physical factors, along with their interactions, that cause pigs to develop swine respiratory disease. We consider pigs to always be in an SRD environment under common indoor farming conditions because of the major reasons indicated below. to 50 kg, or <20 kg, respectively (Chinese Academy of Agricultural Sciences, 1958). Contagious respiratory disease is one of the common diseases plaguing local Chinese pigs during rearing. Third, we observed that most of the Bamaxiang pigs (93.2%) in our study suffered from asthma and that lesions occurred in 98.4% of their lungs.
Respiratory symptoms appeared throughout the farming period and were most severe at an age of ~160 days, during the hottest and wettest season (June and July), at which time the highest frequency of cough occurred in this population.
Furthermore, we believe that human-driven selection pressure and swine respiratory disease have shaped the evolutionary trajectory of these Bamaxiang pigs. In the breeding of local Chinese pigs, breeders prefer to choose healthy, normal, or fast-growing pigs to reproduce subsequent generations. Disease-resistant pigs would have a greater opportunity to be selected as sires or sows to generate offspring. In addition, during the rearing period, some pigs suffering from severe lung lesions will be easy to die if simultaneously infected by secondary or opportunistic pathogens. In the present study, some Bamaxiang pigs died of swine respiratory disease in the farming process. Unfortunately, these pigs were not included in this study because they had no available genotypic data and no phenotypic records.
Considering these findings together, we believe that the Bamaxiang pigs were under selection pressure in the SRD environment. These pigs generally experienced a long-term mutual struggle with respiratory disease pathogens under common indoor feeding conditions. In the SRD environment, these pigs could be genetically classified into two subgroups: one resistant to SRD and another susceptible to swine respiratory disease. The related genes may have differentiated during pig adaptation to respiratory disease, and their allele frequencies might also have changed.

| SRD phenotypes and their effects on weight gain
Previously, there were two major methods for determining the extent of swine respiratory disease through genetic analyses: one based on the count or frequency of asthma cases  and the other based on the scoring of lung lesions (Okamura et al., 2012).
Here, we employed both methods to estimate the extent of SRD suffering in a Chinese Bamaxiang pig herd. Compared to previous findings in Erhualian pigs, the EPL scores of the Bamaxiang pigs (average value: 22.3) were generally higher than those of Erhualian pigs (average value: 3.5), possibly due to differences in germplasms and/or statistical batches, or other factors. This observation suggested that Bamaxiang pigs might be more susceptible to swine respiratory disease than Erhualian pigs under common indoor farming conditions.
Unexpectedly, we found that the phenotypic correlation between EPL scores and LL levels was as low as 0.152 in the Bamaxiang pig population, which indicated that these two traits might reflect different impacts of swine respiratory disease on pigs. Similarly, a low correlation between clinical signs of respiratory system disorders and lung lesions at slaughter was previously observed in veal calves (Leruste et al., 2012). Both traits exerted a negative effect on the weight gain of Bamaxiang pigs, which was consistent with previous reports from other pig breeds (Fang et al., 2015;Huang et al., 2017).
We considered the EPL scores to exert a more adverse effect on weight gain than LL levels, due to greater differences between the disease-resistant and disease-susceptible pigs when grouped by EPL scores. Furthermore, the superposed effect of EPL scores and LL levels together was larger than each single effect alone.

| Candidate genes for host adaptation to SRD
Under general indoor farming conditions, pigs usually live in a microbe-rich environment, and interactions will occur between the F I G U R E 4 Expression of the six highlighted genes in healthy lung tissues and affected lung lesions. The differences in the expression of TRAF6, CD44, CD22, and SNRPA between healthy lung tissues and affected lung lesions were extremely significant, with corrected P-values of 2.09 × 10 -5 , 5.24 × 10 -15 , 4.71 × 10 -10 , and 6.48 × 10 -5 , respectively. TGFB1 and CYP2B6 were significantly differentially expressed in healthy lung tissues and affected lung lesions, with corrected p-values of 1.89 × 10 -2 and 4.88 × 10 -2 , respectively host and microbes. The host immune system will be triggered to cope with the invasion of microorganisms, and the microbes in turn influence the host, such as through adverse effects on host growth due to SRD infection. In this study, we used genetic differentiation analyses to identify candidate genes responding to host adaptation to SRD. We identified 107 common candidate genes using three animal grouping methods based on EPL scores, LL levels, and both parameters. Interestingly, two major groups of biological processes, related to host immune responses and fatty acid metabolism, were enriched, which was consistent with the interaction between the host and microbes. These results also indicated that the genes involved in these biological processes might play important roles in host adaptation to SRD.
For host immune responses, two enriched biological processes, myeloid dendritic cell differentiation (GO:0043011) and myeloid dendritic cell activation (GO:0001773), overrepresented by RELB, TGFB1, and TRAF6, might play important roles in the initiation of immunity in response to swine respiratory disease, because it has been reported that myeloid dendritic cells pick up antigens at the periphery and travel to T-cell areas to initiate immunity (Steinman & Inaba, 1999). Additionally, the biological processes of cytokine production and regulation of cytokine production involved in the immune response were enriched by the overrepresented genes CLC, FFAR3, TGFB1, and TRAF6. For fatty acid metabolism, six related GO terms were enriched by the cytochrome P450 family 2 genes CYP2A6, CYP2A7, CYP2A13, CYP2B6, CYP2F1, and CYP2S1; these terms included the epoxygenase P450 pathway, arachidonic acid metabolic process, unsaturated fatty acid metabolic process, long-chain fatty acid metabolic process, icosanoid metabolic process, and fatty acid derivative metabolic process. In addition to their involvement in fatty acid metabolism, cytochrome P450 (CYP) enzymes often play a dominant role in target tissue metabolic activation of xenobiotic compounds in the respiratory and gastrointestinal tracts (Ding & Kaminsky, 2003). Furthermore, the enriched biological process of regulation of phospholipid metabolic process, overrepresented by APOC1, APOC2, and TGFB1, is related to the host immune response to microbial infection, as phospholipids have proven to be essential in inflammation (Nixon, 2009 TRAF6 encodes tumor necrosis factor receptor (TNFR)-associated factor protein 6, which functions as a signal transducer in the NF-kappa B pathway that activates I kappa B kinase in response to proinflammatory cytokines (Cao, Xiong, Takeuchi, Kurama, & Goeddel, 1996). TRAF6-mediated signals have proven critical for the development, homeostasis, and/or activation of B cells, T cells, and myeloid cells, as well as for organogenesis of thymic and secondary lymphoid tissues (Walsh, Lee, & Choi, 2015). In addition, TRAF6 was demonstrated to be related to the pathways of bacterial infections in cystic fibrosis airways (Walsh et al., 2001). CD44 is a transmembrane adhesion receptor and serves as the major cell-surface receptor for the nonsulfated glycosaminoglycan hyaluronan (Aruffo, Stamenkovic, Melnick, Underhill, & Seed, 1990).
CD44 is involved in recruiting T cells to inflammatory sites and reg- ulates T-cell-mediated endothelial injury (Degrendele, Kosfiszer, Estess, & Siegelman, 1997;Rafi-Janajreh et al., 1999). Teder et al. (2002) demonstrated a role of CD44 in resolving the inflammatory response following lung injury. CD22 is also a transmembrane sialo-adhesion protein expressed by nearly all mature B cells as well as a member of the immunoglobulin superfamily (Engel, Wagner, Miller, & Tedder, 1995). In humans, CD22 has been observed to positively associate with regional emphysema severity of the lungs and is regarded as a potential causal gene for airflow obstruction (Lamontagne et al., 2014). TGFB1 encodes a secreted ligand of the transforming growth factor-beta superfamily of proteins. The TGFB1 gene product is a secreted protein with numerous functions, including roles in lung development and disease (Bartram & Speer, 2004); TGFB1 appears to be the predominant isoform involved in pulmonary fibrosis (Bartram & Speer, 2004). In mice, an overdose of TGFB1 transmission resulted in disruption of epithelial differentiation and inhibitory synthesis of the Clara cell secretory protein, phospholipids, and surface proteins A, B, and C. In humans, the T allele of the SNP rs1800469 in the promoter region of the TGFB1 gene has been associated with phenotypes related to asthma (Silverman et al., 2004), and TGFB1 levels in the airways of asthmatics are reported to be higher than in normal subjects (Nakao, 2001).
Similarly, the expression of TGFB1 in affected pig lungs was found to be significantly higher than that in healthy lungs in the present study. SNRPA encodes small nuclear ribonucleoprotein polypeptide A, which has been reported to function in the binding stem loop II of the U1 snRNA and might be involved in coupled prem-RNA splicing and polyadenylation (Jessen, Nagai, Oubridge, Teo, & Pritchard, 1991). Recently, many reports have indicated that SNRPA is highly associated with tumors and apoptosis (Tewari, Beidler, & Dixit, 1995). CYP2B6 encodes a member of the cytochrome P450 superfamily of enzymes. The cytochrome P450 proteins are monooxygenases that catalyze many reactions involved in drug metabolism and the synthesis of cholesterol, steroids, and other lipids.
Altogether, these six genes could be considered as strong candidates associated with pig adaptation to SRD.

| CON CLUS IONS
In this study, we first investigated the relationship of two SRD-related traits, enzootic pneumonia-like (EPL) scores and lung lesion (LL) levels in a Bamaxiang population. We found that there was no strong relationship between EPL scores and LL levels in the Bamaxiang population. Swine respiratory disease has a persistent influence on weight gain in disease-susceptible pigs during the growth period, and the joint effect of EPL scores and LL levels on weight gain was found to be greater than each single effect alone.
By performing genetic differentiation analyses, we identified 107 common candidate genes, which are mainly involved in biological processes related to the immune response, fatty acid metabolism, lipid metabolism, and growth factor signaling pathways. Among these candidate genes, six genes (TRAF6, CD44, CD22, TGFB1, CYP2B6, and SNRPA) were highlighted as the most essential genes related to pig adaptation to SRD due to their central regulatory roles in the host immune response or fat metabolism and their differential expression between healthy lung tissues and lung lesion tissues.
In the present study, we attempted to reveal the molecular mechanism of mammalian adaption to a chronic disease through genetic differentiation and RNA sequencing analyses. Our results advance the understanding of the genetic mechanisms of host resistance or susceptibility to swine respiratory disease.

ACK N OWLED G EM ENTS
This study was funded by grants from National Swine Industry and

CO N FLI C T O F I NTE R E S T
None Declared.

DATA ACCE SS I B I LIT Y
RNA-seq data for this study are available at the NCBI database with the access number PRJNA496272.