The animals used for genetic variation identification were provided by the Lanxi breeding farm (Lanxi, Heilongjiang, China). A total of 226 Min and 186 Landrace pigs' genomic DNA were used for the SNPs genotyping and association analysis and the phenotype traits of all these animals have been described previously (14). In brief, the piglets' feces were observed and assigned a daily diarrhea score according to the following standard: normal solid feces 0, slight diarrhea with soft and loose feces 1, moderate diarrhea with semi-liquid feces 2, and severe diarrhea with liquid and unformed feces 3 (37). The piglet diarrhea score was calculated as the summation of the daily diarrhea score during the experiment. Additionally, the body weight at birth, 21-day, 28-day, and 35 days were recorded and the average daily gain (ADG) was calculated as follows, ADG = (weight at 35-day - weight at birth)/35.

Animal experimentation was conducted in line with the guidelines for the care and use of animals approved by the Laboratory Animal Management Committee of Northeast Agricultural University (Harbin, Heilongjiang, China).

The core promoters of the porcine CISH gene were analyzed using the online neural network promoter prediction (BDGP) software (http://www.fruitfly.org/seqtools/promoter.html). JASPAR software (http://jaspar.binf.ku.dk/cgi-bin/jaspar_db.pl) and PROMO (http://www.urogene.org/cgi-bin/methprimer/methprimer.cgi) were used to search the putative transcription factor binding sites. Functional SNPs in the coding sequence (CDS) of porcine STAT5A were predicted using the in silico method. The information on nsSNPs in the porcine STAT5A gene was first collected from the Ensembl SNP database (http://www.ensembl.org/). Five bioinformatics tools including PhD-SNP (https://snps.biofold.org/phd-snp/), SIFT (http://sift.jcvi.org/), SNAP (http://www.rostlab.org/services/SNAP), Meta-SNP (http://snps.biofold.org/meta-snp/), and PolyPhen-2.0 (http://genetics.bwh.harvard.edu/pph2/uses), which were described in more detail by Emadi et al. (38) and Li et al. (33), were then employed to predict the most deleterious nsSNPs based on incorporating the scores of all the servers. SNPs in the 3'UTR of porcine STAT5A were analyzed by miRBase (http://www.mirbase.org) to screen the putative functional SNPs which might alter the microRNA (miRNA) response elements. Similarly, SNPs in the 5' UTR and intron 1 were analyzed by JASPAR and PROMO with a threshold score of 80.

To obtain the porcine CISH promoter, the 2.5 kb upstream genomic sequence of porcine CISH was retrieved from Ensembl, and primer pair CISH-C was designed using Primer 5 software (Supplementary Table S1). The PCR reaction contained 100 ng DNA templates, 0.5 μM of each primer as in Supplementary Table S1, and 10 μL 2 × Taq Master Mix (TaKaRa, Dalian, China). The PCR conditions were 4 min at 94 °C; 35 cycles of 30 s at 94 °C, 30 s at the annealing temperature (Supplementary Table S1), 2 min at 72 °C and a final step at 72 °C for 8 min. PCR products were purified, sequenced, and ligated into PMD18-T (TaKaRa, Dalian, China) to produce PMD18-T-CISH. Then, using PMD18-T-CISH as a template, a series of 5' deletion fragments of CISH promoter were amplified using forward primers CISH-C (C1-C5), CISH-P (P1-P3), and the reverse primer CISH-R (Supplementary Table S1). To explore the effect of GATA1 on CISH transcription, primer CISH-G was designed and a specific fragment containing putative GATA1 binding element was produced by PCR (Supplementary Table S1). All these PCR products were purified, digested with Xho I and Kpn I or Hind III (TaKaRa, Dalian, China), and inserted into the pGL3-Basic vector (Promega, Madison, Wisconsin, USA) to produce luciferase reporter plasmids pGL3-CISH-C (C1-C5), pGL3-CISH-P (P1-P3) and pGL3-CISH-G.

Using pGL3-CISH-P1 plasmid as a template and the mutagenic primers CISH-Mut1 and CISH-Mut2 (Supplementary Table S1), three mutants of predicted STATx binding sites were generated by PCR. Then these PCR products were purified, digested by restriction endonucleases, and inserted between Xho I and Kpn I sites of the pGL3-basic vector. All the reconstructed plasmids were confirmed by DNA sequencing and named pGL3-CISH-mut1, pGL3-CISH-mut2, and pGL3-CISH-mut3.

To construct pGL3-STAT5A-Luc and pGL3-STAT5A-SNP1, specific fragments containing different variations in intron 1 of porcine STAT5A were amplified with the primer pair STAT5A-Luc (Supplementary Table S1) and the genotyped genomic DNA, inserted into the pGL3-basic vector. Then, using pGL3-STAT5A- Luc plasmid as the template and the mutagenic primers STAT5A-Mut (Supplementary Table S1), the third fragment was produced and the corresponding pGL3-STAT5A-SNP2 was constructed as described above. All the positive plasmids were confirmed by endonuclease digestion and DNA sequencing.

Hela and IPEC-J2 were selected as in vitro model systems to explore porcine CISH core promoters. The cells were cultured in 24-well plates with DMEM supplemented with 10% FBS (Gibco, Carlsbad, CA, USA), maintained at 37 °C in the atmosphere of 5% CO₂ for 18–24 h. Then, the cells were transfected with 0.5 μg of specific promoter-luciferase plasmid, 0.005 μg of pRL-TK used as an internal control (Promega, Madison, Wisconsin, USA), and 1.5 μL of X-treme GENE HP DNA Transfection Reagent (Roche, USA). Simultaneously, the plasmid pGL3-Basic mixed with pRL-TK was transfected into the corresponding cells as the negative control of the luciferase report assay. After 24–48 h, all the cells were harvested and lysed with the luciferase assay buffer (Promega, Madison, Wisconsin, USA), and the enzymatic activity of firefly and renilla was examined according to Dual-Luciferase Reporter Assay System (Promega, Madison, Wisconsin, USA). Transfection was performed in triplicate and repeated two or three times, and the luciferase activity was calculated as the ratio of firefly to renilla.

The complete coding region of porcine STAT3 (NM_001044580.1), STAT5A (NM_214290.1), and GATA1 (NM_001278767.1) were amplified using porcine cDNA and primer pairs (Supplementary Table S1) respectively. The purified PCR products were inserted between the Xho I and EcoR I or Not I sites of the pCMV-HA vector (Promega, Madison, WI, USA) to construct pCMV-HA-STAT3, pCMV-HA-STAT5A, and pCMV-HA-GATA1, respectively. These expression plasmids were confirmed by DNA sequencing and Western blot analysis. Then, IPEC-J2 cells were transfected with 0.25 μg of CISH promoter luciferase reporter plasmid (pGL3-CISH-C4, pGL3-CISH-P1, or pGL3-CISH-G), 0.25 μg of corresponding expression plasmid (pCMV-HA-STAT3, pCMV-HA-STAT5A, pCMV-HA-GATA1, or pCMV-HA), and 0.005 μg of pRL-TK. Relative luciferase activity was calculated as described above.

To validate the role of transcriptional factors on porcine CISH mRNA expression, IPEC-J2 cells were transfected with 2 μg of pCMV-HA-STAT3, pCMV-HA-STAT5A, pCMV-HA-GATA1, or pCMV-HA, respectively. All the cells were collected, and RT-qPCR was performed.

For the Western blot, IPEC-J2 cells transfected with pCMV-HA-STAT3, pCMV-HA-STAT5A, pCMV-HA-GATA1, or pCMV-HA were lysed with RIPA buffer (SEVEN, Beijing, China), boiled with 5 × denaturing loading buffer, resolved in 12% SDS-PAGE, and transferred to Immuno-Blot polyvinylidene fluoride (PVDF) membrane (Millipore, Billerica, MA, USA). Then, the membrane was blocked with 5% BSA, washed with TBST, and incubated with HA-tag antibody (ABclonal, Wuhan, China) in 1:1000 dilution for 12 h at 4 °C, and a secondary antibody at room temperature for 2 h. Lastly, this incubated membrane was washed with TBST and Super ECL kit (SEVEN, Beijing, China) was used to show the bands.

For the RNA extraction and RT-qPCR, total RNAs were extracted from IPEC-J2 cells using TRIzol (Takara, Dalian, China) and reverse transcribed into cDNA using PrimerScript RT Master Mix (Takara, Dalian, China). According to SYBR Premix Ex Taq (Takara, Dalian, China), RT-qPCR was conducted on an ABI 7500 system (Applied Biosystems, Foster City, CA, USA) with the following reaction: 100 ng cDNA, 0.2 μM of each primer (Supplementary Table S1), and 10 μL SYBR mix (Takara, Dalian, China). The qPCR conditions were 95 °C for 30 s, 40 cycles with 95 °C for 5 s, and 60 °C for 35 s, with the melting curves constructed at the same time. In line with the 2^−ΔΔCT method (35), Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) (36, 37) was selected as a housekeeping gene to normalize the expression of target genes and relative mRNA expression of a specific gene in cells was calculated.

The 5'UTR and 3'UTR of porcine STAT5A were amplified using the primer pair STAT5A-5'UTR or STAT5A-3'UTR (Supplementary Table S1) and the mixed genomic DNA pools which consisted of five healthy Min, five diarrheal Min, five healthy Landrace, and five diarrheal Landrace pigs described above. All the PCR products were purified, sequenced commercially (Sangon, Shanghai, China), and aligned with Clustal Omega. The predicted functional SNPs by in silico tools were genotyped in Min and Landrace populations through PCR-based single-strand conformation polymorphism (PCR-SSCP). In brief, fragments containing SNPs were amplified by PCR using primer pair STAT5A-SNP (Supplementary Table S1). Then 1 μL PCR products mixed with 9 μL denaturation buffer were denatured for 10 min at 98 °C, placed on iced water for 5 min, separated on a 14% PAGE gel, and resolved by silver staining. PCR products with diverse genotypes were purified and sequenced commercially (Sangon, Shanghai, China).

The significant difference was evaluated using the t-test in GraphPad Prism 5 software (GraphPad, La Jolla, CA, USA) or one-way ANOVA. The observed heterozygosity (Ho), expected heterozygosity (He), effective allele numbers (Ne), and chi-square test for Hardy-Weinberg equilibrium of the genetic variation polymorphisms were calculated using Popgene software (version 1.32). The Haploview 4.2 software was used to analyze the linkage disequilibrium (LD) between g.508A>C and g.566C>T and the haplotype was constructed.

A mixed linear model procedure of SAS version 8.0 was used to perform association analysis between SNP or haplotype and the phenotype trait in Min pigs. The following statistical model was adopted:

Where Y_ij is the observed phenotype traits, μ is the population mean, G_i is the fixed effect of genetics, S_j is the random effect of the sow, and e_ij is the random residual.

A 2071-bp fragment containing a 5'-flanking region and 106-bp coding sequence was produced by PCR, and four putative core promoters in the region of −1,807/−1,758, −1,358/−1,309, −460/−411, and −201/−152 were predicted by BDGP online software (Figure 1A). Subsequently, five 5'-deletions were amplified and inserted into the pGL3 plasmid to produce the CISH promoter luciferase reporter gene plasmids. As shown in Figure 1B, all these plasmids had luciferase activity in both IPEC-J2 and Hela cells. Further analysis showed the CISH-C4 fragment (sequence from −412 to 106 bp) had higher activity than the CISH-C5 fragment (sequence from −258 to 106 bp), suggesting the existence of important positive regulatory transcription elements in −412/−258 bp. However, the activity of porcine CISH-C2 (sequence from −1,498 to 106 bp) was lower than CISH-C3 (from −528 to 106 bp), indicating there were elements in the region of −1,498/−528 that could negatively regulate the transcription of porcine CISH.

In silico analysis showed the CISH-C4 fragment contained one core promoter and tetramer STATx transcription factor binding sites (Figure 2A). To verify this prediction, three 5'-deletion plasmids were constructed. Luciferase activity analysis showed the deletion of the first two putative STATx binding sites (−412/−338) did not alter the luciferase activity of fragment CISH-C4 (Figure 2B). However, the deletion of the third and fourth STATx binding sites (−338/−227) significantly reduced the luciferase activity of CISH-C4 (Figure 2B). Additionally, as shown in Figure 2B, the disappeared luciferase activity of the CISH-P3 fragment indicated the existence of a minimal core promoter in the region of −227/−106. To verify the existence of STATx binding sites, the predicted third or fourth STATx binding sites within the CISH-P1 fragment were mutated (Figure 2C). Luciferase activity analysis revealed that the mutations reduced the luciferase activity of CISH-P1 significantly (Figure 2D).

In this study, based on the literature, porcine STAT3 and STAT5A were selected to verify their effect on CISH transcriptional activity and mRNA expression. The expression plasmid pCMV-HA-STAT3 or pCMV-HA-STAT5A, together with CISH promoter luciferase reporter plasmids pGL3-CISH-C4 or pGL3-CISH-P1, were transfected into IPEC-J2. The luciferase activity analysis showed that the overexpression of STAT3 could significantly increase CISH transcriptional activity (P < 0.01) (Figures 3A–C). Similarly, the overexpression of STAT5A up-regulated CISH transcriptional activity (P < 0.01) (Figures 3D–F). The RT-qPCR revealed that the overexpression of STAT3 did not affect CISH mRNA expression (Figures 3G, H). However, STAT5A promoted the CISH expression at 48 h after transfection (P < 0.01) (Figures 3I, J).

In silico analysis showed the CISH-G fragment containing GATA1 binding sites upstream of the tetramer STATx transcription factor binding sites (Figure 4A). The reconstruction plasmid pCMV-HA-GATA1 and pGL3-CISH-G were transfected into IPEC-J2. The luciferase activity analysis showed that GATA1 overexpression significantly repressed CISH transcriptional activity (P < 0.01) (Figures 4B, C). The RT-qPCR revealed that the overexpression of GATA1 reduced the CISH mRNA expression at 24 and 48 h after transfection (Figures 4D, E).

In this study, a total of 49 SNPs in the porcine STAT5A coding region were retrieved from the Ensembl database. Among the 14 nsSNPs, only c.870G>T was the predicted deleterious nsSNP with Polyphen-2 (Supplementary Table S2). Then, a 1618-bp fragment containing a 210-bp intron17, 1384-bp exon18 (containing 1221-bp 3'UTR), and 24-bp downstream region of porcine STAT5A was amplified. Although 9 SNPs were identified in this fragment, none of these SNPs was predicted to alter the miRNA binding sites (data not given). Similarly, a 909-bp fragment containing 231-bp 5'UTR and 678-bp intron1 of porcine STAT5A was amplified and three SNPs (g.373C>G, g.508A>C, and g.566C>T) were identified in intron 1. According to JASPAR, g.373C>G caused a zinc finger protein 454 (ZNF454) polymorphism; g.508A>C was predicted to alter the binding site of E2F transcription factor 4 (E2F4) or transcription factor Dp-1 (TFDP1), where E2F4 or TFDP1would bind to the C allele but not the A allele; g.566C>T altered the binding sites of PR domain containing 4 (PRDM4) or RELB proto-oncogene, and NF-kB subunit (RELB) (Supplementary Table S3). However, among these SNPs, PROMO revealed only g.566C>T altered the binding of ETS transcription factor ELK1 (Elk-1) (Supplementary Table S2). Collectively, the in silico analysis indicated that g.508A>C or g.566C>T were the putative functional SNPs.

g.508A>C and g.566C>T were genotyped by PCR-SSCP in Min and Landrace populations. As shown in Figure 4, three genotypes (AA, AC, and CC) of g.508A>C (Figures 5A) were observed, while only CC and CT genotypes of g.566C>T (Figures 5A, D, E) were found. The frequency of the A allele (g.508A>C) was 0.61 in Min pigs, and the C allele (g.566C>T) was 0.94 in this population. However, these two alleles were fixed in the Landrace population. In Min pigs, the Ho and He of g.508A>C were 0.49 and 0.48, while Ho and He of g.566C>T were 0.11 and 0.11, respectively; the effective allele number (Ne) of g.508A>C and g.566C>T was 1.91 and 1.12, respectively; and the polymorphism information content (PIC) values were 0.36 and 0.10, respectively (Table 1). Statistical analysis showed g.508A>C was not associated with Min piglets' phenotype traits, while animals with CC genotype at g.566C>T had more 28-day body weight, 35-day body weight, and ADG than those of CT genotype (P < 0.05) (Table 2). Additionally, three haplotypes were constructed and Hap1:AC was the major haplotype with a frequency of 0.61, while the minor Hap3:CT was 0.06 (Table 3). Association analysis indicated Min piglets with Hap1:AC had higher 28-day body weight, 35-day body weight, and ADG when compared to Hap3:CT (P < 0.05) (Table 3).

As shown in Figure 6A, the STAT5A promoter luciferase reporter pGL3-STAT5A-Luc (g.508A and g.566C), pGL3-STAT5A-SNP1 (g.508C and g. 566C), and pGL3-STAT5A-SNP2 (g.508A and g.566T) were constructed to verify the putative function of g.508A>C or g.566C>T. After transfection, the luciferase activity of these reporter vectors was calculated and compared. The statistical analysis indicated the luciferase activity of g.508A>C A allele was higher than that of the C allele (P < 0.05), while the difference between the C and T allele of g.566C>T was not significant (Figure 6B).