Genetic Features of Plasmid- and Chromosome-Mediated mcr-1 in Escherichia coli Isolates From Animal Organs With Lesions

The genomic context of the mcr-1 gene in Escherichia coli from animal feces has been widely reported. However, less is known about the mcr-1-carrying plasmid characteristics and other functional regions of Escherichia coli isolates from animal organs with lesions. The present study investigated the antimicrobial resistance, population structure, and genetic features of mcr-1-positive Escherichia coli strains isolated from animal organs with lesions. The antimicrobial susceptibility testing indicated that 24 mcr-1-positive Escherichia coli isolates were resistant to at least three or all antimicrobial categories. MLST analysis suggested that the dominant clone complexes (CC) were mainly CC156, CC448, and CC10. In addition, ST10596, a newly discovered sequence type in swine, failed to be classified. Meanwhile, the mcr-1 gene located on the different plasmids was successfully transferred to the recipients, and whole-genome sequencing indicated the mcr-1 gene was embedded in mcr-1-pap2 cassette but not flanked by ISApl1. The mcr-1 gene is located on the chromosome and embedded in Tn6330. Furthermore, NDM-5 was found on the IncX3-type plasmid of J-8. The virB6 and traI gene of type IV secretion system (T4SS) were truncated by IS2 and IS100 and located on the IncX4- and the IncHI2/HI2A/N-type plasmids, respectively. The multidrug-resistant (MDR) region of IncHI2/HI2A/N-type plasmids contained two class 1 integrons (In0, In640) and four composite transposons (Tn4352, Tn6010, cn_4692_IS26, cn_6354_IS26). Overall, 24 mcr-1-positive Escherichia coli isolates in our study showed MDR, or even extensively drug resistant (XDR), and exhibited population diversity. The T4SS gene truncation by the insertion sequence may affect the efficiency of plasmid conjugative transfer. Furthermore, the class 1 integrons and composite transposons in the MDR region of IncHI2/HI2A/n-type plasmid contributed to the multireplicon plasmid formation, the acquisition, and transfer of antimicrobial resistance genes (ARGs).


INTRODUCTION
The overuse and misuse of antimicrobials have led to the rapid development of resistance in Gram-negative bacteria, and this has now become a serious public health crisis (Martin et al., 2020). Among these bacteria, Escherichia coli is the most common pathogen which causes animal and human diseases on a global scale while being a major cause of multiple acquired infections (Aslam et al., 2018). Over the past few years, the emergence and common occurrence of extended-spectrum beta-lactamase (ESBL) -and carbapenemase-producing E. coli have result in reduced effectiveness of antimicrobial treatment (Hordijk et al., 2020), such that doctors are now considering alternative antimicrobials, such as colistin.
For many years, colistin (polymyxin E) has been used in animal production for the treatment and prevention of infectious diseases. It was subsequently identified as a form of last-resort treatment against severe human infections caused by multidrugand carbapenem-resistant Gram-negative bacteria (Kempf et al., 2016;Poirel et al., 2017). Colistin resistance was previously thought to be chromosome mediated, which mainly related to the point mutations of two-component systems (TCSs) PhoP/PhoQ and PmrA/PmrB and can transmit colistin resistance through vertical transmission (Olaitan et al., 2014). It was only in 2015 that the plasmid-mediated mcr-1 gene was first described in E. coli and Klebsiella pneumoniae. Both species were, respectively, isolated from pigs and humans in China and showed the horizontal transfer of colistin resistance . The emergence of this gene has further reduced the effectiveness of the alternative antimicrobials, thereby posing new challenges to public health. Since then, colistin-resistant mcr-1-positive E. coli isolates have been reported in the environment, humans, animals, and food worldwide (Falgenhauer et al., 2016;Luo et al., 2017;Shen et al., 2018). Numerous researches have shown that the high clonal diversity was observed among mcr-1-positive E. coli (MPECs) and mainly belonged to the common dominant clone complexes, such as ST10, ST48, ST165, etc. (Matamoros, 2017;Bai et al., 2018). Among these, ST10 is the most common transmission clone group and frequently disseminates through different means (Sellera et al., 2017;Patiño-Navarrete et al., 2020).
The mcr-1 gene is found in multiple plasmid types, including IncI2, IncX4, IncX3-X4, IncHI1, IncHI2, IncP, IncY, IncK2, IncFI, IncFI, IncFII, and IncFIB-I2 (Malhotra-Kumar et al., 2016;McGann et al., 2016;Poirel et al., 2016;Sun et al., 2016;Xavier et al., 2016;Zurfluh et al., 2016;Donà et al., 2017;Zhang et al., 2017). This was consistent with multiple mcr-1 mobilization events which tend to potentially facilitate the association of mcr-1 with other resistance mechanisms, thereby creating multidrug-resistant bacterial hosts (Zhong et al., 2018). Of these plasmids, IncI2-, IncX4-, and IncHI2-type ones are the most common vectors of mcr-1, with the insertion sequence (IS) ISApl1 mostly appearing upor downstream of the gene (Matamoros, 2017). Few studies reported mcr-1-carrying hybrid plasmid and their mechanism of resistance. The hybrid plasmid IncX3-X4 is based on homologous recombination or the transposase-mediated formation of a cointegrate between IS26 and the nic site of oriT. This also suggested that even if the insert sequence on the plasmid carrying the mcr-1 gene does not mobilize the transfer of mcr-1, it can still be recombined with other plasmids containing the same insert sequence to form a more complex hybrid plasmid . These researches indicated that IS may contribute to plasmid-plasmid and plasmid-chromosome rearrangements as well as gene mobilization. Therefore, it might be significantly valuable to study IS on plasmids carrying the mcr-1 gene.
This study aimed to investigate the population structure of MPECs and the features of mcr-1-carrying plasmids in MPECs from animals which showed reduced effectiveness to antimicrobial. The present study may improve the understanding of the epidemiology, likelihood of spread and the potential, and evolutionary mechanisms of the mcr-1 gene in veterinary medicine.

Isolation of mcr-1-Positive E. coli Strains
Samples were obtained from lesion organs of swine, poultry (chicken), and bovine (dairy cow) between 2016 and 2019 in Yangling, Shaanxi, China. A total of 109 samples were collected from the organs with lesions in swine (n = 37) and chickens (n = 53) suffering from colibacillosis and the uterine rinsing fluid of dairy cows (n = 19) suffering from endometritis. The samples were collected in different farms, and each sample was from a different animal. The samples were stored aseptically at 4 • C and transported to the laboratory in time for isolation and identification of strains. All animal procedures and study design were approved by the animal ethics committee of Northwest A&F University (Approval No.: 2016012).
All samples were cultivated on MacConkey agar (MAC) plates containing 2 µg/ml colistin and incubated at 37 • C for 18 h. The species were further identified by the amplification and sequencing of 16S rRNA. Then, MPECs were detected by amplification and sequencing of mcr-1 .

Multilocus Sequence Typing of MPECs
MPECs were sequenced using Illumina HiSeq TM2000 platform at the Novogene Bioinformatics (Beijing, China). For Illumina sequencing, at least 1 µg genomic DNA was used for each strain in sequencing library construction. DNA samples were sheared into 400-500 bp fragments using a Covaris M220 Focused Acoustic Shearer following manufacture's protocol. Illumina sequencing libraries were prepared from the sheared fragments using the NextFlex TM Rapid DNA-Seq Kit. The prepared libraries then were used for paired-end Illumina sequencing (2 × 150 bp) on an Illumina HiSeq TM2000 machine. The raw data generated from Illumina platform were submitted to Enterobase 1 , and sequence types were automatically assigned according to the Achtman scheme. A minimum spanning tree was generated using GrapeTree software to analyze the distribution of sequence type of MPECs (Zhou et al., 2018) characterized by multilocus sequence typing (MLST) analysis.

Mating Experiment of MPECs
To evaluate the transferability of mcr-1 gene, broth-mating experiments were performed with six representative MPECs (donor) and streptomycin-resistant E. coli C600 (recipient). Briefly, donor strain was mixed with recipient strain in LB broth (1:3) and the mixture was incubated at 37 • C for 4 h in a shaker incubator (180 r/min). Subsequently, mixture was spread on MAC plates containing 2,000 µg/ml streptomycin (Solarbio, Beijing, China) and 2 µg/ml colistin and incubated at 37 • C for 24 h to select transconjugants. The susceptibility of transconjugants for antimicrobials was determined based on the MICs.

Whole-Genome Sequencing Analysis of MPECs
The total genomic DNA was extracted from representative MPECs with TIANamp Bacteria DNA Kit (Tiangen, Beijing, China) following the manufacturer's instructions and sequenced (MiSeq, Illumina, San Diego, CA, United States), producing 2 × 150-bp paired-end reads (Majorbio Company, Shanghai, China) and long-read (Pacific Biosciences, Menlo Park, CA, United States) technology. The short and long reads were assembled using Unicycler (Wick et al., 2017). Gene prediction and annotation were done with RAST genome annotation server version 2.0 2 and BLASTp. Plasmid replicon types were identified by the PlasmidFinder 2.1 3 from the Center for Genomic Epidemiology (CGE). Acquired resistance genes and chromosomal mutation-mediated antimicrobial resistance genes (ARGs) were identified with ResFinder 4.1 4 from CGE and resistance gene identifier (RGI) 5 from the Comprehensive Antibiotic Resistance Database (CARD). Conjugal transfer components and type IV secretion systems (T4SS) of the plasmids were performed using oriTfinder 6 . Insertion sequence (IS) elements, transposon (Tn), and integron (In) were identified using ISfinder 7 , INTEGRALL 8 , and MobileElement Finder v1.0.3 9 . Virulence genes and FimH type were identified with VirulenceFinder 2.0 10 and FimTyper 1.0 11 , respectively, available at the Center for Genomic Epidemiology. Phylogroups were predicted using the ClermonTyping tool at the IAME Research Center web 12 .

Antimicrobial Susceptibility Profiles of 24 MPECs
After isolating bacterial strains from 109 samples, 24 MPECs were screened based on mcr-1 gene amplification and sequencing (data not shown). Among the 24 MPECs, nine strains were from poultry, 14 strains were from swine, and one strain was from bovine (Supplementary Table 2). The MICs of the 24 MPECs against 14 antimicrobials are as shown in Supplementary Table 1. All of MPECs were resistant to at least three or more antimicrobial categories; nine of 24 MPECs, namely J-1, J-2, J-3, J-4, J-5, J-6, J-7, J-8, and N-14, were resistant to all the 14 antimicrobials tested, including β-lactams (penicillins, cephems, monobactams, carbapenems), tetracyclines, fluoroquinolones, aminoglycosides, sulfonamides, phenicols, fosfomycins, and polymyxins. In addition, five of these (B-1, B-2, B-3, B-9, and18-16) were also sensitive only to meropenem (Figure 1). The above results indicated that all the 24 MPECs were resistant to at least three or all antimicrobial categories and exhibited either MDR or XDR for the common antimicrobials in Enterobacteriaceae.

Conjugative Transfer of mcr-1 Gene in Six MPECs
Based on the antimicrobial resistance spectrum and the ST distribution of the 24 MPECs, six MPECs (S-4, A-3, J-8, J-9, 19-5, and B-2) were selected for the mating assay. Transconjugants were obtained from all donors except J-8 and B-2. These transconjugants exhibited the same level of colistin resistance (4 µg/ml) as the donors, but at the same time, both showed a 16-fold increase in resistance when compared with the recipient E. coli C600 (0.25 µg/ml) ( Table 1). In the case of three E. coli transconjugants (Tc-S-4, Tc-J-9, and Tc-19-5), no additional resistance phenotype was cotransferred along with the mcr-1 gene. However, for Tc-A-3, showed β-lactam, tetracycline, fluoroquinolone, sulfonamide, fosfomycin, and florfenicol resistance phenotypes were also observed in addition to that of colistin resistance. The results suggested that mcr-1 genes of S-4, A-3, J-9, and 19-5 were located on the conjugative transfer plasmid and could be successfully transferred to a recipient strain to confer the colistin resistance phenotype. Although the mcr-1 gene of both J-8 and B-2 conjugated failure, J-8 was resistant to all tested antimicrobials. Therefore, this strain was chosen, along with the other successfully conjugated ones for the next experiment.

Features of Resistance Genes and Plasmids in Five MPECs
Five MPECs (S-4, A-3, J-8, J-9, and 19-5) were selected for whole-genome sequencing (WGS) analysis to identify the gene and plasmid features which are responsible for resistance. Beside the mcr-1 gene, all the five MPECs contained more than three different types of acquired resistance genes as well as at least two resistant plasmids ( Table 2). In addition, locus mutations of gyrA (S83L, D87N) and parC (S80I) were also found on the chromosome. In the case of J-8, the mcr-1 gene was located on the chromosome and flanked by two ISApl1 (Figure 3), with an additional bla NDM−5 gene located on the IncX3-type plasmid. Conversely, the mcr-1 genes of the other four MPECs were located on three different transferable plasmids IncHI2/HI2A/N (A-3), IncX4 (S-4), and IncI2 (J-9 and 19-5).
mcr-1 gene. Besides, on pS4-mcr-1.1, virB6 gene in T4SS was truncated into two fragments by IS2. In the case of pOW3E1, IS2 was found upstream of the mcr-1 gene. In addition, a type II toxin-antitoxin system (TAs) hicA/hicB was also found between upstream of T4SS and downstream of the chaperone-encoding dnaJ.
The mcr-1 in strain 19-5 and J-9 was, respectively, carried by the plasmid pmcr-1 (60,734 bp; GC content of 42.50%) and the plasmid pMCR_J9_7 (60,959 bp; GC content of 42.36%), all belonging to the IncI2 incompatibility group ( Table 2). BLASTn analysis showed that pMCR_J9_7 was similar to the plasmid pD90-2 (CP022452; nucleotide coverage 99%; identity FIGURE 3 | Genetic context of mcr-1 gene in chromosome and plasmids from this study compared with plasmid pHNSHP45. Black arrows represent nikA (plasmid mobilization relaxosome protein), nikB (relaxase), and hypothetical protein (hp). Green arrows denote insertion sequence ISApl1. Red arow denotes colistin resistance gene mcr-1. Regions of >99% identity are marked by gray shading. 99.99%) and pHNSHP45 (KP347127; nucleotide coverage 90%; identity 99.63%), while pmcr-1 was similar to the plasmid PN21 (MG557851; nucleotide coverage 100%; identity 99.97%) and pHNSHP45 (KP347127; nucleotide coverage 99%; identity 99.97%) (Figure 6). A sequence comparison of the two plasmids indicated that they shared 99.63% identity with 91% nucleotide coverage. The gaps shown in the circle diagram represented the differences between the five mcr-1-bearing plasmids and the absence of IS683 and ISApl1. For colinear analysis with the plasmids pMCR_J9_7 and pmcr-1, pHNSHP45 was selected as reference one. As shown in Figure 7, ISApl1 located upstream of the mcr-1 gene and IS683 located between repA and parA were absent in plasmid pMCR_J9_7 and pmcr-1. In addition, stbD and stbE genes, both of which are located downstream of virB5 and encoding stability proteins, were absent from the plasmid pMCR_J9_7. A similar observation was made for dnaJ, a chaperone-encoding gene, which was absent between parA and traL in plasmid pmcr-1. For both identified IncI2-type plasmids, only the mcr-1 gene, embedded within the mcr-1-pap2 cassette, was associated with colistin resistance (Figure 3). Moreover, a type II toxin-antitoxin system (TAs) hicA/hicB was found to be present on the plasmid pMCR_J9_7, with the hicA gene, encoding antitoxin proteins, being absent from the plasmid pmcr-1.
FIGURE 7 | Colinear sequence comparison of IncI2 reference plasmid pHNSHP45 (GenBank accession No. KP347127) with plasmids including pmcr-1 (this study) and pMCR_J9_7 (this study). Boxed arrows represent the position and transcriptional direction of ORFs. Dark gray shading denotes regions of homology (>99% identity nucleotide identity). Blue colors represent the plasmid replication-and maintenance/stability-associated genes; teal, plasmid conjugative transfer genes and T4SS genes; green, IS and Tn; red, antimicrobial resistance genes; pink, TAs genes; gray, other genes.
FIGURE 9 | Colinear sequence comparison of IncHI2/HI2A/N plasmid pMCR_A3_2 (this study) with other plasmids including pMCR1_025943 (GenBank accession No. CP027202) and pHNSHP45-2 (GenBank accession No. KU341381). Boxed arrows represent the position and transcriptional direction of ORFs. Dark gray shading denotes regions of homology (>99% identity nucleotide identity). Blue colors represent the plasmid replication-and maintenance/stability-associated genes; teal, plasmid conjugative transfer (tra-/trb-gene cluster) and T4SS genes; green, IS and Tn; red, antimicrobial resistance genes; pink, TA genes; gray, other gene. The light green shade represents the composition transposon and integron (In); MDR, multidrug resistant.

DISCUSSION
MPECs from different sources have been isolated worldwide. In the present study, 24 MPECs were isolated from organs having lesions obtained from animals, which did not respond to antimicrobial treatment and which were, therefore, considered to show MDR, or even XDR based on the definition provided by Magiorakos et al. (2012). MLST analysis classified the 24 MPECs into nine sequence types (six dominant CCs and two singletons), among which ST10596 has been identified for the first time. The main CCs were CC156 (ST156), CC448 (ST448), and CC10 (ST48 and ST10).
Escherichia coli ST10 and ST48 (single locus variant of ST10), belonging to CC10, are considered widely disseminated clones worldwide. CC10 is highly disseminated via foods, the environment, animals, and humans (Matamoros, 2017), and they are usually associated with the presence of colistin (mcr-1) (Papa-Ezdra et al., 2020), ESBL (bla CTX−M−14/−15 ), and fluroquinolone [aac(6 )-Ib-cr and qnrB] (Oteo et al., 2009;Fam et al., 2011) resistance genes. However, in this study, in addition to those genes, MPEC ST10 (S-4) and 19-5 (ST48) were also found to contain aminoglycoside, tetracycline, sulfonamide, macrolide, lincosamide, and phenicol, as well as peroxides resistance genes. These results suggested that CC10 possessed a more extensive set of resistance genes and it could, therefore, become a reservoir for the transmission of resistance, and CC156 also could be a reservoir for mcr-1 genes in a similar way to CC10. Furthermore, MPECs of ST448 are generally sporadic and not dominant (Bai et al., 2018), but results showed that six MPECs of ST448 were members of CC448, hence suggesting that CC448 could become a new dominant clone complex carrying the mcr-1 gene in pig farms.
The mcr-1 gene along with other resistance genes are found on conjugative transfer plasmids or mobile element transposons. As such, there is a risk for the horizontal transmission of colistin as well as other drug resistance between bacteria in the infected guts of animals, thereby further exacerbating the risk of spreading resistance (Guenther et al., 2017;Bai et al., 2018;. The results of this study showed that the mcr-1 genes of strains S-4, A-3, J-9, and 19-5 located on the three different transferable plasmids IncHI2/HI2A/N (A-3), IncX4 (S-4), and IncI2 (J-9, 19-5), could be successfully transferred. Despite the fact that IncI2-, IncX4-, and IncHI2-type plasmids were the most common carriers of the mcr-1 gene, yet their mobile elements (IS, Tn, In) were very different. Therefore, the study and description of mobile elements on plasmids can greatly contribute to the understanding of the coupling transfer of plasmids, the formation of hybrid plasmids, and the translocation of ARGs.
Self-transmissible IncX4-type plasmid, a key vector for the transmission of the mcr-1 gene between Enterobacteriaceae worldwide, is the most common one which carries the resistance gene due to its high conjugative transfer ability (Biswas et al., 2020;Migura-Garcia et al., 2020). The structure of mcr-1carrying IncX4-type plasmids are highly conserved, and they consist of mcr-1-pap2 gene cassette, T4SS systems, plasmid replication/maintenance genes, TA modules, and an insertion sequence IS26 (Bai et al., 2018). The IncX4-type plasmid pOW3E1, in addition to the previously mentioned modules, contained IS2 close to the mcr-1 gene (Zurfluh et al., 2016), but the insertion sequence had no effect on either the gene or the conjugative transfer of plasmid. In the present study, IS26 and IS2 were found on the IncX4-type plasmid pS4-mcr-1.1, with IS2 truncating the virB6 gene of T4SS. Since the virB6 product contributes to DNA transfer as a component of the endosomal translocation channel protein (Álvarez-Rodríguez et al., 2020), its truncation may affect the efficiency of plasmid conjugative transfer, although the process through which this occurs remains unknown.
The first reported mcr-1-carrying IncI2-type plasmid pHNSHP45 was identified in E. coli of swine origin in China and ISApl1 was located upstream of the resistance gene . However, in this study, the two IncI2-type plasmids (pmcr-1 and pMCR_J9_7) did not have any ISApl1 adjunct to the mcr-1 gene, with IS683 also lacking between repA and parA. Furthermore, colinear analysis of the plasmid pHNSHP45 along with those of pMCR_J9_7 and pmcr-1 indicated the absence of stbD and stbE genes, both encoding plasmid stability proteins, from pMCR_J9_7. Similarly, the dnaJ gene, encoding for the chaperone protein, was absent from the pmcr-1 plasmid. These results implied that the IncI2-type plasmid could lose some insertion sequences during the evolution process to maintain the conservative plasmid backbone structure, which is conducive to the stability of the mcr-1 gene and the conjugative transfer of the plasmid.
IncHI2/HIA-type plasmids are usually > 250 kb, unlike IncX4-and IncI2-type plasmids. In addition, they contain an MDR region with various ISs, transposons, and integrons that capture various exogenous resistance genes, thereby conferring multidrug resistance and metal tolerance to E. coli . However, despite these differences, the backbone of IncHI2/HI2A-type plasmids carrying mcr-1 gene are also conserved, in a similar way to those of the IncX4 and IncI2 plasmids. Hence, they were found to include plasmid replication/maintenance genes, T4SS systems, tellurium ion resistance gene clusters, and TA modules. Studies on the IncHI2type plasmids carrying the mcr-1 gene mostly focused on the description of the genetic environment of the gene, while the mobile elements of the MDR region of the IncHI2 type were overlooked. In this context, it is to be noted that the IncHI2/HI2A-type plasmid pEGY1-MCR-1 carrying the mcr-1 gene, is flanked by two ISApl1 elements as well as In641 (Hammad et al., 2019). In this study, pMCR_A3_2, a IncHI2/HI2A/Ntype hybrid plasmid with an extra IncN-type replicon, was flanked by two IS26 present in the MDR region. This implied that the replicon could have been transposed by a transposon composed of IS26. Alternatively, this could have been the result of homologous recombination or the transposase-mediated formation of a cointegrate from IS26 and the nic site of oriT. Additional results showed that the MDR region of the pMCR_A3_2 plasmid, identified in this study, included two class 1 integrons, namely, In0 and In640. However, that of the IncHI2-type plasmid pHNSHP45-2 carrying the mcr-1 gene has been reported to contain two different class 1 integrons, namely, In1185 and In641 (Zhi et al., 2016). Furthermore, In1519 was also identified from plasmid pMCR1_025943 (CP027202). Apart from those class 1 integrons, a large number of insertion sequences and transposons such as IS26/26-like (15DI) and Tn3 families, were found to be present in the MDR region. The presence of these mobile integration elements facilitates the capture of numerous classes of resistance genes, thereby conferring various resistance phenotypes to E. coli containing the IncHI2-type plasmids.
Generally, ISApl1 is detected upstream of the mcr-1 gene and is believed to be involved in the mobilization of the mcr-1 cassette (Snesrud et al., 2016). In this study, the ISApl1 was absent from both upstream and downstream of the mcr-1 gene on IncX4, IncI2, and IncHI2/HI2A/N plasmids, even though mcr-1 on the chromosome was embedded in Tn6330. Previous studies have shown that the transposon Tn6330 was immobilized on several plasmid backgrounds, following the loss of the flanking ISApl1 elements due to the high instability of being spread through plasmid transfer Wang et al., 2018). The loss of the ISApl1 transposon element in mcr-1-carrying plasmids has recently been proposed to be essential for the maintenance of the mcr-1 gene. This is believed to be more beneficial to the host bacteria for adapting to environmental changes, especially when adapting from the pressure of antimicrobial agents to a pressure-free environment (Porse et al., 2016).

CONCLUSION
In conclusion, 24 MPECs in our study showed MDR or even XDR and exhibited the population diversity. The T4SS genes located on IncX4-and IncHI2/HI2A/N-type plasmids were truncated by the insertion sequence. The truncation may affect the efficiency of plasmid conjugative transfer. Additionally, the MDR region of IncHI2/HI2A/N-type plasmids contained two class 1 integrons and four composite transposons. These mobile elements contributed to the multireplicon plasmid formation and the acquisition and transfer of ARGs.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.

ETHICS STATEMENT
The animal study was reviewed and approved by all animal procedures and study design were approved by the animal Ethics Committee of Northwest A&F University (Approval No: 2016012). Written informed consent was obtained from the owners for the participation of their animals in this study.

AUTHOR CONTRIBUTIONS
XS and ZL: conceptualization. ZL: methodology, software, formal analysis, data curation, and visualization. XS, YL, YF, and WM: validation. JL, HM, SL, BC, and HH: investigation. WZ: resources. ZL and YL: writing -original draft preparation. XS and YL: writing, review, and editing. XS: supervision, project administration, and funding acquisition. All authors have read and agreed to the published version of the manuscript.