Three Novel Players: PTK2B, SYK, and TNFRSF21 Were Identified to Be Involved in the Regulation of Bovine Mastitis Susceptibility via GWAS and Post-transcriptional Analysis

Bovine mastitis is a common inflammatory disease caused by multiple factors in early lactation or dry period. Genome wide association studies (GWAS) can provide a convenient and effective strategy for understanding the biological basis of mastitis and better prevention. 2b-RADseq is a high-throughput sequencing technique that offers a powerful method for genome-wide genetic marker development and genotyping. In this study, single nucleotide polymorphisms (SNPs) of the immune-regulated gene correlative with mastitis were screened and identified by two stage association analysis via GWAS-2b-RADseq in Chinese Holstein cows. We have screened 10,058 high quality SNPs from 7,957,920 tags and calculated their allele frequencies. Twenty-seven significant SNPs were co-labeled in two GWAS analysis models [Bayesian (P < 0.001) and Logistic regression (P < 0.01)], and only three SNPs (rs75762330, C > T, PIC = 0.2999; rs88640083, A > G, PIC = 0.1676; rs20438858, G > A, PIC = 0.3366) were annotated to immune-regulated genes (PTK2B, SYK, and TNFRSF21). Identified three SNPs are located in non-coding regions with low or moderate genetic polymorphisms. However, independent sample population validation (Case-control study) data showed that three important SNPs (rs75762330, P < 0.025, OR > 1; rs88640083, P < 0.005, OR > 1; rs20438858, P < 0.001, OR < 1) were significantly associated with clinical mastitis trait. Importantly, PTK2B and SYK expression was down-regulated in both peripheral blood leukocytes (PBLs) of clinical mastitis cows and in vitro LPS (E. coli)–stimulated bovine mammary epithelial cells, while TNFRSF21 was up-regulated. Under the same conditions, expression of Toll-like receptor 4 (TLR4), AKT1, and pro-inflammatory factors (IL-1β and IL-8) were also up-regulated. Interestingly, network analysis indicated that PTK2B and SYK are co-expressed in innate immune signaling pathway of Chinese Holstein. Taken together, these results provided strong evidence for the study of SNPs in bovine mastitis, and revealed the role of SYK, PTK2B, and TNFRSF21 in bovine mastitis susceptibility/tolerance.

Bovine mastitis is a common inflammatory disease caused by multiple factors in early lactation or dry period. Genome wide association studies (GWAS) can provide a convenient and effective strategy for understanding the biological basis of mastitis and better prevention. 2b-RADseq is a high-throughput sequencing technique that offers a powerful method for genome-wide genetic marker development and genotyping. In this study, single nucleotide polymorphisms (SNPs) of the immune-regulated gene correlative with mastitis were screened and identified by two stage association analysis via GWAS-2b-RADseq in Chinese Holstein cows. We have screened 10,058 high quality SNPs from 7,957,920 tags and calculated their allele frequencies. Twenty-seven significant SNPs were co-labeled in two GWAS analysis models [Bayesian (P < 0.001) and Logistic regression (P < 0.01)], and only three SNPs (rs75762330, C > T, PIC = 0.2999; rs88640083, A > G, PIC = 0.1676; rs20438858, G > A, PIC = 0.3366) were annotated to immune-regulated genes (PTK2B, SYK, and TNFRSF21). Identified three SNPs are located in non-coding regions with low or moderate genetic polymorphisms. However, independent sample population validation (Case-control study) data showed that three important SNPs (rs75762330, P < 0.025, OR > 1; rs88640083, P < 0.005, OR > 1; rs20438858, P < 0.001, OR < 1) were significantly associated with clinical mastitis trait. Importantly, PTK2B and SYK expression was down-regulated in both peripheral blood leukocytes (PBLs) of clinical mastitis cows and in vitro LPS (E. coli)-stimulated bovine mammary epithelial cells, while TNFRSF21 was up-regulated. Under the same conditions, expression of Toll-like receptor 4 (TLR4), AKT1, and pro-inflammatory factors (IL-1β and IL-8) were also up-regulated. Interestingly, network analysis indicated that INTRODUCTION Bovine mastitis is the most complex and costly inflammatory disease with high incidence, which seriously affects developing dairy industry worldwide (1)(2)(3). Although the breeding conditions of dairy cows have improved, mastitis is still a concern. Dairy cow mastitis causes direct economic losses in many ways, including decreased milk production, increased treatment cost, antibiotic residues, bacterial contamination, and even increased elimination rate due to cow death (4)(5)(6)(7). Previous studies have shown that cow mastitis has complex traits (clinical and subclinical mastitis) and is affected by multiple factors, including genetic features, nutrition, and hygiene standard, pathogen infections, and seasonal changes. Among them, pathogen infection (Gram-negative, such as Escherichia coli; Gram-positive, such as Staphylococcus aureus and Streptococcus uberis) is the main cause of mastitis (8)(9)(10)(11)(12)(13). Therefore, rapid identification of pathogens is critical to determine the symptoms of infection (14). However, mastitis-related immune response is a complex biological process involving immune cells, mammary epithelial cells, and endothelial cells (8). It is well-known that mastitis causes a dramatic increase in bovine milk somatic cell counts (SCC), mainly neutrophils (8,15,16). It is also known that somatic cell score (SCS) is the primary trait for detection of mastitis with high hereditary capacity (17). Thus, screening and identifying susceptible or resistant genes associated with SCC or SCS will provide novel strategies to reduce the incidence of mastitis and improve the quality of dairy cows populations (17)(18)(19).
It has been recognized that better understanding of the genetic and biological basis of a complex environment will facilitate genome prediction and the development of appropriate control strategies (20)(21)(22)(23)(24). In the past decade, different research strategies, including single nucleotide polymorphism (SNP) in a candidate gene, quantitative trait loci (QTL) and GWAS (12,25,26), were successfully used to identify the genes significantly associated with the mastitis traits. GWAS has achieved unprecedented success in identifying gene regions and candidate gene variations that are closely related to clinical phenotypes and disease susceptibility (chromosome and gene level, in term of the association between SNPs and traits) (27)(28)(29). In addition, identification of associated gene mutations can help reveal the pathogenesis of disease, provide an entry point for treatment, and analyze common genetic variations to identify multiple risk loci for complex diseases (30,31). GWAS is widely regarded as a potent method to identify SNPs in dairy cattle mastitis traits (17,32). GWAS data showed that Bos Taurus autosome 2, 4, 6, 10, 14, 18, and 20 associated with clinical mastitis were significantly correlated with SCS in cows (33)(34)(35). Additionally, two clinical mastitis candidate genes [vitamin D-binding protein precursor (GC) and neuropeptide FF receptor 2 (NPFFR2)] were detected using high-density single nucleotide polymorphic array and genomic sequencing (18,26). Wu et al. detected five mastitis susceptibility genes (NPFFR2, SLC4A4, DCK, LIFR, and EDN3) in Danish Holsteins (36). In 2015, Wang et al. identified another two mastitis susceptibility genes (TRAPPC9 and ARHGAP39) in Chinese Holstein (17). In the same year, Wu et al. confirmed the differential expression of TLR4/NF-κB signaling pathway-related genes (up-regulation genes: TLR4, MyD88, IL-6, and IL-10; down-regulation genes: CD14, TNFα, MD-2, IL-1β, NF-κB, and IL-12) in the mammary gland of Chinese Holstein mastitis cows by Gene-Chip microarray (37). Genetic variations in immune response, specific pathogen (LY75, DPP4, ITGB6, and NR4A2) and lymphocyte antigen-6 complex genes (LY6K, LY6D, LYNX1, LYPD2, SLURP1, and PSCA) might lead to clinical mastitis in American Holstein cows (38). Additionally, single gene polymorphisms (CXCR1, MAP4K4) and their signaling pathways (TLR4/NF-κB) served as genetic markers for mastitis in different cow populations (12,39). However, there is some inconsistency in genetic variation or polymorphism in genes associated with mastitis traits. Therefore, genetic mutations and polymorphisms in mastitis-associated genes should be screened and validated in different populations.
Type IIB restriction site-associated DNA sequencing (2b-RADseq) is a high-throughput genotyping-by-sequencing (GBS) method based on type IIB restriction endonuclease (for instance, BsaXI and AlfI), which provides a powerful method for identifying gene SNP in the population genome (40,41). It has strong technical repeatability, uniform depth of sequencing, high cost-effectiveness and genome wide coverage (40,42). Furthermore, 2b-RADseq technique successfully predicted multi-locus sequence typing (MLST) and provided more details on the population information than MLST technique (41). In addition, this method is also suitable for creating high-density genetic or linkage maps of genomic region or locus markers and revealing the regions associated with related traits by QTL mapping and association analysis (43,44). More importantly, 2b-RAD can identify key SNPs associated with traits by deep sequencing in fewer samples, thereby preliminary labeling candidate genes (45). Therefore, 2b-RAD may be an ideal genomic screening and labeling platform for detecting mastitis resistance or susceptibility genes in dairy cattle.
In this study, GWAS-2b-RADseq was used to screen and identify mastitis susceptibility or resistance SNPs in Chinese Holstein, and to validate the reliability of key significant SNPs in another independent sample population. Then, immuneregulated or inflammatory associated genes were marked based on the results of significant differential SNPs Gene Ontology (GO) annotation. Finally, the expression levels of candidate genes were evaluated both in vitro (peripheral blood leukocytes, PBLs) and in vivo (bovine mammary epithelial cells, bMECs). Understanding the genetic and biological pathways of bovine mastitis susceptibility or resistance candidate genes is a worthwhile strategy for anti-mastitis research and application.

Sample Libraries and Preparation
The procedures involving animals were approved by the Animal Welfare Committee of Nanjing Agricultural University, Nanjing, China, and approval No.20160615. And all animal experiments were carried out in strict accordance with the guidelines and rules established by the committee.
Chinese Holstein cows used for 2b-RADseq library and independent sample population validation came from two different pastures of the Nanjing Weigang Dairy Co., Ltd. At 2b-RADseq stage, forty dairy cows, daughters of ten bulls, were selected from 596 lactating Chinese Holstein in the third earlylactation period (15-60 d), which were further divided into two subgroups according to their clinical mastitis phenotypes [mammary gland (visual observation): red, swollen, fever, pain, and milk flocculation, etc., and laboratory diagnosis (mainly SCC)] (46,47): case group [20 cows, clinical mastitis appeared in all three lactation periods, and SCC values [(2100.3 ± 891.1) × 10 4 ) > 100 × 10 4 cells/ml] and control group [20 cows, no mastitis during the three lactation periods, and SCC values [(22.4 ± 6.8) × 10 4 ) < 30 × 10 4 cells/ml]. At the independent sample population validation stage, 383 cows (73 in case group and 310 in control group), daughters of thirty-one bulls, were selected from 886 lactating cows in the third early-lactation period (15-60 d) in another pasture for population genetics verification. All experimental cows' mastitis clinical traits were diagnosed by pasture veterinarians, and SCC were obtained from pasture Dairy Herd Improvement (DHI). In their respective pastures, all animals have the same growth, feeding environment, and similar production levels (305 days correction of milk yield, 305D).
Blood samples were aseptically collected from the bovine coccygeal vein to a disposable (EDTA-K2) vacuum anticoagulation tube. And blood samples were collected from the case group during the presence of clinical mastitis based on SCC and clinical findings. Genomic DNA was extracted from whole blood using TIANamp Genomic DNA Kit (Cat#DP304-03, TIANGEN Biotech). The quality of genomic DNA was detected by NanoDrop and agarose gel methods (extracted 3 µL of genomic DNA, loaded on 1% agarose gel, 100 V CV 25 min, viewed under ultraviolet light and photographed). Genomic DNA with a quality (OD 260 /OD 280 = 1.7-1.8) that meets the experimental requirements was used for subsequent 2b-RADseq analysis.
2b-RAD Library, Sequencing, and Raw Reads Quality Control Forty sample libraries set up met a protocol developed by 2b-RAD sequencing needs with a little change and five-label tandem technique (40,48). The Bos Taurus genomes (Bos taurus UMD 3.1.1) was used as the reference for predicting electronic-enzyme-cut digestion of genomic DNA. Finally, Bael restriction endonuclease, a commonly used restriction endonuclease that specifically recognize the AC(N) 4 GTAYC (N NNNNNNNNNACNNNNGTAYCNNNNNNNNNNNN) and GRTAC(N) 4 GT (NNNNNNNNNNNNNNNTGNNNNCATR GNNNNNNN) sites, and has the advantages of high stability, good repeatability, and uniform distribution in the genome, was selected to digest genomic DNA. The restriction enzyme digested DNA fragment tags of each sample were linked by standard 5 ′ -NNN-3 ′ connector. Paired-end sequencing was carried out on the Illumina Hiseq Xten (https://support.illumina.com/ downloads/sequencing-analysis-viewer-software-v2-1-8.html) platform. Construction the library and steps of the modeling followed Wang et al. with a minor change (41); Supplemental statistical model ( Figure S1).

SNPs RAD Typing, Linkage Disequilibrium (LD) and Genetic Diversity Analysis
Single nucleotide polymorphism marker typing (RAD typing) was performed on Enzyme Reads using the maximum likelihood (ML) method in SOAP software. The statistic SNP typing results, using R Package cluster analysis of the differences between sample SNPs. SNPs were annotated using SnpEff software (Version: 4.1g) (http://snpeff.sourceforge.net/). Plink software was used to calculate the r 2 value of the pairwise SNPs, with the main parameters set to: -r 2 -Id-window-kb 1000-Id-window 50-Id-window-r 2 0.2. According to the median of the software, the work F (x) = 1/(log 10 ((x+10 (7−C) )/10 7 ) + C) was used to fit, and map by chromosome/grouping. To find the LD block in the casecontrol group, we added the parameters "block output GAB-pair wise Tagging' to R package runner.
Principal component analysis (PCA) method was used to evaluate the population. Then a correlation test for each subgroup, including the first five PC selected as covariate analysis signs for population division adjustment (Table S2). And the first two influential feature vectors were selected to draw the correlation between the samples. To assess experimental samples' genetic diversity, polymorphic information content (PIC), observed heterozygosity (Ho) and expected heterozygosity (He) values were also calculated for each SNP loci. In addition, the genetic differentiation coefficient between the subgroups was statistical.

Association Analysis Between SNPs and Clinical Mastitis Traits in Chinese Holstein
We considered a GWAS for the quality traits of dairy cow mastitis, and 2b-RAD markers genotype for each individual SNP loci. To ensure the accuracy of the analysis, we used multistage GWAS for identification of important SNPs associated with mastitis traits. GWAS analysis was performed using Bayesian and Logistic regression analysis model to compare significant SNPs between case-control. Quantile-Quantile Plot (QQ-plot) evaluated the rationality of the two analysis models. GO enrichment analysis performed on all genes with SNPs, and their functions described in conjunction with GO annotations. Hypergeometric Distribution Test (Cytoscape software 3.6) used to calculate significant gene enrichment in each GO entry. Protein and protein interaction networks were predicted and constructed in the string 10.5 database (Bos Taurus) (49).

Independent Sample Population Validated Significant SNPs
The independent sample population validation had the same screening conditions as the 2b-RADseq stage. The casecontrol study was used for the population validation of significant SNPs. The genomic DNA from peripheral blood of cows was extracted by Rapid Blood Genomic DNA Isolation Kit (#B518223, Sangon Biotech). The whole blood genomic DNA quantified (OD260/OD280=1.7-1.8) by the SparkTM multimode microplate reader (Tecan, Swizerland).

Isolation and Primary Culture of Bovine Mammary Epithelial Cells (bMECs) From Control Group
The steps of isolation of bMECs were as follows (50,51): (1) Three healthy Chinese Holstein cows in the third earlylactation stage were selected and slaughtered. Then breast tissues (including mammary acinar tissue) were collected and cleaned with 75% ethanol for two times (5 min/times), and placed in Dulbecco's Modified Eagle's Medium/Nutrient Mixture F-12 Ham (DMEM/F12) complete medium supplemented with 15% Fetal Bovine Serum (FBS, Lot#42F1376K, Gibco, USA), 1% Penicillin-Streptomycin Solution (Double-antibiotics, Cat#C0222, Beyotime, China), 1% Hydrocortisone (Cat#CS-2226, MedChemExpress, USA), and 1% Estradiol (Cat#HY-B0141, MedChemExpress, USA). The breast tissues were brought back to the aseptic processing room within 2 h. (2) Rinse the breast tissue with 1 × PBS buffer (containing double-antibiotics) under sterile conditions until clarified. (3) The mammary acinar tissues were isolated under aseptic conditions, and washed with DMEM/F12 medium until clarification, then centrifuged at 800 rpm for 2 min. (4) Mammary acinar tissues were cut into 1-2 mm 3 tissue pieces with sterile scissors, rinsed to clear with DMEM/F12 medium, and centrifuged at 800 rpm for 2 min. (5) The tissue pieces were uniformly inoculated into a cell culture flask, placed in a 37 • C, 5% CO 2 incubator, inverted for 30-60 min, and the culture solution was removed. (6) The cell culture flask was placed upright and DMEM/F12 complete medium was added. The culture medium was changed every 48 h, and the bMECs were purified when the cell adhering density reached 90%. (7) Primary culture of bMECs: The purified cell suspension was inoculated into new culture flask, placed in 37 • C, 5% CO 2 cell incubator, and the culture medium was changed every 48 h. (8) Identification of bMECs: Epithelial cells are sensitive to Cytokeratin-18 (C-04, Cat#sc-51582, SANTA CRUZ, USA), and immunofluorescence staining is used to identify bMECs ( Figure S7B). bMECs were sub-cultured (primary cells within 20 generations).

Statistical Model Analysis
Linear models are commonly methods for analyzing phenotypic and genotype correlations. Strict quality control is used to remove SNP sites that are not performing well in RAD typing. In this study, Bayesian and Logistic regression models were used to detect significant SNP loci in clinical mastitis traits in dairy cows (Supplemental statistical models). Case-control study validated the population genetics of important SNP loci associated with mastitis traits in dairy cows. Differential  expressions of genes relative quantitative Analysis (2 − Ct ) are based on qRT-PCR results via Student's t-test and one-way ANOVA (and non-parametric).

High Quality SNPs Acquisition, RAD Typing and Population Genetic Differentiation
In this project, the genomic DNAs of 40 Chinese Holstein cows digested with Bael restriction endonuclease and the total of 7,957,920 unique tags were obtained. On Clean Reads, the proportion of A/C/G/T/N bases appeared in each position were shown in Figures 1A,B was quality values of the base sequencing results at each position and the base recognition accuracy rate were 99.99%. The average distance between adjacent tags was 9,589 bp ( Figure 1C), and the unique tag alignment ratio for all samples was 59.69%−72.71%. 10,058 high quality SNPs from 7,957,920 unique tags were selected for RAD typing, and the distribution of SNPs on chromosomes was analyzed by sliding window (Figure S2). All SNPs and their average sequencing depth (17.43×) were shown in Figure 1D (Bayesian analysis) and Figure 1E (logistic regression analysis). Differential SNP cluster analysis showed that data pattern similarity within group was higher; while the similarity between groups was lower. The genetic relationship between the samples was also shown in Figure 1F. PCA was performed on the obtained SNP markers to obtain the two most influential feature vectors ( Figure 1G). The LD coefficients between pairwise markers of SNPs in the genome were calculated, and average LD coefficients between the molecular markers of the case and control group were shown in Figure 1H, respectively. Their LD decay rates also appeared to be the same, and C value were 3.098 (case) and 3.081 (control). Here, since heterozygote risk assessment was between two homozygotes, this line fit the data reasonably which matched to additive genotype risk. In this case there was no deviation, and the test was convincing ( Figure S3A, Table S3). We also calculated the genetic differentiation coefficient (Wright's fixation index, Fst) between the two groups. The value was 0.01869, indicating that the genetic differentiation between case-control groups was smaller.

Genome Wide Association Analysis
The Quantile-Quantile Plot (QQ-plot) data showed that the P-value was consistent with expected values at all SNP site, indicated that the two analysis models matched well (Figures S3B,C). Bayesian analysis has identified 42 significant SNPs with P < 0.001 (Table S4), while logistic regression analysis model has selected 51 SNPs with P < 0.01 ( Table S5). As expected, significant SNPs screened by two analytical models would vary. Total of 27 significant SNPs appeared simultaneously in both analytical models, except these loci with odds ratio (OR) and 95% Confidence Interval (CI) value were not "Na, " and Upper bounds of 95% confidence interval (U95) value were not "Nan" ( Table 1). We also noticed that the OR values of the eight SNPs (rs114843903, rs5881560, rs17514753, rs17518215, rs22015301, rs9704351, rs20438858, and rs5088452) were <1, indicating that these SNPs are associated with mastitis susceptibility traits. However, the OR values of other 19 SNPs, including rs75762330 and rs88640083, were >1, indicating that these SNPs were associated with mastitis resistance.
We annotated all 27 significant SNPs to determine their location in the chromosomal genome (Figure 2A, Table 2):14 SNPs located in the intergenic region, 10 in intron, and one in 3 ′ -UTR, one in upstream and one in downstream. Go annotation for 27 significant SNPs revealed that only three SNPs were annotated to immune-related genes, namely SNPs rs75762330, rs88640083, and rs20438858. SNP rs75762330 (C > T, OR > 1, PIC = 0.2999 > 0.25) was annotated to PTK2B gene located on BTA 8 with moderate polymorphism. SNP rs88640083 (A > G, OR > 1, PIC = 0.1676 < 0.25) annotated to SYK gene located on BTA eight with low polymorphism. SNP rs20438858 (G > A, OR < 1, PIC = 0.3366 > 0.25) was annotated to TNFRSF21 located on BTA 23 and was moderate polymorphism. The Fst values of the rs75762330, rs88640083, and rs20438858 were 0.2165, 0.1297 and 0.2459, respectively (Table S9), suggesting that the genetic differentiation of the three SNPs in case-control group was small.

Independent Sample Population Verification Three Significant SNPs
Correlation analyses were performed on three significant SNPs in another larger independent Chinese Holstein dairy population via direct sequencing (Figure 3). The results showed that the p-values of the three SNPs were all <0.05, indicating that they have significant association with Chinese Holstein mastitis traits ( Table 4). The correlation between rs20438858 and risk of mastitis was still statistically significant, with the adjustment allele OR = 0.359 < 1, suggesting that rs20430083 is associated with mastitis susceptibility traits. While, the other two significant SNPs (rs75762330 and rs88640083) with the adjustment allele OR = 2.416 and 1.879 > 1, respectively, were associated with mastitis resistance traits. Population verification data were consistent with the GWAS analysis results. We also noted that AFe (Attributable Fraction) value for rs75762330 and rs88640083 was 0.2489 and 0.2426 > 0, respectively, while rs20438858 AFe value was −1.786 < 0, also indicating that rs75762330 and rs88640083 are associated with mastitis resistance traits, and rs20430083 is associated with mastitis susceptibility traits.

Network Analysis of Three Immune-Regulated Genes and Their Proteins
Gene ontology enrichment analysis (biological process, cellular component, and molecular function) indicated that three important genes are associated with adaptive or innate immune response in Chinese Holstein (Figure 4 and Figure S4). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis showed that seven important pathways related to immune response were detected (False Discovery Rate (FDR) < 0.05, Table 3), suggesting that multiple immune response pathways are involved in Chinese Holstein mastitis ( Figure S5A). Network protein interaction analysis showed that SYK and PTK2B are involved in natural killer cell mediated cytotoxicity pathway. SYK and AKT1 are in the PI3K-AKT and B cell receptor signaling pathways. PTK2B and AKT1 are in the chemokine signaling pathway. Co-expression scores of the above two proteins are between 0.074 and 0.314 ( Figure S5B). TNFRSF21 is involved in the regulation of cytokine-cytokine receptor interaction (FDR = 1.48E-3).
Expression Variation of PTK2B, SYK, TNFRSF21, TLR4, AKT1, NF-κB, IL-1β, IL-8, and IL-10 in Case-Control Group Total of nine genes (three candidate genes: SYK, PTK2B, and TNFRSF21; TLR4, Toll-like receptor; AKT1, Protein kinase B; NF-κB, Nuclear factor-kappa B; IL-8 and IL-1β, pro-inflammatory factor; IL-10, inflammation and immunosuppressive factor) were selected for qRT-PCR analsysis (Figure 5). It showed that expression of SYK and PTK2B was down-regulated significantly (P < 0.001) in mastitis cows in comparison to the control group, while TNFRSF21 was up-regulated. AKT1 expression level had no significant difference (P > 0.05), while NF-κB down-regulated (P < 0.05). TLR4 (P < 0.05) and two other pro-inflammatory factors (IL-8, FIGURE 2 | Significant SNP association analysis and chromosomal localization in Chinese Holstein mastitis. (A) Manhattan map showed significant SNPs associated with Chinese Holstein mastitis traits (co-marked by both models: Bayesian (left), Logistic analysis (right). Triangular markers are significant SNPs that appear in both analytical models (red triangles are SNPs annotated to immune-related genes: rs75762330, rs88640083, and rs20438858). (B) Partial LD block of three significant SNPs, respectively, with a distance interval of 1 Mb. Important SNPs associated with potential candidate functional genes. The redder the LD blocks color, the stronger correlation (black circles). P > 0.05; IL-1β, P < 0.01) were also up-regulated. However, IL-8 expression level had no significant difference in peripheral white blood cells (P > 0.05). The expression level of IL-10 were up-regulated (P < 0.05). The mean C t value with SD of the β-actin in case-control samples were 20 ± 0.2665 and 20.43 ± 0.2789, respectively ( Figure S7A).
Expression Variation of PTK2B, SYK, TNFRSF21, TLR4, AKT1, NF-κB, IL-8, IL-1β, and IL-10 in bMECs via LPS Stimulation The expression level of candidate genes and related proinflammatory factors was detected in LPS-stimulated bMECs in vitro. The optimal LPS concentration was 25 µg/ml ( Figure 6A). We then performed the time course experiment with this concentration (Figures 6B-I). Figure 6B showed that SYK expression was down-regulated, and then rebounded from 30 to 120 min. PTK2B expression level was down-regulated at 0-45 min, and then up from 45 to 120 min ( Figure 6C). TNFRSF21 expression level was up-regulated from 15 to 60 min, and then down during 60-120 min ( Figure 6D). The expression of TLR4 was down-regulated in 0-45 min and up-regulated in 60-120 min ( Figure 6E). IL-1β expression levels were downregulated within 0-60 min and up-regulated at the 120 min ( Figure 6F). IL-8 expression up-regulated and had an upward trend, which was almost the same at 60-120 min (Figure 6G). AKT1 expression was up-regulated, and showed an upward trend at 0-30 min and a downward trend at 30-120 min (Figure 6H). NF-κB expression was down-regulated within 0-45 min, while up-regulated after 60 min ( Figure 6I). The expression of IL-10 was up-regulated, and showed an upward trend at 0-15 and 30-45 min and a downward trend at 15-30 and 45-120 min ( Figure S6). The mean C t value with SD of the β-actin in LPS concentration stimulated and time course were shown in Figures S7C,D, respectively.

DISCUSSION
Clinical mastitis is a common unavoidable disease that seriously affects the health benefits and breeding profits of dairy cows (3,38). In the production process, mastitis is still not fully controlled despite thorough hygiene, antibiotic treatment, vaccination, and other measures (52). Therefore, the association between dairy cows' genetic variation and mastitis is receiving more and more attention (18,35,53,54). This study sought to screen for significant differential SNPs associated with clinical mastitis via GWAS-2bRADseq, to identify potential clinical mastitis susceptibility or resistance genes and to verify for their authenticity and role in mastitis resistance/tolerance in vitro and in vivo.
Two-stage correlation analysis for three mastitis significant SNPs and Population genetic verification.
In stage I, a total of 10,058 high quality SNPs were screened by 2b-RADseq. The average LD coefficient of 100 kb on the genome of case-control two Chinese Holstein cows was about 0.5. However, the corresponding LD coefficient was still above 0.3 when the distance was 1,000 kb. It is suggested that the genetic diversity of the population of Chinese Holstein cows decreased after acclimatization and selection, and the linkage degree (LD) between the loci were strengthened. Bayesian model screened out 42 significant SNPs when P < 0.001, while the Logical regression analysis model had no SNPs; however, the Bayesian model had 257 significant SNPs when P < 0.01, while Logical regression analysis model labeled 51 significant SNPs. Therefore, the accuracy of molecular breeding value (MBV) prediction is important for the successful application of genome selection (GS) (55). Therefore, selecting the appropriate analytical model based on genotypic and phenotypic data is critical to achieving accurate results (56). This study, two association analysis models (Bayesian model and logical regression model) were used to process 2b-RAD sequencing data to improve accuracy and reduce false positives. It's obvious that the thresholds values of the two analytical models were significantly different. Giving us hints that when processing data, it must be statistically analyzed according to their respective thresholds. Comparison of results from two analytical models, we reduced the dimensions when processing data and only considered the SNPs associated with immunerelated genes and their signaling pathways. Finally, we identified three (rs75762330, rs88640083, and rs20438858) novel bovine mastitis traits significant SNPs in Chinese Holstein cows. With regarding to stage II, we used a case-control study to verify the association between three significant SNPs and bovine mastitis. We compared the exposure ratios of important SNPs in case and control groups, and only considered the relationship between SNP and mastitis traits under the exclusion of external matching factors (57)(58)(59). According to the Pitman efficiency increment formula (2R/(R+1)), we determined the appropriate sample size and gained higher test efficiency (1/4). Here, we validated the association of three important SNPs with cows mastitis. SNP rs75762330 and rs88640083 (OR > 1, AFe > 0) are positive, and SNP rs20438858 is negative (OR < 1, AFe < 0) associated with mastitis susceptibility, respectively. Three significant mastitis-associated SNPs are located in genomic non-coding sequences and in low or moderate genetic polymorphisms.
Previous researchers found that conserved non-coding regions (CNCs) in introns and near genes show large allelic frequency shifts, similar in magnitude to missense variations, suggesting that CNCs are critical for gene function regulation and evolution in many species, including yeast, fruit flies, and vertebrates (60-64). Our GWAS data provided a statistical list of SNPs associated with mastitis traits in dairy cows, where three significant SNPs (rs75762330, rs88640083, and rs20438858) are annotated to immune-related genes are located in noncoding regions (intron and intergenic) of the genome. GWAS studies also showed that the mastitis traits were controlled by multiple loci in the genome, and the genetic effects of each locus were relatively small (18,26,34,38). SNP rs88640083 (PIC = 0.1676 < 0.25) had low polymorphism, while SNPs rs75762330 (0.25 < PIC = 0.2999 < 0.5) and rs20438858 (0.25 < PIC = 0.3366 < 0.5) were moderately polymorphic. Despite the low or moderately heritability of clinical mastitis, identifying significant SNPs associated with clinical mastitis are extremely important for breeding programs that reduce the incidence of mastitis (36,65).

LPS Stimulates bMEC in a Dose-and Time-Dependent Manner in vitro
Pathogen-specific mastitis traits are direct indicators of cow mastitis infection. E. coli is the main pathogen of clinical mastitis, while Staphylococcus aureus and Streptococcus uberis mainly cause sub-clinical mastitis (8,66). Once infected, mammary gland innate immunity is initiated; bMECs act as the first line of defense against exogenous pathogens and secrete several cytokines (e.g., IL-1β, IL-8, and TNF-α) to recruit neutrophils into the mammary gland (8,67,68), causing a sharp rise in SCC intra-mammary infection in early-and late-lactation (69)(70)(71) and the corresponding mastitis phenotype. Therefore, bMEC is the key for rapid elimination of bacteria and prevention of mastitis (inflammation of bovine mammary gland) (72). Exogenous stimulation of bMEC activity was dose-and time-dependent manner as well (73)(74)(75). SYK is a non-receptor tyrosine kinase that plays an essential role in various biological functions (76). SYK-Pyk2 (PTK2B) axis is participated in regulation of TNF activation in human neutrophils (77). Therefore, the optimal LPS stimulation concentration and time was determined with reference to the SYK expression level in bMECs. Then we counted and compared the changes of TLR4, IL-1β, IL-8, and NF-κB expression in LPS-stimulated bMECs time course, and their expected parameters and fold differences serve as positive controls (67), and confirming the reliability of time course bMECs study results: the down-regulation of PTKTB2 and SYK, respectively, and the up-regulation of TNFRSF21.

With Mastitis Traits in Chinese Holstein Cows
Gene Ontology function analysis annotated three significant SNPs to three immune regulatory genes PTK2B, SYK, and TNFRSF21, respectively. The LD block data also showed that there were strong LD relationship between rs75762330 and PTK2B, rs88640083 and SYK, rs20438858 and TNFRSF21, suggesting that these three immune-regulated genes are important candidate genes associated with mastitis traits in Chinese Holstein cows. PTK2B is a protein tyrosine kinase that regulates humoral immune response and homeostasis of innate immune cells (78)(79)(80)(81). It is involved in regulating TLR4 cascade and IL-10 production in macrophages and affects the migration of dendritic cells (DCs) (78,80,82). PTK2B participates in the innate immune response through interaction with IkappaB kinase β (IKKβ) and TANK-binding kinase 1 (TBK 1) (83), and is critical for neutrophil degranulation and host defense against bacterial infection (84). PTK2B expression levels in bPBLs of mastitis cows were significantly down-regulated (P < 0.001) and had a time course difference in LPS-stimulated bEMCs, suggesting that PTK2B plays an important role in clinical mastitis.
SYK is an important regulator factor for TLR4 signaling pathway (85,86), and regulates the secretion of IL-1β through CARD 9 in dendritic cells (87). In macrophages, SYK, as a downstream regulatory molecule of TLR4 and IL10, plays a dual role in pro-inflammatory and anti-inflammatory responses (88,89). SYK is also involved in regulating the proliferation of dairy mammary epithelial cells, milking cycles, and milk production (90). The expression levels of SYK in bPBLs and LPS-stimulated bMECs of mastitis cows were significantly downregulated, indicating that SYK is involved in inflammatory and immune response of clinical mastitis. As shown in Figures 6B,C, the temporal sequences of SYK, PTK2B, TLR4, and IL-1β expression are similar, suggesting that there might be an immune response pathway in bMECs participating in bovine mastitis inflammatory response. However, as the stimulation is prolonged, TLR4, NF-κB, and IL1β are up-regulation at later time points, suggesting that different immune response mechanisms may be initiated, and the phenomenon which needs further study.
TNFRSF21 is a member of the TNF/TNFR family, and plays a critical role in immune response and inflammation (91)(92)(93). It regulates the degeneration of the mammary gland and providing protection against infection (94). TNFRSF21 is also involves in NF-κB signaling pathway (95), which is the key pathways to identify exogenous pathogens and induce inflammation and immune response. TNFRSF21 was significantly up-regulated in bPBLs of mastitis cows and LPS-stimulated bMECs, indicating that TNFRSF21 plays an important role in bovine mastitis inflammation (96).
Three Immune-Regulated Genes (PTK2B, SYK, and TNFRSF21) Initiate Different Immune Response in bPBLs and bMECs to Cope With Mastitis In bPBLs and bMECs, PTK2B, SYK, and TNFRSF21 are involved in different mechanisms of immune responses associated with mastitis (Figure 7). There was no significant difference in the expression level of AKT1 in bPBLs. But interestingly, they were significantly up-regulated in bMECs, and the trends were completely opposite to SYK and PTK2B, suggesting that PTK2B, SYK, and AKT1 participate in the immune cascade response of bMECs. PTK2B-AKT1 is involved in chemokine signaling pathway, while SYK-AKT1 is involved in B cell receptor signaling pathway and PI3K-AKT signaling pathway. However, AKT1, as a key player in the Jak2/STAT5 pathway, plays an important role in the regulation of differentiation, secretion, survival, and proliferation of mammary epithelial cells as well as mammary remodeling and lactation sustainability in dairy cows (90,(97)(98)(99)(100). JAK2 and STAT5, two key players in the immune response, are associated with cow mastitis susceptibility (11,101). Mastitis caused by bacterial infection of the mammary gland induces IL8 expression (102), and IL-8 plays an important role in the recruitment and degranulation of neutrophils in mammary gland (13). However, there were no significantly different in the expression levels of IL-8 in bPBLs of Chinese Holstein, while there were significant differences in LPS-stimulated bMECs, indicating that IL-8 has vital function on innate immunity of bMECs. The expression of IL-1β was significantly upregulated in both clinical mastitis bPBLs and LPS-stimulated bMECs, suggesting that it is essential in bovine mastitis. IL-10, an anti-inflammatory cytokine involved in adaptive and innate immune regulation, plays a critical role in infection by limiting excess immune response to the pathogens and preventing host damage (103-106). It's regulated by Pyk2 (PTK2B) in LPS-stimulated mouse macrophages, and negatively regulated by SYK in dendritic cells (82,89). Therefore, upregulation of IL-10 expression in bPBLs, suggesting that IL-10 may be regulated by PTK2B and SYK, and thus participate in the regulation of inflammatory response in dairy cow mastitis. NF-κB pathway plays a central role in the host innate and acquired immune responses to microbial pathogen infection and is essential for maintaining immune homeostasis (107,108). NF-κB expression was down-regulated in mastitis cows' bPBLs and time-dependent in LPS-stimulated bMECs. Therefore, we hypothesized that SYK and TNFRSF21 may be involved in the NF-κB pathway, stimulating IL-1β secretion and promote mastitis inflammation in both bPBLs and bMECs. SYK and PTK2B may be involved in the SYK/PTK2B/AKT1/JAK2 pathway in bMECs that stimulates IL-8 secretion and recruits neutrophils to kill pathogenic microorganisms, resulting in the resistance to mastitis inflammation.

CONCLUSIONS
In this study, we were committed to improve understanding of biogenetic variation of mastitis, and providing a genetic basis for constructing and improving anti-mastitis characteristic in Chinese Holstein cows. First, we used GWAS-2b-RADseq to identify significant SNPs. We used two-stage correlation analysis to find 27 significant SNPs associated with risk of mastitis, among which three SNPs (rs75762330, rs88640083, and rs20438858) were annotated to immune-regulated genes (PTK2B, SYK, and TNFRSF21). Next, population genetic verification of three important SNPs was performed in another individual cow population, confirming that SNPs rs75762330, rs88640083, and rs20438858 were associated with mastitis susceptibility. Transcriptional analysis showed that PTK2B and SYK expression level was down-regulated in both bPBLs of clinical mastitis cows and in vitro LPS (E. coli)-stimulated bMECs indicating that PTK2B and SYK are important candidate genes associated with mastitis resistance traits. TNFRSF21 was up-regulated suggesting that TNFRSF21 is an important candidate gene associated with mastitis susceptibility traits. Network analysis showed that different immune response mechanisms are activated to deal with pathogen (E. coli) invasion. However, their molecular mechanisms for response to bovine mastitis remain unclear. Do SYK and TNFRSF21 participate in the NF-κB pathway? What are their statuses in bMECs? What are the co-expression mechanisms of SYK and PTK2B in bMECs in response to pathogens? Is there a SYK/PTK2B/AKT1/JAK2 pathway in bMECs or bovine mammary gland? Answering these questions will keep scientists interested in the immune response mechanisms caused by bovine mastitis.

DATA AVAILABILITY
The publicly available datasets generated for this study can be found in the NCBI SRA database -PRJNA556499 -https://www. ncbi.nlm.nih.gov/sra/PRJNA556499.

ETHICS STATEMENT
The procedures involving animals were approved by the Animal Welfare Committee of Nanjing Agricultural University, Nanjing, China, and approval No.20160615. And all animal experiments were carried out in strict accordance with the guidelines and rules established by the committee.

ACKNOWLEDGMENTS
We thank Nanjing Weigang Dairy Co., Ltd. for providing experimental Chinese Holstein blood samples and Shanghai Oe Biotech Co., Ltd. for providing 2b-RAD genome sequencing technology support.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu. 2019.01579/full#supplementary-material   (1) Genotyping and genotype score of 10058 SNPs (A): 0 (two bases of the type were different from the reference genome); 1 (two bases of typing were the same as one of the reference genomes); 2 (two bases of typing were the same as those of the reference genome). (2) The consistency of Bayesian (B) and logistic regression (C) analysis for SNPs observed and expected value -log10 (P), respectively. The P-value observation is almost the same as the expected, indicated that the analysis model was reasonable. There were several SNPs P-values exceeded expected, which suggested that these locus might be significantly associated with dairy cows' mastitis traits.          Table S8 | Primers of internal reference gene (β-actin), three immune-regulate genes, TLR4, AKT1, NF-κB, and three interleukin genes for RT-qPCR analysis. Data Sheet 1 | A supplement to the "Materials and Methods" section.