• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    Rumen microbiome structure and metabolites activity in dairy cows with clinical and subclinical mastitis

    2021-09-19 13:14:16YueWangXuemeiNanYiguangZhaoLinshuJiangMenglingWangHuiWangFanZhangFuguangXueDengkeHuaJunLiuJunhuYaoandBenhaiXiong

    Yue Wang ,Xuemei Nan ,Yiguang Zhao ,Linshu Jiang ,Mengling Wang ,Hui Wang ,Fan Zhang,Fuguang Xue,4,Dengke Hua,Jun Liu,Junhu Yao and Benhai Xiong*

    Abstract Background:Due to the high prevalence and complex etiology,bovine mastitis (BM) is one of the most important diseases to compromise dairy cow health and milk quality.The shift in milk compositions has been widely investigated during mastitis,but recent studies suggested that gastrointestinal microorganism also has a crucial effect on the inflammation of other peripheral tissues and organs,including the mammary gland.However,research focused on the variation of rumen inner-environment during mastitis is still limited.Therefore,the ruminal microbial profiles,metabolites,and milk compositions in cows with different udder health conditions were compared in the present study.Furthermore,the correlations between udder health status and ruminal conditions were investigated.Based on the somatic cell counts (SCC),California mastitis test (CMT) parameters and clinical symptoms of mastitis,60 lactating Holstein dairy cows with similar body conditions (excepted for the udder health condition) were randomly divided into 3 groups(n=20 per group) including the healthy (H) group,the subclinical mastitis (SM) group and the clinical mastitis (CM) group.Lactation performance and rumen fermentation parameters were recorded.And rumen microbiota and metabolites were also analyzed via 16S rRNA amplicon sequencing and untargeted metabolomics,respectively.Results: As the degree of mastitis increased,rumen lactic acid(LA) (P<0.01),acetate,propionate,butyrate,valerate(P<0.001),and total volatile fatty acids (TVFAs) (P<0.01) concentrations were significantly decreased.In the rumen of CM cows,the significantly increased bacteria related to intestinal and oral inflammation,such as Lachnospiraceae(FDR-adjusted P=0.039),Moraxella (FDR-adjusted P=0.011) and Neisseriaceae (FDR-adjusted P=0.036),etc.,were accompanied by a significant increase in 12-oxo-20-dihydroxy-leukotriene B4 (FDR-adjusted P=5.97×10?9) and 10beta-hydroxy-6beta-isobutyrylfuranoeremophilane(FDR-adjusted P=3.88×10?10).Meanwhile,in the rumen of SM cows,the Ruminiclostridium_9 (FDR-adjusted P=0.042) and Enterorhabdus (FDR-adjusted P=0.043)were increased along with increasing methenamine (FDR-adjusted P=6.95×10?6),5-hydroxymethyl-2-furancarboxaldehyde (5-HMF) (FDR-adjusted P=2.02×10?6) and 6-methoxymellein (FDR-adjusted P=2.57×10?5).The short-chain fatty acids (SCFAs)-producing bacteria and probiotics in rumen,including Prevoterotoella_1 (FDRadjusted P=0.045)and Bifidobacterium (FDR-adjusted P=0.035),etc.,were significantly reduced,with decreasing 2-phenylbutyric acid (2-PBA) (FDR-adjusted P=4.37×10?6).Conclusion: The results indicated that there was a significant shift in the ruminal microflora and metabolites associated with inflammation and immune responses during CM.Moreover,in the rumen of cows affected by SM,the relative abundance of several opportunistic pathogens and the level of metabolites which could produce antibacterial compounds or had a competitive inhibitory effect were all increased.

    Keywords:Dairy cows,Lactation performance,Mastitis,Rumen fermentation,Ruminal metabolisms,Ruminal microbiota

    Background

    Bovine mastitis (BM) is commonly recognized as an intramammary infection (IMI) caused by pathogens.The main manifestations of clinical mastitis (CM) are increased body temperature,udder redness and pain,inflammatory cell infiltration,acinar duct edema,as well as interstitial bleeding,which result in elevated milk somatic cell counts (SCC) and decreased milk yield [1].In contrast,the subclinical mastitis (SM) has high concealment and long incubation period.Although the udders of cows with SM do not have visible changes in the appearance,the increase of milk SCC does occur [1].Usually,mastitis treated by antibiotics is easy to relapse,leading to a continuous decline in milk yield and quality.In addition,nipple necrosis accelerates the culling of cows [2].Although SCC has been widely applied to diagnose the udder health,multitudinous factors can affect SCC [3].Therefore,integrated approaches are used to improve the accuracy of IMI diagnosis,for instance,the combination of California mastitis test (CMT),milk pH value,conductivity,and enzyme test etc.[4].

    Notably,recent studies have shown that gastrointestinal microbiota played a vital role in inflammation of tissues outside the gut,such as mammary gland [5–14].First of all,dysbiosis in gastrointestinal microbiota may be one of the causes of BM.Intestinal bacteria may be derived from milk,and in turn affect the microorganism in the udder after the formation of intestinal flora [6].Ma et al.[7] transplanted fecal bacteria from CM cows to sterile mice,and consequently led to mastitis in mice.Secondly,immune regulation and metabolism of endogenous substances are both potential related factors of mastitis,which also associated gastrointestinal microbiome [8].Several studies in human have confirmed the direct or indirect relationships between gastrointestinal microorganisms and breast inflammation and immune regulation[9–11].Besides,a considerable number of symbiotic bacteria distributed in the gastrointestinal tract are closely related to the immune regulation in dairy cows[12].Probiotics can penetrate the intestinal epithelium,and reach the udder tissue through internal circulation[13],acting as an effective immune enhancer [7].Thirdly,diet was also reported to play a vital role in determining the structure of microbiome and bioactivity of bacterial metabolites in the mammary gland.For example,Mediterranean diet can increase the abundance of Lactobacillus,bile acids and bacterial-modified metabolites and decrease reactive oxygen species and pro-inflammatory metabolites in the mammary gland of monkeys [14].The previous studies have triggered the speculation on the relationship between the gastrointestinal environment and the health status of mammary gland,which should be verified by further investigation.

    Unlike monogastric animals,ruminal microorganisms are the largest microbial group in the gastrointestinal tract.Ruminal bacteria fermentation in dairy cows produces large amounts of volatile fatty acids (VFAs),also called short-chain fatty acids (SCFAs),which are not only the energy source,but also the precursors of milk compositions [15].Milk fat,lactose and protein are formed through acetate,butyrate,propionate and microbial proteins,which are produced by microbial fermentation of feed in the rumen [16].Therefore,milk compositions are indirectly derived from the fermentation products of ruminal microorganisms.Therefore,rumen inner-environment has a decisive influence on milk quality.However,the correlation between the changes in rumen microorganisms and metabolites and the udder health of dairy cows is still in short of knowledge.

    The present study hypothesized that the rumen innerenvironment,including ruminal microorganism profile,fermentation parameters and metabolites activities,varied during IMI.Accordingly,we identified the ruminal microorganisms and metabolites in cows with SM,CM and healthy (H) udder by 16S rRNA sequencing and liquid chromatography-mass spectrometry (LC-MS),respectively.The current study probed into the possible shift in the profile of ruminal microflora and metabolites activities during mastitis,thereby improved the understanding of the relationship between rumen innerenvironment and mastitis.

    Materials and methods

    Ethics statement

    All experimental designs and protocols were approved by the Animal Ethics Committee of the Chinese Academy of Agricultural Sciences (Beijing,China;approval number:IAS-2019-6) and were in accordance with the recommendations of the academy’s guidelines for animal research.

    Animals,diets and experimental design

    This study was carried out in a suburban dairy farm in Beijing,China.Total mixed rations (TMR) with a concentrate to forage ratio of 4:6 were offered to the Holstein cows three times a day at 07:00,13:00,and 19:00,respectively.Cows were fed ad libitum and had free access to water.The TMR ingredients and nutritional compositions are shown in Additional file 1:Table S1.In the current study,the udder health status of cows was comprehensively judged according to the degree of inflammation in each quarter sample from each cow based on the clinical manifestations of the udder and the results of the CMT and SCC in milk [1].CMT detection methods and judgment basis are listed in Additional file 1:Table S2.According to the results,60 Holstein cows were selected and divided into 3 groups:20 cows with healthy udder (H group) (SCC <100,000 cells/mL;no clinical symptoms in udders;CMT results were negative);20 cows with SM (SM group) (600,0003000,000 cells/mL,with obvious signs of inflammation in udders,including udder swelling,redness and milk clots,etc.,as well as positive or strongly positive CMT results).The basic information of cows including parity,days in milk,average daily milk yield,clinical symptoms in udders and milk SCC are displayed in Additional file 1:Table S3.

    Milk and rumen fluid sampling

    Throughout the test period,the milk yield of selected cows was recorded for 3 consecutive days after confirmation of udder health status.Milking was performed 3 times a day using an automatic milking system(Ruishengyuan Machinery Assembly Co.,Ltd.,Hebei,China).Milk samples from H group were mixed from 4 quarters,while milk samples from cows with IMI (SM and CM groups) were only collected from inflammatory quarters.Milk samples from each cow were collected 3 times a day in the morning,noon and evening,respectively,with 50 mL each time,and mixed in a ratio of 4:3:3 [17].Each milk sample was added 0.6 mg/mL potassium dichromate as a preservative and were then stored at 4°C for analysis of milk compositions.Milk fat,protein,lactose and SCC were analyzed by a MilkoScan FT2 milk composition analyzer (FOSS,Copenhagen,Denmark) [17].

    Before the morning feeding on the sampling day,a stomach tube (GCYQ-1-A,Guidi scientific instrument Co.,Ltd.,Shanghai,China) and a 200-mL syringe were used to collect rumen fluid samples from each cow.The first tube of rumen fluid was discarded to minimize saliva contamination and a 50-mL of rumen fluid was collected per cow.The pH value of the rumen fluid was immediately measured by a pH meter (8362sc Portable pH meter;Biyuntian biotechnology Co.,Ltd.,Shanghai,China) [18].The remaining rumen fluid samples were poured into sterile test tubes and immediately stored in liquid nitrogen.These samples were then used to detect VFAs (GC-6890 N,Agilent,Palo Alto,USA),ammonia nitrogen (NH3-N) (HBS-1101 ELIASA,Huawei Delang Instrument Co.,Ltd.,Wuxi,China),ruminal fluid urea nitrogen (RUN) (GF-D200 semi-automatic biochemical analyzer;Gaomi Rainbow Biological Co.,Ltd.,Shandong,China),lactic acid (LA) (MultiskanAsc ELIASA,Thermo,USA) [15],microorganisms composition and metabolites.

    DNA extraction,PCR amplification and sequencing

    Microbial community genomic DNA was extracted from rumen fluid samples using the FastDNA? Spin Kit for Soil (MP Biomedicals,California,USA) according to the manufacturer’s instructions.In brief,1 mL rumen fluid,978 μL sodium phosphate buffer and 122 μL methyltransferase (MT) buffer were added to a lysing matrix E tube,and mixed thoroughly.The tube was then centrifuged at 14,000 r/min for 10 min.Subsequently,the supernatant was transferred to a 1.5-mL centrifuge tube with a supplementation of 250 μL PPS.After centrifuging at 14,000 r/min at room temperature for 5 min,the supernatant was discarded.Then 500 μL 5.5 mol/L guanidine isothiocyanate solution was added and the mixture was transferred to the Spin Filter.After the addition of 500 μL salt/ethanol wash solution (SEWSM),the Spin Filter was centrifuged at 14,000 r/min for 1 min and the filtrate was discarded.After a 3-min drying process,100 μL DNA elution solution-ultra pure water(DES),which was preheated at 55°C for 5 min,was added to the Spin Filter.Finally,total DNA was obtained after centrifuging at 14,000 r/min for 2 min [19].The DNA extract was checked on 1% agarose gel,and DNA concentration and purity were determined by a Nano-Drop 2000 UV-vis spectrophotometer (Thermo Scientific,Wilmington,USA).A total of 10 ng DNA were then used in the following sequencing analysis.The hypervariable V3-V4 region of the bacterial 16S rRNA gene was amplified with primer pairs 338F (5′-ACTCCT ACGGGAGGCAGCAG-3′) and 806R (5′-GGACTA CHVGGGTWTCTAAT-3′) by an ABI GeneAmp? 9700 PCR thermocycler (Thermo Fisher Scientific,Waltham,Mass,USA).PCR reactions were performed in triplicate.The PCR product was extracted from 2% agarose gel and purified using an AxyPrep DNA gel extraction kit(Axygen Biosciences,Union City,CA,USA) according to the manufacturer’s instructions and quantified using a Quantus?Fluorometer (Promega,Beijing,China) [19].The purified amplicons were pooled in equimolar and paired-end sequenced (2×300) on an Illumina MiSeq platform (Illumina,San Diego,USA) [20].The raw sequences have been submitted to NCBI Sequence Read Archive (SRA) database (Accession Number:PRJNA669201).

    Processing,annotation,and statistical analysis of sequencing data

    The raw 16S rRNA gene sequencing reads were demultiplexed,quality-filtered by Trimmomatic and merged by FLASH.Operational taxonomic units (OTUs) with a 97%similarity cutoff were clustered using UPARSE(version 7.1,http://drive5.com/uparse/),and the chimeric sequences were identified and removed [21].The taxonomy of each OTU representative sequence was analyzed by RDP Classifier (http://rdp.cme.msu.edu/) against the 16S rRNA database (Silva SSU128) using a confidence threshold of 0.7[21].The alpha diversity (Shannon,Simpson,Ace,Chao index) of the rumen microbiota based on OTU level was analyzed through MOTHUR(version 1.30.2).Beta diversity(principal co-ordinates analysis,PCoA) and other analyses were also conducted through QIIME (version 1.9.1).According to the abundance of microbial community,Kruskal-Wallis H test was used to conduct hypothesis testing for the microflora in the three groups,to evaluate the effects of ruminal bacterial abundance on the udder health status.False discovery rate(FDR)was used to conduct multiple testing for the correction of P-value.Hierarchical cluster analysis(HCA)and intergroup correlation analysis were performed using the OmicShare tools (http://www.omicshare.com/tools).Phylogenetic Investigation of Communities by Reconstruction of Unobserved States(PICRUSt),a technique that uses evolutionary modeling to predict metagenomes from 16S data and a reference genome database,was conducted to perform functional prediction analysis of differential microbiota.The analysis was operated using the PICRUSt software package (https://picrust.github.io/picrust/install.html#install) [21].The accuracy of the prediction was assessed by the value of the nearest sequenced taxon index (NSTI) presenting as mean±standard error of mean (SEM).The NSTI represents the average phylogenetic distance between the sequenced genomes with the closest relationship between OTU and all microorganisms in a sample.Therefore,the smaller of this value is,the more reliable of the prediction result will be[21].

    Metabolite extraction and LC-MS/MS analysis

    A 100-μL rumen liquid sample was measured,and mixed with a 400-μL methanol:water (4:1,v/v) solution to extract the metabolites.The mixture was settled at ?20°C and treated by a high throughput tissue crusher Wonbio-96c (Wanbo Biotechnology Co.,Ltd.,Shanghai,China) at 50 Hz for 6 min,followed by a 30-s vortex,and then ultrasound at 40 kHz for 30 min at 5°C.The samples were then placed at ?20°C for 30 min to precipitate proteins.After centrifugation under 13,000 ×g at 4°C for 15 min,the supernatant was carefully transferred to sample vials for LC-MS/MS analysis.Chromatographic separation of the metabolites was performed on an ExionLC?AD system(AB Sciex,USA) equipped with an ACQUITY UPLC BEH C18 column (100 mm×2.1 mm i.d.,1.7 μm;Waters,Milford,USA).The mobile phases consisted of 0.1% formic acid in water with formic acid (0.1%)(solvent A) and 0.1% formic acid in acetonitrile:isopropanol (1:1,v/v) (solvent B).The solvent gradient processing was in accordance with the study of Ogunade et al.[22].The sample injection volume was 20 μL and the flow rate was set to 0.4 mL/min.The column temperature was maintained at 40°C.During the period of analysis,all samples were stored at 4°C.The UPLC system was coupled to a quadrupole-timeof-flight mass spectrometer (Triple TOF?5600+,AB Sciex,USA) equipped with an electrospray ionization(ESI) source operating in positive mode and negative mode.The optimal conditions were set as follows:source temperature,500°C;curtain gas (CUR),30 psi;both Ion Source GS1 and GS2,50 psi;ion-spray voltage floating (ISVF),?4000 V in negative mode and 5000 V in positive mode,respectively;declustering potential,80 V;a collision energy (CE),20–60 V rolling for MS/MS.Data acquisition was performed with the Data Dependent Acquisition (DDA) mode.The detection was carried out over a mass range of 50–1000 m/z [22].

    Metabolomics data preprocessing and annotation

    After LC-MS/MS analyses,the raw data were imported into the Progenesis QI 2.3(Nonlinear Dynamics,Waters,USA) for peak detection and alignment.The preprocessing results generated a data matrix that consisted of the retention time (RT),mass-to-charge ratio (m/z) values,and peak intensity.Metabolic features detected at least 50% in any set of samples were retained.After filtering,minimum metabolite values were imputed for specific samples,in which the metabolite levels fell below the lower limit of quantitation and each metabolic features were normalized by sum.The internal standard was used for data quality control (QC) (reproducibility).Metabolic features that the relative standard deviation of QC>30%were discarded.Following normalization procedures and imputation,statistical analysis was performed by log transformed data to identify significant differences in metabolite levels among the three groups.Mass spectra of these metabolic features were identified using the accurate mass,MS/MS fragments spectra and isotope ratio difference with searching in reliable biochemical databases as Human metabolome database (HMDB) (http://www.hmdb.ca/) and Metlin database (https://metlin.scripps.edu/).Concretely,the mass tolerance between the measured m/z values and the exact mass of the components of interest was ±10 ppm.Metabolites with MS/MS fragments score above 30 were considered as confidently identified.Otherwise,metabolites had only tentative assignments [23].

    Data analysis of metabolomics

    A multivariate statistical analysis was performed using ropls (R packages) (version 1.6.2) (http://bioconductor.org/packages/release/bioc/html/ropls.html) from Bioconductor on Majorbio Cloud Platform (https://cloud.majorbio.com).Principle component analysis (PCA)using an unsupervised method was applied to obtain an overview of the metabolic data,and general clustering,trends,or outliers were visualized.All the variables of metabolites were scaled to unit-variances prior to conducting the PCA.Orthogonal partial least squares discriminate analysis (OPLS-DA) was used to determine the metabolic changes among the three groups.All of the metabolite variables were scaled to Pareto scaling prior to conducting the OPLS-DA.The model validity was evaluated by the model parameters R2and Q2to avoid the risk of over-fitting.Variable importance in the projection (VIP) was calculated in OPLS-DA model.P-values were estimated with paired Student’s t-test on single dimensional statistical analysis.Significance was declared at VIP value >1 and FDR-adjusted P-value <0.05.Differential metabolites between 2 groups were summarized,and mapped into their biochemical pathways through metabolic enrichment and pathway analysis based on the KEGG database search (http://www.genome.jp/kegg/).These metabolites were classified according to the involved pathways or the performed functions.Scipy (Python packages,https://docs.scipy.org/doc/scipy/) was used to identify significantly enriched pathways using Fisher’s exact test [22].Receiver operator characteristic (ROC) curve of SPSS Statistics 22 (IBM,Chicago,USA) software was used to conduct the ROC analysis,which was used to evaluate the performance of biomarkers.Area under curve(AUC) was used to indicate the accuracy of the prediction.The higher the AUC value was and the closer to the upper left corner of the curve represented,the higher the prediction accuracy would be [23].The Spearman correlation coefficient was calculated using the OmicShare tools,which were used to analyze the correlation between rumen microbes,metabolites and rumen fermentation parameters.The coefficients range from ?1.0 to 1.0,with variables close to 0 being irrelevant,and variables close to 1 or ?1 being strongly correlated.

    Statistical analysis

    Data of the basic information of cows,milk yield and compositions,rumen fermentation variables and alpha diversity indices of the three groups were analyzed using one-way ANOVA and Student’s t-test of SPSS Statistics 22 (IBM,Chicago,USA) software.Significance was declared at P<0.05,and 0.05

    Results

    Lactation and rumen fermentation

    As the degree of mastitis increased,milk yield,milk fat and lactose contents were significantly reduced (P<0.01),while SCC and milk protein were significantly increased (P<0.01).There were no significant differences in pH-values (P=0.214),the concentrations of RUN(P=0.851),isobutyrate (P=0.325),isovalerate (P=0.213) and acetate to propionate ratio (A/P) (P=0.111)in rumen fluid among the 3 groups.However,the LA,TVFAs (P<0.01),acetate,propionate,butyrate and valerate (P<0.001) contents in the CM group were significantly lower than those of the H and SM groups.In addition,the NH3-N concentration in rumen fluid tended to decline as the degree of mastitis increased(P=0.061) (Table 1).

    Rumen microbial diversity

    A total of 3,072,603 high-quality sequences were generated from the 60 rumen fluid samples.We obtained 2973 OTUs ≥97% identity,with the Good’s coverage of H,SM and CM groups all above 98% (98.9%±0.0013,98.7%±0.0020 and 98.6%±0.0011;mean±standard deviation).Pan OTU was used to observe the increase of the total number of OTU with the increasing number of samples.The pan species curve exhibited that the curve gradually flattened as the number of samples increased,indicating that the sample size of this sequencing was sufficient (Additional file 1:Fig.S1).Similarly,the rarefaction curve showed that the rate of increase in OTU number slowed down with the increasing reads per sample,illustrating that the sequencing depth was adequate(Additional file 1:Fig.S2).Alpha diversity indices analysis suggested that there were significant differences in all the indexes at the OTU level,except for the Simpson index (P=0.10).The Shannon (P=0.02),Ace (P=0.01)and Chao (P=0.01) indexes in the SM and CM groups were significantly lower than those in the H group(Table 2).The PCoA based on unweighted Unifrac distance (P=0.003) (Fig.1a),weighted normalized Unifrac distance (P=0.022) (Fig.1b) and OTU classification level showed that ruminal microbiota in the CM group was distinguishable from the H and SM groups.The ruminal microbiota could also be slightly distinguished between H and SM groups but there were still some intersections.

    Table 1 Milk yield,compositions and rumen fermentation parameters in cows with different udder health status

    Rumen microorganism composition and structure

    In total,21 phyla and 356 genera of bacteria were obtained (Additional file 1:Table S4 and S5).Bacteria with relative abundance <1% of all samples were classified as others.The dominating bacteria at phylum level in the 3 groups were Firmicutes (49.69±0.57%),Proteobacteria (39.14±0.42%),Bacteroidetes (4.65±0.28%),Actinobacteria (3.22±0.06%) and Cyanobacteria (1.29±0.04%) (Additional file 1:Fig.S3a).At the genus level,Prevotella_1 (23.89±0.43%),Succiniclasticum (10.74±0.11%),Mollicutes_RF39 (8.77±0.10%),Pseudobutyrivibrio (4.67±0.27%),Succinivibrionaceae_UCG-001 (3.82±0.16%) and Ruminococcaceae_NK4A214_group (3.08±0.25%) were the top 6 in all 3 groups (Additional file 1:Fig.S3b).

    Table 2 Rumen microbial alpha diversity indices in cows with different udder health status

    Hierarchical cluster analysis conducted for taxa at genus level (top 55 species in total relative abundance)further demonstrated the distribution of ruminal bacteria identified in different udder health status.Through average-linkage clustering techniques,the rumen microorganisms in the 60 dairy cows were aggregated into 3 clusters.Consistent with the PCoA results,ruminal bacteria in the CM group gathered together and separated from the H and SM groups.Although the H and the SM groups were generally differentiated,they still had some intersections,such as SM13,SM4,SM20,SM18 and SM9.There were 25 major bacteria in H group,including Ruminococcus_2,Pseudoscardovia,Acetitomaculum,Bacteroidales,Bifidobacterium,etc.Twelve bacteria,including Succinicl asticum,Ruminiclostridium_9,vadinBE97,Prevotellaceae_UCG-001 etc.,were detected in SM group.The CM group contained 17 genera of bacteria,such as Herbaspirillum,Moraxella,Neisseriaceae,etc.(Additional file 1:Fig.S4).

    Fig. 1 Principal co-ordinates analysis(PCoA)of ruminal microbiota based on a unweighted unifrac and b weighted normalized unifrac of OTUs among H,SM and CM cows.H,healthy;SM,subclinical mastitis;CM,clinical mastitis

    Differentially abundant ruminal microbiota

    Kruskal-Wallis H test was used to perform the analysis of significant differences of ruminal bacteria among the 3 groups.The significantly diverse microbes at genus level were assembled in 6 phylum,incorporating in Firmicutes,Bacteroidetes,Proteobacteria,and Cyanobacteria,etc.(Table 3 and Additional file 1:Fig.S5).The ruminal bacteria with high abundance in the CM group were mainly Pseudobutyrivibrio,Gastranaerophilales,Moraxella,Neisseriaceae,etc.In SM group,Ruminiclostridium_9,Bacteroidales_UCG-001 and Enterorhabdus,etc.,were most abundant.Comparing with the other two groups,Lachnospiraceae,Prevotella_1,Mollicutes_RF39,Bifidobacterium,etc.,differed significantly in H group.

    Prediction of rumen microbial function

    To further speculate the contribution of these key flora,we attempted to predict the function of rumen microbiota in H (average NSTI=0.33±0.005),SM (average NSTI=0.33±0.007) and CM (average NSTI=0.35±0.005) groups by PICRUSt.At superclass level,42.97% of the genes participated in metabolism,21.79% belonged to genetic information processing,and 11.55% were involved with environment information processing.No significant difference was detected in the functions of rumen microorganisms among the 3 groups at level 1(Fig.2a).While at class level,a total of 39 KEGG pathways were noted,and 26 of them were different among H,SM and CM samples (FDR-adjusted P<0.05).The abundance of flora related to amino acids coenzymes(FDR-adjusted P=0.008) and vitamins (FDR-adjusted P=0.004),nucleotides (FDR-adjusted P=0.001) and lipid metabolism (FDR-adjusted P=0.014),etc.,were prominent in H group.Conversely,flora correlated with energy metabolism (FDR-adjusted P=0.001),infection diseases (FDR-adjusted P=0.001),immunity system diseases (FDR-adjusted P=0.001),etc.,accounted for a greater proportion in the CM group.In the SM group,the significantly different functional flora was less than H and CM groups,mainly focusing on translation (FDRadjusted P=0.003) and genetic information processing(FDR-adjusted P=0.002) (Fig.2 and Additional file 1:Table S6).At subclass level,among the 268 KEGG characteristics,33 were different in abundance among the 3 groups (FDR-adjusted P<0.05).Flora related to methane metabolism (FDR-adjusted P=0.003),bacterial invasion and infection (FDR-adjusted P=0.04),primary immunodeficiency (FDR-adjusted P=0.037) and lipopolysaccharide biosynthesis (FDR-adjusted P=0.043),etc.,were significantly enriched in the CM group (Fig.2c and Additional file 1:Table S7).

    Correlation between ruminal flora,and rumen fermentation and lactation performance

    According to Spearman correlation coefficients,the rumen microbiota was associated with several rumen fermentation parameters,and milk yield and compositions(Fig.3).Lachnospiraceae_UCG-006 (r=0.409,FDRadjusted P=0.039) was positively correlated with milk protein.Neisseriaceae (r=0.434,FDR-adjusted P=0.008),Moraxella(r=0.402,FDR-adjusted P=0.005)and Gastranaerophilales (r=0.402,FDR-adjusted P=0.004) were positively correlated with milk SCC.Bifidobacterium (r=0.349,FDR-adjusted P=0.036) was positively associated with lactose,while Pasteurellaceae (r=?0.354,FDRadjusted P=0.039),Moraxella (r=?0.360,FDR-adjusted P=0.045) and Gastranaerophilales (r=?0.360,FDRadjusted P=0.043)were negatively associated with lactose.Prevotella_1 (r=0.384,FDR-adjusted P=0.047) was positively correlated with acetate,but Moraxella (r=?0.396,FDR-adjusted P=0.045) and Gastranaerophilales (r=?0.396,FDR-adjusted P=0.044) were negatively correlated to acetate.Herbaspirillum (r=?0.370,FDR-adjusted P=0.043) and Gastranaerophilales (r=?0.346,FDR-adjusted P=0.043) were negatively associated with butyrate.In addition,Pseudobutyrivibrio (r=?0.400,FDR-adjusted P=0.02),Herbaspirillum (r=?0.347,FDR-adjusted P=0.04) and Lachnobacterium (r=?0.333,FDR-adjusted P=0.042) were negatively associated with LA (Additional file 1:Table S8).

    Table 3 The significantly differential bacteria(relative abundance)among cows with different udder health status

    Comparative analysis of ruminal metabolites

    The total ion chromatogram (TIC) of QC samples in positive and negative ion mode showed an overlap in the response intensity and retention time of each chromatographic peak (Additional file 1:Fig.S6),which indicated the high reliability of data quality.To exhibit the overall difference and the degree of the variability of rumen samples among and within the groups,unsupervised PCA was performed under positive and negative ionization modes,respectively (Fig.4).Symbols representing ruminal metabolites in H,SM and CM groups could be basically separated.Further supervised OPLSDA analysis of the data was performed in positive(Additional file 1:Fig.S7a,c and e) and negative ion modes (Additional file 1:Fig.S8a,c and e).Significant differences were detected between CM and H groups,SM and H groups,and CM and SM groups.Model parameter R2Y (cum) were all above 0.7 and Q2(cum) were all above 0.5,indicating a good ability of model prediction.Model evaluation parameters of response permutation testing,Q2,were all below 0,suggesting that there was no overfitting of the model.The results illustrated that there were significant differences in the rumen metabolism under different udder health conditions of dairy cows (Additional file 1:Fig.S7b,d and f and Fig.S8b,d and f).

    Fig. 2 Phylogenetic investigation of communities by reconstruction of unobserved states(PICRUSt)functional prediction analysis for rumen microorganisms from H,SM and CM cows.a Superclass level.b Class level.c Subclass level.H,healthy;SM,subclinical mastitis;CM,clinical mastitis

    Fig. 3 Correlation analysis between the differentially abundant bacteria and parameters of rumen fermentation as well as lactation performance.Each row represents a genus,each column represents a lactation performance or rumen fermentation parameter,and each cell represents a correlation coefficient between a lactation or rumen fermentation parameter and a bacterial genus.Red indicates the positive correlation,while the blue indicates the negative correlation

    Fig. 4 Principal component analysis(PCA)of the rumen metabolome of H,SM and CM cows in a positive and b negative ion mode.H,healthy;SM,subclinical mastitis;CM,clinical mastitis

    Fig. 5 Significantly differential metabolites(top 15)in the rumen of dairy cows between a CM and H cows, b SM and H cows and c CM and SM cows based on variable importance in the projection(VIP)value under merger mode of positive and negative ion.H,healthy;SM,subclinical mastitis;CM,clinical mastitis.Each column represents a sample,each row represents a metabolite,and the color indicates the relative expression level of the metabolites in the samples.The right side is the metabolite VIP bar graph,the length of bar represents the contribution of the metabolites to the difference between the two groups(VIP>1),the color of the bar means the FDR-adjusted P-value (*0.01

    Differential metabolites identification and evaluation

    A total of 642 differential metabolites were selected including 384 and 258 in positive and negative ion mode,respectively,from the 60 rumen fluid samples.Under positive and negative ionization modes,the significantly differential metabolites (VIP>1 and P-value <0.05)among the 3 groups were shown in Additional file 1:Table S9-S14.Significantly differential metabolites (top 15) in the rumen between CM and H groups,SM and H groups,and CM and SM groups are shown in Fig.5a,b and c respectively.Comparing to H group,10betahydroxy-6beta-isobutyrylfuranoeremophilane (FDR-adjusted P<0.001) and 12-oxo-20-dihydroxy-leukotriene B4 (FDR-adjusted P<0.001) in CM group increased 12.23-and 16.36-fold,respectively.Meanwhile,these two in CM group increased 2.88 and 3.64-fold,respectively,compared with SM group (Table 4).Methenamine(FDR-adjusted P<0.001),5-hydroxymethyl-2-furancarboxaldehyde (5-HMF) (FDR-adjusted P<0.001) and 6-methoxymellein (FDR-adjusted P<0.001) in CM and SM groups were both higher than those in H group.Besides,xestoaminol C (FDR-adjusted P=0.016),lentialexin (FDR-adjusted P=0.031) and cinnamic acid (FDRadjusted P=0.042) in SM group were increased compared to H group.Less abundance of 2-phenylbutyric acid (2-PBA) (FDR-adjusted P<0.001) was detected in CM group compared with H group.In CM group,Nacetylcadaverine (FDR-adjusted P=0.020) and (3R,5Z)-5-octene-1,3-diol (FDR-adjusted P=0.011) were reduced compared with SM group (Table 4 and Fig.5).

    Table 5 Metabolic pathway enrichment analysis of significantly differential metabolites among health,subclinical mastitis and clinical mastitis cows

    Fig. 6 Receiver operating characteristic(ROC)curves of metabolites in the prediction and diagnosis of mastitis.ROC curve reflects the relationship between sensitivity and specificity.The X-axis is false positive rate(specificity).The closer the X-axis to zero,the higher the accuracy will be.The Y-axis is true positive rate(sensitivity).The larger the Y-axis is,the better the accuracy is.The area under curve(AUC)is used to indicate the accuracy of prediction.AUC value is between 0 and 1,the higher the AUC value is,the higher the accuracy of prediction is.a 10beta-Hydroxy-6beta-isobutyrylfuranoeremophilane;b 12-Oxo-20-dihydroxy-leukotriene B4;c Methenamine;d 5-Hydroxymethyl-2-furancarboxaldehyde;e 6-Methoxymellein

    Table 4 HMDB compound classification of differential metabolites in cows with different udder health status

    To evaluate the diagnostic value for BM of these key differential metabolites in the rumen,receiver operator characteristic (ROC) analysis was performed,which reflected the relationship between false positive rate (Xaxis) and true positive rate (Y-axis) (Fig.6).In the current study,the AUC values of the 5 key differential metabolites,including 10beta-hydroxy-6beta-isobutyrylfuranoeremophilane,12-Oxo-20-dihydroxy-leukotriene B4,methenamine,5-HMF and 6-methoxymellein,were all above 0.9,and the ROC curves were all close to the upper left corner.This indicated that these 5 differential metabolites in the rumen screened in the present study were of high value to be served as the biomarkers for bovine mastitis.

    Differential metabolites classification and metabolic pathways enrichment analysis

    Compound classification of differential metabolites in the rumen fluid samples were processed based on comparisons with the HMDB.Differentially clustered metabolites belonging to lipid-like molecules,organ heterocyclic compounds,and organic acids and derivatives compounds accounted for the top 3 among the 3 groups(Table 4 and Additional file 1:Fig.S9).

    Fig. 7 Correlation analysis of differential metabolites and index of lactation and rumen fermentation.Each row represents a metabolite,each column represents a lactation or rumen fermentation parameter.Red means positive correlation,while blue means negative correlation

    Fig. 8 Correlation analysis of differential bacteria genus and differential metabolites in rumen samples.Each row represents a metabolite,each column represents a bacterium.Red means positive correlation,while blue means negative correlation

    Furthermore,the enrichment of metabolic pathway analysis showed that,10beta-hydroxy-6beta-isobutyrylfuranoeremophilane was enriched in limonene and pinene degradation pathway (FDR-adjusted P=0.008).12-Oxo-20-dihydroxy-leukotriene B4 was enriched in arachidonic acid metabolism pathway (FDR-adjusted P=0.026).Besides,methenamine was the intermediate product of formaldehyde biosynthesis (FDR-adjusted P=0.036).5-Hydroxymethyl-2-furancarboxaldehyde participated in the furfural degradation (FDR-adjusted P=0.032).6-Methoxymellein was enriched in the biosynthesis of plant secondary metabolites pathway(FDR-adjusted P=0.045).And 2-phenylbutyric acid was enriched in the butanoate metabolism pathway(FDR-adjusted P=0.042) (Table 5).

    Correlations between ruminal metabolites,and lactation and rumen fermentation parameters

    Based on Spearman correlation coefficients and Euclidean distance matrix,we observed significant correlations between significantly differential metabolites,and lactation as well as rumen fermentation parameters(Fig.7).10beta-Hydroxy-6beta-isobutyrylfuranoeremophilane (r=?0.584,FDR-adjusted P<0.001) and 12-Oxo-20-dihydroxy-leukotriene B4 (r=0.601,FDR-adjusted P<0.001) were negatively correlated to acetate,whereas the 2-PBA (r=0.520,FDR-adjusted P<0.001) was positively correlated to acetate.Meanwhile,10beta-hydroxy-6beta-isobutyrylfuranoeremophilane (r=?0.472,FDR-adjusted P=0.046) and 12-oxo-20-dihydroxy-leukotriene B4 (r=?0.504,FDRadjusted P=0.002) were also negatively associated with butyrate.Besides,(3R,5Z)-5-octene-1,3-diol (r=0.612,FDR-adjusted P<0.001) and N-acetylcadaverine (r=0.670,FDR-adjusted P<0.001) were positively associated with LA.10beta-Hydroxy-6beta-isobutyrylfuranoeremophilane (r=?0.700,FDR-adjusted P<0.001),12-Oxo-20-dihydroxy-leukotriene B4 (r=?0.688,FDR-adjusted P<0.001),5-HMF (r=?0.537,FDR-adjusted P<0.001) and 6-methoxymellein (r=?0.542,FDR-adjusted P<0.001)were negatively associated with lactose,while 2-PBA was positively associated with lactose (r=0.547,FDR-adjusted P<0.001).Furthermore,5-HMF was negatively associated with milk fat (r=?0.604,FDR-adjusted P<0.001).(3R,5Z)-5-Octene-1,3-diol (r=?0.606,FDR adjusted P<0.001) and N-acetylcadaverine (r=?0.582,FDR-adjusted P<0.001) were negatively correlated to milk protein.10beta-Hydroxy-6beta-isobutyrylfuranoeremophilane (r=0.716,FDR-adjusted P<0.001),12-oxo-20-dihydroxy-leukotriene B4 (r=0.712,FDR-adjusted P<0.001),methenamine (r=0.603,FDR-adjusted P<0.001),5-HMF (r=0.550,FDR-adjusted P<0.001) and 6-methoxymellein(r=0.600,FDR adjusted P<0.001) were positively associated with milk SCC,but 2-PBA was negatively associated with milk SCC (r=?0.573,FDR-adjusted P<0.001)(Additional file 1:Table S15).

    Associations between BM-linked rumen microbes and metabolites

    The Spearman correlation coefficients were calculated to reflect the correlations between significantly differential bacteria and metabolites via heat map (Fig.8).In CM group,Lachnospiraceae_UCG-006 (r=0.493,FDRadjusted P=0.038),Pasteurellaceae (r=0.565,FDRadjusted P=0.043),Moraxella (r=0.599,FDR-adjusted P=0.033),and Gastranaerophilales (r=0.589,FDRadjusted P=0.031) were positively related to 10betahydroxy-6beta-isobutyrylfuranoeremophilane.Lachnospiraceae_UCG-006 (r=0.542,FDR-adjusted P=0.042),Moraxella (r=0.574,FDR-adjusted P=0.038)and Gastranaerophilales (r=0.561,FDR-adjusted P=0.037)were also positively related to 12-oxo-20-dihydroxy-leukotriene B4.Pasteurellaceae (r=?0.556,FDR-adjusted P=0.036),Pseudobutyrivibrio (r=?0.588,FDR-adjusted P=0.01),Moraxella(r=?0.595,FDR-adjusted P=0.025)and Gastranaerophilales (r=?0.585,FDR-adjusted P=0.025)were negatively correlated with N-acetylcadaverine.In SM group,Enterorhabdus (r=0.527,FDR-adjusted P=0.047)were positively associated with 5-HMF.The preponderant strains in H group,including Prevotella_1(r=0.572,FDRadjusted P=0.040),Lachnospiraceae (r=0.578,FDRadjusted P=0.042),and Izimaplasmatales(r=0.545,FDRadjusted P=0.042) were positively correlated with 2-PBA,while Izimaplasmatales was negatively associated with methenamine (r=?0.577,FDR-adjusted P=0.041) and 6-methoxymellein (r=?0.499,FDR-adjusted P=0.002)(Additional file 1:Table S16).

    Discussion

    Milk SCC,CMT method and udder clinical symptoms were used to determine the udder health status in the current study.Screening out cows with H,SM and CM through comprehensive analysis improved the prediction accuracy of disease incidence and progression [24,25].Besides,unlike the work of Zhong et al.[26],we increased the difference gap of SCC levels among the 3 groups and considered udder clinical symptoms when grouping to further investigate the correlation between BM and changes in rumen microflora and metabolites.

    Consistent with expectations,the severity of BM was accompanied by a decline in lactation performance.The significant reduction in milk yield,milk fat and lactose contents was a reflection of the malfunction of the milk synthesis system in the mammary gland [1,3].However,the content of total milk protein in CM group was the highest among the 3 groups.Generally,milk yield is inversely proportional to milk protein [26].In line with the report by Rhone et al.[3],decreased milk yield resulted in increased dry matter content and the total protein rate in milk.In addition,raised SCC would increase the permeability of the blood-milk barrier that bring about serum proteins,immunoglobulins,and lactoferrin into the milk,which consequently increased the total amount of milk proteins [1].Furthermore,the inner-environment of rumen,including microbial compositions and their fermented products,such as VFAs,also played decisive roles in milk yield and compositions[15].In the current study,the content of ruminal VFAs and LA in CM group were significantly decreased,illustrating that the rumen environment was abnormal during acute BM.Besides,the reduced milk fat synthesis in CM may be due to the decreased acetate and butyrate which could affect the expression of acetyl-coenzyme A carboxylase and fatty acid synthase and thus changed the synthesis of milk fat [27].The less propionate in the rumen of CM cows limited the production of lactose,due to the restricted supply of glucose [12,28].The above results suggested that the alteration of rumen fermentation pattern in cows affected by CM was related to the loss of milk yield and poor milk quality,which may be correlated with the changes in rumen microbial structure and metabolic activity.

    Except for fermentation products,rumen microorganisms also have the ability to induce systemic inflammation [5,29].Shifts in microbial populations associated with dysbiosis in rumen cause translocation and systemic dissemination of increased concentrations of potentially toxic and inflammatory substances,such as lipopolysaccharide [29].In the present study,the main phyla of rumen microorganisms,Firmicutes,Proteobacteria and Bacteroidetes in cows with H,SM and CM were not significantly different,which suggested that the bacterial community at phyla level with crucial function in rumen was relatively stable in cows under different udder health status [26].Subsequent microbial HCA and PCoA analysis verified that only partial differences existed in rumen microbes from H and SM cows,which consisted with the results of Zhong et al.[26].However,the rumen bacteria structure of CM cows was significantly distinguished from the other two groups.In the present study,the rumen microflora imbalance during IMI was characterized by a large decrease in symbiotic bacteria and substantially increase in pathogens,which may be associated with the intestinal or oral inflammation.In the changed rumen flora during CM,Pasteurellaceae has been proved to cause infection through direct or indirect contact,which colonized in the oral cavity of livestock [30].Another changed pathogen from oral cavity was Neisseriaceae which was associated with a high risk of gastric cancer [31].Lachnospiraceae_UCG-006 has been reported with close correlation with intestinal inflammation and diabetes [32],and so is Gastranaerophilales [33].Besides,Ferreira et al.[34] also found that Gastranaerophilales was significantly enriched in milk from CM cows.Additionally,Moraxella is a main pathogen causing respiratory infections,laryngitis and enteritis,and also had strong correlation with rumen wall damage [29].Ralstonia was reported to induce infected hemoperitoneum in human [35],while Herbaspirillum could invade human and livestock through the digestive tract,and further cause infection in immunocompromised individuals [36,37].However,in CM group,the abundance of Pseudobutyrivibrio which was a probiotic in rumen [38] was elevated.It is speculated that,when inflammation occurs,a large number of Pseudobutyrivibrio in rumen produce butyrate to enhance the function of intestinal mucosal immune barrier,further preventing bacteria and their metabolites to enter the bloodstream and inhabiting the associated inflammation,which might exert a beneficial role [38].On the other hand,Ornithinibacillus and Enterorhabdus were identified in SM group.Although partial Bacillus in gastrointestinal tract was harmless,a mass of Ornithinibacillus in water could cause gastroenteritis [39].Enterorhabdus is a common conditional pathogen that can induce intestinal infections.Meanwhile,it is also significantly enriched in milk from mastitis cows [40].Overall,the above microorganisms,which significantly increased in the cows’ rumen during CM,have been reported with certain pathogenicity in humans or other animals,however their pathogenicity in the rumen of cows still needs further investigation.

    In contrast,the high-abundant ruminal microbiota in the H and SM groups mostly belonged to nonpathogenic symbiotic bacteria and intestinal probiotics,such as Prevotella_1,Bifidobacterium,Lachnospiraceae and Mollicutes_RF39 in H group,and Ruminiclostridium_9 in SM group.Jami and Mizrahi [41] found that Mollicutes was one of the most abundant phyla in rumen,just after Bacteroides and Pachyphyta.Similarly,Lachnospiraceae in the bovine milk samples obtained from the healthy quarters,was considered as part of the healthy milk core microbiota [42].Ruminiclostridium_9 could degrade the mucus protein of intestinal epithelial cell,produce SCFAs,and maintain the stability of the intestinal environment [43],although its functions in rumen were still unknown.

    Commensal dysbiosis in mammary gland promoted early inflammation within the udder.Further,intestinal flora dysbiosis induced by fecal bacteria transplantation could also trigger the inflammatory cells dissemination into the udder [44].Thus,we attempted to preliminarily predict the potential function of differential ruminal bacteria among the three groups in the current study.The PICRUSt is a computational approach to predict the functional composition of a metagenome using marker gene data and a database of reference genomes [21].However,the accuracy of PICRUSt (NSTI value) in the three groups was not optimal in the present study.Langille et al.[21] reported that PICRUSt could only predict a limited number of taxa based on the 16S gene copy numbers.Furthermore,the NSTI could make accurate predictions for non-animal related samples,but poor for animal related samples.

    In the present study,numerous increased ruminal microorganisms in the CM group had positive association with milk SCC.These microorganisms,as discussed above,were mostly related to oral cavity and intestinal inflammation.Usually,the elevated SCC in milk was due to pathogens invaded and colonized in udder from external environment.However,increasing studies confirmed the existence of the entero-mammary pathway [26].In fact,gastrointestinal microorganisms can be transmitted to udder through lymphatic vessels and peripheral blood circulation during lactation [45].Addis et al.[45] found that the abundance of Ruminococcus in milk from mastitis cows was significantly higher than that of healthy cows.The destruction of the integrity of rumen epithelial-vascular endothelial barrier due to ruminal damage facilitates the entry of pathogens or bacterial antigens into the portal vein and thereby causes systemic inflammation in other organs,including mammary gland [29].The above studies together with our results triggered speculation on the relationships between microbial endogenous pathways and the occurrence of mastitis.

    Besides,various systemic inflammatory disease(arthritis,spondylitis,mastitis,etc.) are related with the absence of SCFAs-producing bacteria in gastrointestinal tract [5,14].In the current study,abundant Prevotella_1,Lachnospiraceae and Bifidobacterium in H group were observed which could ferment starch and glucose to produce acetate,propionate and butyrate.Meanwhile,Ruminiclostridium_9(cellobiose-degrading bacteria) [43] and Ornithinibacillus(lactic acid-producing bacteria) in the SM group can also produce SCFAs and LA[39].This may explain their positive relationships with VFAs and LA in rumen.However,they were significantly reduced in the CM group.The SCFAs in rumen provide energy substrates for lactation,meanwhile,inhibit the proliferation of pathogens in the gastrointestinal tract [46].Propionate can enhance the resistance to colonization of Salmonella by disrupting the pH stability of bacteria[46].Butyrate can limit the growth of pathogenic Enterobacteriaceae[47].

    Metabolomics data showed that,the significantly changed metabolites in the rumen were mostly associated with inflammation and antibacterial activity.Five significantly increased ruminal metabolites were identified during mastitis,including 12-oxo-20-dihydroxy-leukotriene B4,10beta-hydroxy-6beta-isobutyrylfuranoeremophilane,5-HMF,methenamine,and 6-Methoxymellein.12-oxo-20-dihydroxy-leukotriene B4 is an isomer of leukotriene B4 (LTB4),and in this study,it was significantly increased in the CM group compared to the H group,which was consistent with the elevated LTB4 found in CM cows by Boutet et al.[48].Correlation analysis showed that ruminal bacteria,including Gastranaerophilales,Lachnospiraceae_UCG-006 and Moraxella etc.,were positively correlated with LTB4 during the CM.The lipoxygenase secreted by the pathogens initiates the conversion of arachidonic acid to LTB4,which in turn triggers an inflammatory response.Furthermore,Polymorphonuclear leukocytes (PMNLs) was the main cell type in the increased SCC during mastitis,and the LTB4 was the key chemokine for PMNLs in inflamed tissue[48],which agrees with the strongly positive correlation between 12-oxo-20-dihydroxy-leukotriene B4 and milk SCC in the current study.10beta-Hydroxy-6beta-isobutyrylfuranoeremophilane is a kind of sesquiterpenoids,participating in the degradation of limonene and pinene,both of which were volatile substances with antiinflammatory and antibacterial effects [49].The significantly increased 10beta-hydroxy-6beta-isobutyrylfuranoeremophilane in the CM group compared to the H group may thus accelerate the limonene and pinene degradation and weaken their antibacterial and antiinflammatory activity.However,reports on 10betahydroxy-6beta-isobutyrylfuranoeremophilane in mastitis are quite limited.High concentrations of 5-HMF are cytotoxic by inducing DNA damage [50].Additionally,5-HMF could also induce the production of proinflammatory factors,such as TNF-α and IL-1β [51].Both 6-methoxymellein [52] and methenamine [53] were reported to inhibit the growth of tissue cells or produce harmful metabolites,but also had certain antibacterial activities.The correlation analysis showed that methenamine and 6-methoxymellein were positively associated with milk SCC and ruminal pathogens.Therefore,the increase of methenamine and 6-Methoxymellein during SM is possibly due to either the increase of harmful metabolites produced by pathogens,or a large number of symbiotic microbiota in rumen with resistance to potential pathogens through competitive rejection and production of antibacterial compounds.Besides,compared with the H group,2-PBA,which is a derivative of butyrate,was significantly down-regulated in the CM group[54].Correlation analysis showed that it was positively correlated with butyrate and butyrate-producing bacteria in the rumen and negatively associated with milk SCC.2-PBA has a protective effect on host mucosal defense during infection,which could increase the number of intestinal Lactobacillus and reduce the induction of the pro-inflammatory cytokine IL-23 in macrophage-like cells [55].The anti-infective effect of 2-PBA may explain its negative correlation with milk SCC during BM.

    This study identified the difference of rumen microbiota structure and metabolites among CM,SM and healthy cows.The dysregulation of rumen flora and the significant changes in metabolites related to inflammation were observed during BM.However,the mechanism of the relationships between mastitis and ruminal internal environment still need further investigation.In addition,the composition of the rumen flora and metabolites in dairy cows are also affected by the environment and diet composition.Future research should further investigate the variations of rumen microbes and metabolites during BM with different dietary types,animal bedding types and seasons etc.

    Conclusion

    The current study analyzed the differences in rumen microbiome structure and metabolites among the SM,CM and healthy cows.In the rumen of CM cows,the significantly increased bacteria,such as Pseudobutyrivibrio,Gastranaerophilales and Moraxella,etc.,were found accompanied by the increase of 12-oxo-20-dihydroxy-leukotriene B4 and 10beta-hydroxy-6beta-isobutyrylfuranoeremophilane.The Ruminiclostridium_9 and Enterorhabdus in the rumen of SM cows were found increased with increasing methenamine,5-HMF and 6-methoxymellein.In contrast,the SCFAs-producing bacteria and intestinal probiotics in the rumen during IMI,including Prevoterotoella_1,Mollicutes_RF39 and Bifidobacteria,etc.,were significantly reduced with a decrease of 2-PBA,indicating the decline of immunity and antiinfection ability of the dairy cows.In summary,the abundance of microflora and metabolites associated with inflammation were significantly changed in the rumen of mastitis cows.However,the mechanism of the shift in the rumen inner-environment during BM still needs further investigation.

    Supplementary Information

    The online version contains supplementary material available at https://doi.org/10.1186/s40104-020-00543-1.

    Additional file 1 Table S1TMR Ingredient and nutrient component(%of DM).Table S2Criteria to judge California mastitis test(CMT)results.Table S3Basic information and grouping of experimental cows.Table S4Rumen microbial composition and relative abundances of cows with different udder health states(phylum level).Table S5Rumen microbial composition and relative abundances of cows with different udder health states(genus level).Table S6Differentially abundant KEGG function abundance among health,subclinical and clinical mastitis groups at level 2.Table S7Differentially abundant KEGG function abundance among health,subclinical and clinical mastitis groups at level 3.Table S8Correlation analysis between the differentially abundant bacteria and parameters of rumen fermentation as well as lactation performance.Table S9Significantly different metabolites between rumen fluid samples from healthy and clinical mastitis cows,udder positive ion mode.Table S10Significantly different metabolites between rumen fluid samples from healthy and clinical mastitis cows,udder negative ion mode.Table S11Significantly different metabolites between rumen fluid samples from healthy and subclinical mastitis cows,udder positive ion mode.Table S12Significantly different metabolites between rumen fluid samples from healthy and subclinical mastitis cows,udder negative ion mode.Table S13Significantly different metabolites between rumen fluid samples from subclinical and clinical mastitis,udder positive ion mode.Table S14Significantly different metabolites between rumen fluid samples from subclinical and clinical mastitis cows,udder negative ion mode.Table S15Correlation analysis of differential metabolites and index of lactation and rumen fermentation.Table S16Correlation analysis of differential bacteria genus and differential metabolites in rumen samples.Additional file 2 Fig.S1Pan-Species curve of OTU number in the cows with different udder health status.Pan species was the sum of all the species in a sample,which was used to observe the increase in the total number of species as the number of samples increases.H,healthy;SM,subclinical mastitis;CM,clinical mastitis.Fig.S2Rarefaction curve of OUT number in the cows with different udder health status.H,healthy;SM,subclinical mastiti;CM,clinical mastitis.Fig.S3Rumen microbial community composition analysis.aAt phylum level.bAt genus level.H,healthy;SM,subclinical mastitis;CM,clinical mastitis.Fig.S4Hierarchical cluster analysis (HCA) of rumen bacteria at genus level.Each row in the figure represents a sample,each column represents a genus,and the color indicates the relative abundance of bacteria measured in the group.Red indicates the high relative abundance,and the green indicates low relative abundance.H,healthy;SM,subclinical mastitis;CM,clinical mastitis.Fig.S5Linear discriminant analysis effect size (LEfse) analysis of multilevel species differences in ruminal microbiota.aCladogram;bLEfSe Bar graph.H,healthy;SM,subclinical mastitis;CM,clinical mastitis;LDA,linear discriminant analysis.Fig.S6The total ion chromatograms (TIC) plot of quality control(QC) samples inapositive ion mode andbnegative ion mode.Fig.S7Orthogonal partial least squares discriminant analysis(OPLS-DA) (a,c,e) and response permutation testing (RPT) (b,d,f)of rumen metabolites between H,SM and CM groups in positive ion mode.H,healthy;SM,subclinical mastitis;CM,clinical mastitis.R2X and R2Y represent the interpretation rate of the built model to the X and Y matrix,R2X (cum) and R2Y (cum) represent the cumulative interpretation rate;Q2indicates the predictive power of the model.Fig.S8Orthogonal partial least squares discriminant analysis (OPLS-DA) (a,c,e)and response permutation testing (RPT) (b,d,f)plots of rumen metabolites between H,SM and CM groups in negative ion mode.H,healthy;SM,subclinical mastitis;CM,clinical mastitis.R2X and R2Y represent the interpretation rate of the built model to the X and Y matrix,R2X(cum) and R2Y(cum) represent the cumulative interpretation rate;Q2indicates the predictive power of the model.Fig.S9HMDB compound classification (Superclass level) of significantly differential metabolites betweenaCM and H groups,bSM and H groups andcCM and SM groups.HMDB,Human Metabolome Database;H,healthy;SM,subclinical mastitis;CM,clinical mastitis.

    Abbreviations

    IMI:Intramammary infection;BM:Bovine mastitis;CM:Clinical mastitis;SM:Subclinical mastitis;H:Healthy;SCC:Somatic cell counts;CMT:California mastitis test;VFAs:Volatile fatty acids;SCFAs:Short-chain fatty acids;TVFAs:Total volatile fatty acids;Th17:T helper cell 17;LC-MS:Liquid chromatography-mass spectrometry;NH3-N:Ammonia nitrogen;RUN:Ruminal fluid urea nitrogen;LA:Lactic acid;OTUs:Operational taxonomic units;PCoA:Principal co-ordinates analysis;HCA:Hierarchical cluster analysis;PICRUSt:Phylogenetic investigation of communities by reconstruction of unobserved states;QC:Quality control;RT:Retention time;HMDB:Human metabolome database;PCA:Principle component analysis;OPLS-DA:Orthogonal partial least squares discriminate analysis;VIP:Variable importance in the projection;A/P:Acetate to propionate ratio;TIC:Total ion chromatogram;2-PBA:2-Phenylbutyric acid;ROC:Receiver operator characteristic;AUC:Area under curve;LTB4:Leukotriene B4;PMNLs:Polymorphonuclear leukocytes;5-HMF:5-Hydroxymethyl-2-furancarboxaldehyde

    Acknowledgments

    The authors thank the Beijing Key Laboratory for Dairy Cow Nutrition,Beijing University of Agriculture,Beijing,China for providing the experimental equipment.We also acknowledge the members of the State Key Laboratory of Animal Nutrition (Beijing,Chain) for their assistance with sampling and analysis.

    Authors’contributions

    The experimental scheme was designed by YW and MW.FZ,DH and FX participated in the experiment process and assisted in sampling.The analysis of experimental data and the making of charts were completed by YW.YW completed the initial draft,XN and YZ completed the overall modification of the manuscript.YZ and HW improved and polished the language of the article.LJ and JY provided the necessary experimental equipment and key guidance during the experiment process,while BX was responsible for the overall content and framework of the paper.The authors read and approved the final manuscript.

    Funding

    This study was funded by the National Key R&D Program of China (Grant No.2018YFD0500703,2017YFD0701604) and Beijing Dairy Industry Innovation Team (bjcystx-ny-1).

    Availability of data and materials

    All data generated or analyzed during this study are included in this published article (and its supplementary information files).

    Ethics approval and consent to participate

    All experimental designs and protocols were approved by the Animal Ethics Committee of the Chinese Academy of Agricultural Sciences (Beijing,China)(approval number:IAS-2019-6) and were in accordance with the recommendations of the academy’s guidelines for animal research.

    Consent for publication

    Not applicable.

    Competing interests

    The authors declare that they have no competing interests.

    Author details

    1State Key Laboratory of Animal Nutrition,Institute of Animal Science,Chinese Academy of Agricultural Sciences,Beijing 100193,China.2College of Animal Science and Technology,Northwest A&F University,Yangling 712100,China.3Beijing Key Laboratory for Dairy Cow Nutrition,Beijing University of Agriculture,Beijing 102206,China.4Engineering Research Center of Feed Development,Jiangxi Province Key Laboratory of Animal Nutrition,Jiangxi Agricultural University,Nanchang 330045,China.5Langfang Academy of Agriculture and Forestry,Langfang 065000,China.

    亚洲一区二区三区欧美精品 | 婷婷色综合大香蕉| 色5月婷婷丁香| 亚洲欧美中文字幕日韩二区| 久久精品人妻少妇| 国产探花在线观看一区二区| 国模一区二区三区四区视频| freevideosex欧美| 中文乱码字字幕精品一区二区三区| 国产精品无大码| 国产中年淑女户外野战色| 又大又黄又爽视频免费| 久久久久久久亚洲中文字幕| 王馨瑶露胸无遮挡在线观看| 亚洲精品自拍成人| 久久精品熟女亚洲av麻豆精品| 国产在视频线精品| 精品午夜福利在线看| 国产午夜福利久久久久久| 国产综合懂色| 亚洲欧美日韩卡通动漫| 永久免费av网站大全| 免费黄网站久久成人精品| 在线观看美女被高潮喷水网站| 少妇 在线观看| 久久影院123| 久久久久网色| 又爽又黄a免费视频| 亚洲综合色惰| 免费看光身美女| 欧美丝袜亚洲另类| 丰满乱子伦码专区| 少妇人妻一区二区三区视频| 黄色一级大片看看| 国产乱来视频区| 日本爱情动作片www.在线观看| 男女下面进入的视频免费午夜| 成人免费观看视频高清| 国产久久久一区二区三区| 亚洲最大成人手机在线| 嫩草影院精品99| 边亲边吃奶的免费视频| 精品视频人人做人人爽| 国产人妻一区二区三区在| 99九九线精品视频在线观看视频| 夫妻性生交免费视频一级片| 老师上课跳d突然被开到最大视频| 中文精品一卡2卡3卡4更新| 男女无遮挡免费网站观看| 久久久国产一区二区| 日日啪夜夜爽| 一区二区av电影网| 又爽又黄无遮挡网站| 男人舔奶头视频| 久久精品国产亚洲网站| 精品国产露脸久久av麻豆| 国产一区二区亚洲精品在线观看| 成人美女网站在线观看视频| av一本久久久久| av专区在线播放| 亚洲美女视频黄频| 在线亚洲精品国产二区图片欧美 | 中国美白少妇内射xxxbb| 日本一本二区三区精品| 夜夜爽夜夜爽视频| 赤兔流量卡办理| 欧美97在线视频| 少妇的逼水好多| 亚洲人成网站在线观看播放| 亚洲人成网站在线观看播放| 91狼人影院| 精品少妇黑人巨大在线播放| 欧美变态另类bdsm刘玥| 国产高清三级在线| 免费看日本二区| 大香蕉97超碰在线| 欧美另类一区| 九九久久精品国产亚洲av麻豆| 日本色播在线视频| 黄片无遮挡物在线观看| 亚洲第一区二区三区不卡| 亚洲欧美一区二区三区黑人 | 成人免费观看视频高清| 99re6热这里在线精品视频| 波多野结衣巨乳人妻| 中文字幕人妻熟人妻熟丝袜美| 国产69精品久久久久777片| 久久精品国产亚洲网站| 大香蕉97超碰在线| 下体分泌物呈黄色| 国产淫片久久久久久久久| 久久久久九九精品影院| 国内精品宾馆在线| 国内精品宾馆在线| 亚洲色图综合在线观看| 亚洲欧美日韩无卡精品| 人妻 亚洲 视频| 国精品久久久久久国模美| 日本免费在线观看一区| 直男gayav资源| 日韩欧美精品免费久久| 在线观看一区二区三区| 国产有黄有色有爽视频| 精品人妻一区二区三区麻豆| 男男h啪啪无遮挡| 国产精品久久久久久精品电影| 三级男女做爰猛烈吃奶摸视频| av在线app专区| 国产精品女同一区二区软件| 午夜亚洲福利在线播放| 欧美性猛交╳xxx乱大交人| 久久久久久国产a免费观看| 美女国产视频在线观看| 亚洲精品色激情综合| 国产精品精品国产色婷婷| 欧美zozozo另类| av一本久久久久| 嫩草影院新地址| 麻豆成人午夜福利视频| 97精品久久久久久久久久精品| 国产69精品久久久久777片| 蜜桃亚洲精品一区二区三区| 美女视频免费永久观看网站| 建设人人有责人人尽责人人享有的 | 久久鲁丝午夜福利片| 午夜爱爱视频在线播放| 国产精品久久久久久精品电影小说 | 国产真实伦视频高清在线观看| 免费看a级黄色片| 街头女战士在线观看网站| 欧美日韩视频精品一区| 精品一区二区免费观看| 各种免费的搞黄视频| 91久久精品国产一区二区三区| 91在线精品国自产拍蜜月| 亚洲欧美日韩卡通动漫| 国产免费视频播放在线视频| 岛国毛片在线播放| 五月天丁香电影| 精品人妻熟女av久视频| 狠狠精品人妻久久久久久综合| 涩涩av久久男人的天堂| 99久久精品一区二区三区| 免费大片黄手机在线观看| 精品国产露脸久久av麻豆| 欧美精品人与动牲交sv欧美| 欧美日韩在线观看h| 欧美日韩精品成人综合77777| 国产黄a三级三级三级人| 黄色怎么调成土黄色| 热99国产精品久久久久久7| 免费在线观看成人毛片| 亚洲精品乱码久久久v下载方式| 国产精品人妻久久久久久| 九九在线视频观看精品| 亚洲精品日韩在线中文字幕| 深夜a级毛片| 日韩av免费高清视频| 天美传媒精品一区二区| 精品久久久精品久久久| 成年免费大片在线观看| av黄色大香蕉| 亚洲丝袜综合中文字幕| 内射极品少妇av片p| 欧美三级亚洲精品| 中文欧美无线码| 日韩伦理黄色片| eeuss影院久久| 国产精品av视频在线免费观看| 日韩精品有码人妻一区| 国产伦精品一区二区三区视频9| 国产精品一及| 噜噜噜噜噜久久久久久91| 亚洲精品影视一区二区三区av| av女优亚洲男人天堂| 亚洲欧美清纯卡通| 国产中年淑女户外野战色| 亚洲精品456在线播放app| 免费看av在线观看网站| 精品亚洲乱码少妇综合久久| 男女无遮挡免费网站观看| 免费看日本二区| 99热这里只有是精品在线观看| 久久久久网色| 熟女人妻精品中文字幕| 在线亚洲精品国产二区图片欧美 | 国产一区二区在线观看日韩| 免费播放大片免费观看视频在线观看| 亚洲国产欧美在线一区| 日韩av免费高清视频| 嫩草影院新地址| 69人妻影院| 一区二区三区乱码不卡18| 99热全是精品| 国产色婷婷99| kizo精华| 国产av码专区亚洲av| av播播在线观看一区| 成人黄色视频免费在线看| 亚洲怡红院男人天堂| 精品久久久久久久久av| 激情 狠狠 欧美| 精品久久久久久久末码| 97热精品久久久久久| 午夜激情福利司机影院| 青春草国产在线视频| 黄色欧美视频在线观看| 十八禁网站网址无遮挡 | 爱豆传媒免费全集在线观看| 在线观看三级黄色| 91久久精品国产一区二区成人| 大香蕉97超碰在线| 国产视频内射| 精品视频人人做人人爽| 免费看不卡的av| 男人添女人高潮全过程视频| 久久精品国产鲁丝片午夜精品| 亚洲人与动物交配视频| 晚上一个人看的免费电影| 婷婷色av中文字幕| 51国产日韩欧美| 99九九线精品视频在线观看视频| 啦啦啦中文免费视频观看日本| 欧美激情在线99| 欧美少妇被猛烈插入视频| 欧美性感艳星| 国产视频首页在线观看| 国产免费视频播放在线视频| 日日摸夜夜添夜夜爱| 在线免费十八禁| 久久久久久久久久久丰满| 亚洲丝袜综合中文字幕| 肉色欧美久久久久久久蜜桃 | av一本久久久久| 国产精品熟女久久久久浪| 国产视频内射| 亚洲人成网站在线播| 99久久精品热视频| 男女下面进入的视频免费午夜| 免费电影在线观看免费观看| 国产亚洲午夜精品一区二区久久 | 熟女人妻精品中文字幕| videos熟女内射| 在线播放无遮挡| 纵有疾风起免费观看全集完整版| 只有这里有精品99| 日韩av免费高清视频| 欧美成人a在线观看| av女优亚洲男人天堂| 国产成人精品一,二区| 精品少妇黑人巨大在线播放| 亚洲成人av在线免费| 久久精品久久久久久噜噜老黄| 久久久久久国产a免费观看| 男人爽女人下面视频在线观看| 精品亚洲乱码少妇综合久久| 国产毛片a区久久久久| 日韩一本色道免费dvd| 亚洲伊人久久精品综合| 久久精品国产亚洲av涩爱| videos熟女内射| 久久亚洲国产成人精品v| 欧美三级亚洲精品| 精品久久久久久久久亚洲| av在线天堂中文字幕| 在线免费观看不下载黄p国产| 老女人水多毛片| 美女脱内裤让男人舔精品视频| 日韩伦理黄色片| 一二三四中文在线观看免费高清| 国产亚洲一区二区精品| 小蜜桃在线观看免费完整版高清| 身体一侧抽搐| 亚洲av男天堂| 最近手机中文字幕大全| 99久久人妻综合| av天堂中文字幕网| 黄片wwwwww| 一二三四中文在线观看免费高清| 男人狂女人下面高潮的视频| 色5月婷婷丁香| av国产免费在线观看| 亚洲av一区综合| 亚洲精品一区蜜桃| 天天一区二区日本电影三级| 国产精品成人在线| 青春草亚洲视频在线观看| 青春草国产在线视频| 亚洲欧美日韩东京热| 欧美xxⅹ黑人| 永久免费av网站大全| 一区二区三区精品91| 寂寞人妻少妇视频99o| 在线观看av片永久免费下载| 免费看a级黄色片| 免费看a级黄色片| 国产一区二区三区av在线| 久久久久精品久久久久真实原创| 亚洲一级一片aⅴ在线观看| 一级毛片电影观看| 男女边吃奶边做爰视频| 国产精品国产三级国产av玫瑰| 激情五月婷婷亚洲| 亚洲,一卡二卡三卡| 插阴视频在线观看视频| 黄色视频在线播放观看不卡| 亚洲欧美日韩卡通动漫| 国产黄片视频在线免费观看| 亚洲欧洲日产国产| 精品一区在线观看国产| 一级毛片电影观看| 亚洲aⅴ乱码一区二区在线播放| 午夜福利视频精品| 插阴视频在线观看视频| 免费观看无遮挡的男女| 国产一级毛片在线| 国产探花极品一区二区| 久久久亚洲精品成人影院| 久久99热这里只有精品18| 丰满人妻一区二区三区视频av| 久久国内精品自在自线图片| 又爽又黄a免费视频| 你懂的网址亚洲精品在线观看| 精品人妻视频免费看| 18禁在线播放成人免费| 午夜激情福利司机影院| 国产免费又黄又爽又色| 国产精品国产三级专区第一集| 午夜亚洲福利在线播放| 免费看日本二区| 丝袜喷水一区| 岛国毛片在线播放| 亚洲国产av新网站| 高清欧美精品videossex| 国产大屁股一区二区在线视频| 夜夜看夜夜爽夜夜摸| 九九久久精品国产亚洲av麻豆| 日韩欧美精品免费久久| 国产伦精品一区二区三区四那| 自拍偷自拍亚洲精品老妇| 亚洲av一区综合| 如何舔出高潮| 国产欧美亚洲国产| 26uuu在线亚洲综合色| 男插女下体视频免费在线播放| 大又大粗又爽又黄少妇毛片口| 2022亚洲国产成人精品| 久久久久久国产a免费观看| 热re99久久精品国产66热6| 久久女婷五月综合色啪小说 | 国模一区二区三区四区视频| av免费观看日本| 大香蕉97超碰在线| 人妻系列 视频| 岛国毛片在线播放| 最近最新中文字幕免费大全7| 大香蕉97超碰在线| 国产成人午夜福利电影在线观看| 国产真实伦视频高清在线观看| 欧美xxxx黑人xx丫x性爽| 黄色欧美视频在线观看| 亚洲真实伦在线观看| 女人十人毛片免费观看3o分钟| 欧美高清性xxxxhd video| 在线天堂最新版资源| 久久久成人免费电影| 一本一本综合久久| 尾随美女入室| 少妇丰满av| 在现免费观看毛片| 内射极品少妇av片p| 亚洲一级一片aⅴ在线观看| 99热6这里只有精品| 欧美最新免费一区二区三区| 日韩制服骚丝袜av| 男人爽女人下面视频在线观看| 国产精品99久久久久久久久| 亚洲精品乱码久久久久久按摩| 欧美性感艳星| 女的被弄到高潮叫床怎么办| 欧美日韩一区二区视频在线观看视频在线 | 国产片特级美女逼逼视频| 99久久人妻综合| 亚洲人与动物交配视频| 99热这里只有精品一区| 亚洲最大成人中文| 七月丁香在线播放| 午夜激情久久久久久久| 黄色视频在线播放观看不卡| 日日摸夜夜添夜夜添av毛片| 精品国产乱码久久久久久小说| 日韩一区二区三区影片| 国产有黄有色有爽视频| 听说在线观看完整版免费高清| av在线老鸭窝| 国产精品成人在线| 亚洲四区av| 夫妻性生交免费视频一级片| 亚洲av福利一区| 99热这里只有是精品在线观看| 神马国产精品三级电影在线观看| 久久韩国三级中文字幕| 久久久精品免费免费高清| 97超视频在线观看视频| 国产 精品1| 久久鲁丝午夜福利片| 欧美bdsm另类| 成年人午夜在线观看视频| 亚洲av二区三区四区| 黄色怎么调成土黄色| 黄色欧美视频在线观看| 汤姆久久久久久久影院中文字幕| av国产精品久久久久影院| 黄色配什么色好看| 爱豆传媒免费全集在线观看| 超碰97精品在线观看| 中文欧美无线码| 久久6这里有精品| 国产日韩欧美在线精品| 久久久久久伊人网av| 日本一二三区视频观看| 大陆偷拍与自拍| 成人国产av品久久久| 好男人视频免费观看在线| 国产综合懂色| 欧美人与善性xxx| 精品久久久久久久人妻蜜臀av| 男女国产视频网站| 校园人妻丝袜中文字幕| 国产黄色视频一区二区在线观看| 精品99又大又爽又粗少妇毛片| 国产在视频线精品| 久久精品综合一区二区三区| 一区二区av电影网| 精品午夜福利在线看| 老司机影院成人| 国产91av在线免费观看| 亚洲三级黄色毛片| 2022亚洲国产成人精品| 国产男人的电影天堂91| 日韩 亚洲 欧美在线| 80岁老熟妇乱子伦牲交| 久久热精品热| 成人特级av手机在线观看| 一级片'在线观看视频| 69人妻影院| videos熟女内射| 看免费成人av毛片| 噜噜噜噜噜久久久久久91| av卡一久久| 国产伦在线观看视频一区| 亚洲一级一片aⅴ在线观看| 亚洲成人一二三区av| 日本色播在线视频| 一本久久精品| 国产毛片a区久久久久| 亚洲婷婷狠狠爱综合网| 国产毛片a区久久久久| 18禁裸乳无遮挡免费网站照片| 亚洲成人av在线免费| 少妇的逼好多水| 又黄又爽又刺激的免费视频.| 欧美日韩精品成人综合77777| 欧美成人午夜免费资源| 国产 一区精品| 看黄色毛片网站| 国产真实伦视频高清在线观看| 欧美激情在线99| 在线a可以看的网站| 久久97久久精品| 精品久久久久久久人妻蜜臀av| 日本猛色少妇xxxxx猛交久久| 国产91av在线免费观看| 18禁裸乳无遮挡动漫免费视频 | 免费大片黄手机在线观看| 国产一区亚洲一区在线观看| 人妻 亚洲 视频| 日韩大片免费观看网站| 看黄色毛片网站| 久久久久久久午夜电影| 黄色欧美视频在线观看| 舔av片在线| 在线观看一区二区三区激情| 色哟哟·www| 国产色爽女视频免费观看| av播播在线观看一区| 狠狠精品人妻久久久久久综合| 深夜a级毛片| 久久精品国产a三级三级三级| 精品国产三级普通话版| 亚洲人成网站高清观看| 国产欧美另类精品又又久久亚洲欧美| 中文字幕人妻熟人妻熟丝袜美| 一级片'在线观看视频| 黄色欧美视频在线观看| 偷拍熟女少妇极品色| 亚洲人成网站在线观看播放| 国产成人91sexporn| 99久久人妻综合| 街头女战士在线观看网站| 蜜桃久久精品国产亚洲av| eeuss影院久久| 国产成人a区在线观看| 精品久久久久久久久亚洲| 在线免费十八禁| 免费观看av网站的网址| 大片免费播放器 马上看| 人妻制服诱惑在线中文字幕| 一本久久精品| 国产成人精品久久久久久| 亚洲欧美清纯卡通| 亚洲精品亚洲一区二区| 国产精品99久久99久久久不卡 | 高清视频免费观看一区二区| 卡戴珊不雅视频在线播放| 久久精品国产亚洲av涩爱| 精品少妇久久久久久888优播| 男人爽女人下面视频在线观看| 久久久a久久爽久久v久久| 51国产日韩欧美| 亚洲人成网站在线播| 成人综合一区亚洲| 免费看光身美女| 日韩人妻高清精品专区| 欧美精品一区二区大全| kizo精华| 欧美日本视频| 日韩视频在线欧美| 亚洲欧美一区二区三区黑人 | 午夜福利网站1000一区二区三区| 日日啪夜夜撸| 亚洲精品aⅴ在线观看| 成人黄色视频免费在线看| 国产av国产精品国产| 亚洲丝袜综合中文字幕| 男人舔奶头视频| 欧美性感艳星| 99热这里只有精品一区| 在线a可以看的网站| 菩萨蛮人人尽说江南好唐韦庄| 久久久久精品性色| 欧美潮喷喷水| 一区二区三区精品91| 在线观看一区二区三区| 免费看光身美女| 夫妻性生交免费视频一级片| 国产av国产精品国产| 99久久精品一区二区三区| 午夜福利高清视频| 日本av手机在线免费观看| 国产精品国产三级国产专区5o| 国产成人午夜福利电影在线观看| 一级毛片aaaaaa免费看小| 麻豆久久精品国产亚洲av| 97在线视频观看| 免费黄频网站在线观看国产| 日本免费在线观看一区| 欧美bdsm另类| 国产人妻一区二区三区在| 国产美女午夜福利| 国产成人一区二区在线| 久久99蜜桃精品久久| 午夜爱爱视频在线播放| 亚洲精品国产av蜜桃| 亚洲最大成人手机在线| 尾随美女入室| 欧美日韩国产mv在线观看视频 | 黄色一级大片看看| 只有这里有精品99| 高清日韩中文字幕在线| 婷婷色综合www| 狂野欧美激情性bbbbbb| 欧美人与善性xxx| 菩萨蛮人人尽说江南好唐韦庄| 国产精品熟女久久久久浪| 午夜免费男女啪啪视频观看| 日韩av在线免费看完整版不卡| 日本色播在线视频| 成人无遮挡网站| 十八禁网站网址无遮挡 | 精品少妇黑人巨大在线播放| 国产精品无大码| 青春草亚洲视频在线观看| 男人爽女人下面视频在线观看| 亚洲av成人精品一区久久| 成人毛片60女人毛片免费| 久久精品人妻少妇| 亚洲成人中文字幕在线播放| 亚洲av国产av综合av卡| 黄片无遮挡物在线观看| 大片电影免费在线观看免费| 99热6这里只有精品| 成人免费观看视频高清| 久久久久久伊人网av| 亚洲国产色片| 亚洲国产精品999| 少妇的逼水好多| 国产免费一级a男人的天堂| 久久综合国产亚洲精品| 精品久久久久久久人妻蜜臀av| 国产黄色免费在线视频| 777米奇影视久久| 精品一区二区三区视频在线| 亚洲国产精品999| 亚洲av免费在线观看| 亚洲欧美日韩另类电影网站 | 亚洲精品日本国产第一区| 99久久人妻综合| 纵有疾风起免费观看全集完整版| 搡老乐熟女国产| 蜜桃亚洲精品一区二区三区| 精品国产露脸久久av麻豆| 久久影院123| 久久这里有精品视频免费| 超碰av人人做人人爽久久|