s2fieldsofstudy:
Biology
Agricultural And Food Sciences

text:
Genomic and phenotypic biology of a novel Dickeya zeae WH1 isolated from rice in China: Insights into pathogenicity and virulence factors

Soft rot caused by Dickeya zeae is an important bacterial disease affecting rice and other plants worldwide. In this study, Nanopore and Illumina sequencing platforms were used to sequence the high-quality complete genome of a novel D. zeae strain WH1 (size: 4.68 Mb; depth: 322.37x for Nanopore, 243.51x for Illumina; GC content: 53.59%), which was isolated from healthy rice root surface together with Paenibacillus polymyxa, a potential biocontrol bacterium against D. zeae strain WH1. However, the pure WH1 culture presented severe pathogenicity. Multilocus sequence analysis (MLSA) indicated that strains WH1, EC1, and EC2 isolated from rice were grouped into a clade differentiated from other D. zeae strains. The average nucleotide identity (ANI) and DNA-DNA hybridization (DDH) analyses demonstrated that WH1 was phylogenetically closest to EC2. Furthermore, the pathogenicity determinants and virulence factors of WH1 were mainly analyzed through genomic comparison with complete genomes of other D. zeae strains with high virulence (EC1, EC2, MS1, and MS2). The results revealed that plant cell wall-degrading extracellular enzymes (PCWDEs), flagellar and chemotaxis, and quorum sensing were highly conserved in all analyzed genomes, which were confirmed through phenotypic assays. Besides, WH1 harbored type I, II, III, and VI secretion systems (T1SS, T2SS, T3SS, and T6SS), but lost T4SS and T5SS. Like strains MS1 and MS2 isolated from bananas, WH1 harbored genes encoding both capsule polysaccharide (CPS) and exopolysaccharide (EPS) biosynthesis. The results of pathogenicity assays demonstrated that WH1 produced severe soft rot symptoms on potato tubers, carrots, radishes, and Chinese cabbage. Meanwhile, WH1 also produced phytotoxin(s) to inhibit rice seed germination with an 87% inhibitory rate in laboratory conditions. More importantly, we confirmed that phytotoxin(s) produced by WH1 are different from zeamines produced by EC1. Comparative genomics analyses and phenotypic and pathogenicity assays suggested that WH1 likely evolved through a pathway different from the other D. zeae strains from rice, producing a new type of rice foot rot pathogen. These findings highlight the emergence of a new type of D. zeae strain with high virulence, causing soft rot in rice and other plants.

Soft rot caused by Dickeya zeae is an important bacterial disease affecting rice and other plants worldwide. In this study, Nanopore and Illumina sequencing platforms were used to sequence the high-quality complete genome of a novel D. zeae strain WH1 (size: 4.68 Mb; depth: 322.37x for Nanopore, 243.51x for Illumina; GC content: 53.59%), which was isolated from healthy rice root surface together with Paenibacillus polymyxa, a potential biocontrol bacterium against D. zeae strain WH1. However, the pure WH1 culture presented severe pathogenicity. Multilocus sequence analysis (MLSA) indicated that strains WH1, EC1, and EC2 isolated from rice were grouped into a clade differentiated from other D. zeae strains. The average nucleotide identity (ANI) and DNA-DNA hybridization (DDH) analyses demonstrated that WH1 was phylogenetically closest to EC2. Furthermore, the pathogenicity determinants and virulence factors of WH1 were mainly analyzed through genomic comparison with complete genomes of other D. zeae strains with high virulence (EC1, EC2, MS1, and MS2). The results revealed that plant cell wall-degrading extracellular enzymes (PCWDEs), flagellar and chemotaxis, and quorum sensing were highly conserved in all analyzed genomes, which were confirmed through phenotypic assays. Besides, WH1 harbored type I, II, III, and VI secretion systems (T1SS, T2SS, T3SS, and T6SS), but lost T4SS and T5SS. Like strains MS1 and MS2 isolated from bananas, WH1 harbored genes encoding both capsule polysaccharide (CPS) and exopolysaccharide (EPS) biosynthesis. The results of pathogenicity assays demonstrated that WH1 produced severe soft rot symptoms on potato tubers, carrots, radishes, and Chinese cabbage. Meanwhile, WH1 also produced phytotoxin(s) to inhibit rice seed germination with an 87% inhibitory rate in laboratory conditions. More importantly, we confirmed that phytotoxin(s) produced by WH1 are different from zeamines produced by EC1. Comparative genomics analyses and phenotypic and pathogenicity assays suggested that WH1 likely evolved Introduction Dickeya spp. are Gram-negative and causative agents of a wide range of bacterial soft rot and have been listed among the top 10 most important bacterial phytopathogens due to the high costs and economic consequences they incur (Mansfield et al., 2012). Currently, the Dickeya genus encompasses 12 recognized species (Boluk et al., 2021). Of which, D. zeae strains were isolated from a wide and diverse range of hosts, causing serious economic losses in crop yield in different parts of the world, especially for rice, potato, pineapple, and banana (Hu et al., 2018;Feng et al., 2019;Boluk et al., 2021). Infected rice plants by D. zeae have a dark brown decay of the tillers at the site of infection with a bad smell, which later leads to the collapse of the entire plant (Hussain et al., 2008). In addition, D. zeae strains can reside on the surface and/or within the intercellular spaces of plant organs without exhibiting symptoms (Pérombelon, 1992). The reason these bacteria can grow on some plants without symptoms but cause devastating diseases in others remains unclear (De Boer, 2002;Van Gijsegem et al., 2021). However, when in a susceptible host, these strains may persist for months until the environmental conditions are not favorable to disease initiation (De Boer, 2002;Kraepiel and Barny, 2016). Additionally, they can persist overwinter in contaminated plant residue and soil. Under certain conditions, such as decreasing O 2 concentration, high temperature, high humidity, and there being water film on the surface of the plant organs, latency is broken and the bacteria start to grow and produce plant cell wall degrading enzymes (PCWDEs) leading to disease development in susceptible hosts (Bez et al., 2021).
Generally, D. zeae strains infect plants through three main steps: (1) adhesion to plant surface and penetration in the plant tissues, either via wound sites or through natural openings such as stomata; (2) invasion of apoplasts; and (3) plant cell wall degradation (Reverchon and Nasser, 2013). Some typical virulence factors are involved in these processes, such as type I-VI secretion systems (T1SS-T6SS), PCWDEs, flagella, and flagella-mediated motility, quorum sensing signal molecules acylhomoserine lactone (AHL) and Vfm, and biofilms . A previous study showed that the zeamines produced by D. zeae EC1 were important virulence factors for rice, but antagonists for many bacteria and fungi (Liao et al., 2014). However, the zms gene cluster encoding zeamines biosynthesis was absent in the genome of D. zeae strain EC2 . Additionally, other virulence determinants are present among D. zeae strains, such as T3SS, which is present in strain EC1 but absent in strain EC2 (Zhou et al., 2015;Zhang et al., 2020). Therefore, the diversity involved in virulence factors would be present in the genomes of D. zeae strains, even though they were isolated from the same host.
In this study we obtained a complete genome sequence of D. zeae strain WH1, isolated from the rice root surface, using Illumina next-generation sequencing and Nanopore sequencing technologies. The genome sequence was annotated and compared with the representative genome sequences of other D. zeae strains (EC1, EC2, MS1, and MS2), focusing on virulence determinants and potential regulatory mechanisms. Substantial phenotype and pathogenicity assays were performed for strain WH1. The genome-wide comparison and phenotypic and pathogenicity assays could help to clarify the differentiation among D. zeae strains from rice and their genetic adaptions to the host. This information directs attention to this newly emerging rice pathogen to provide better insights into the early prevention of rice foot rot disease.

Strain WH1 isolation and characterization
Three healthy rice roots were collected from the field of Nanling, Wuhu City, Anhui Province, China, in August 2021. The rice roots with few adherent soils were cut into pieces and placed into 10 mL of sterile water. The mixtures were incubated at 25 • C for 30 min with shaking at 150 rpm. The mixture at different dilutions with sterile water was spread onto LB agar plates and incubated at 25 • C for 48 h. Single colonies were picked and sub-cultured three times. Pure cultures were used to screen antagonists against rice pathogenic fungi. One strain, WH1, had the strong ability to inhibit the growth of pathogenic fungi (data not shown). The Frontiers in Microbiology 02 frontiersin.org strain was identified as Dickeya genus through sequencing of the 16S rRNA gene using 27F/1492R primers. Subsequently, the strain was further confirmed through Koch's postulates. Additionally, the screening of antagonistic bacteria against WH1 was performed according to the methods previously described (Li et al., 2020). The screened antagonistic bacteria were identified through sequencing of the 16S rRNA gene. All isolated bacteria were stored in glycerol at -80 • C at Anhui Normal University, China. The strain WH1 was mainly investigated in this study. For subculture, the strain was grown at 28 • C overnight in LB broth (10 g/L tryptone, 5 g/L yeast extract, and 10 g/L NaCl, pH7.0) in an incubator shaker at 150 rpm.

Genomic DNA extraction and genome sequencing
The strain WH1 was incubated to logarithmic phase at 28 • C and its genomic DNA was extracted using the TIANamp Bacterial DNA kit (Tiangen Biotech, Beijing, China) according to the manufacturer's instructions. The quality of the obtained genomic DNA was detected with agarose gel electrophoresis, NanoDrop One Spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA), and Qubit 3.0 Fluorimeter (Life Technologies, Carlsbad, CA, USA), respectively. The DNA libraries for MinION Nanopore sequencing and Illumina sequencing were constructed, respectively. DNA sequencing was respectively carried out on the PromethION sequencer (Oxford Nanopore Technologies, Oxford, United Kingdom) and NovaSeq 6000 (Illumina, USA) platforms at BENAGEN Technologies, Wuhan, China.

Genome assembly and annotation
Unicycler (version 0.5.0) software was used to assemble short reads from the Illumina platform and long reads from the MinION Nanopore platform (Wick et al., 2017). Gene prediction was performed using Prokka (version 1.14.6) software (Seemann, 2014); Aragorn was used to predict the tRNA of the genome (Laslett and Canback, 2004); rRNAmmer was used to predict genomic rRNA sequences (Lagesen et al., 2007). Gene function annotation was performed using the local BLAST method to compare the predicted gene sequences with GO (Ashburner et al., 2000), KEGG (Kanehisa et al., 2004), COG (Tatusov et al., 2000), UniProt (Apweiler et al., 2004), and NR databases, respectively, to obtain the corresponding annotation information. The genomic sequence was submitted to the NCBI GenBank Genome database under the accession number CP101400.

Phylogenetic analysis of Dickeya strains
Sixteen Dickeya strains with complete genome sequences were used to extract the gene sequences of aptD, gyrB, inf B, and rpoB for multilocus sequence analysis (MLSA). The concatenated sequences of these four genes were processed on MEGA5.2 using a neighbor-joining algorithm with 1,000 bootstrapped replications (Tamura et al., 2011). To obtain an estimate of overall genomic similarities, a pairwise comparison of the WH1 genome with 15 genomes of other Dickeya strains was performed using the average nucleotide identity (ANI) based on the nucleotide MUMmmer algorithm (ANIm) in the Jspecies Web server (Richter et al., 2016). In addition, in silico DNA-DNA hybridization (isDDH) values were calculated using the Genome-to-Genome Distance Calculator (GGDC) with formula 2 and BLAST + alignment. 1 The cut-off values of 96 and 70% were assigned as a species-delineation framework for ANI and isDDH, respectively (Goris et al., 2007).

Genomic comparison between WH1 and other Dickeya zeae strains
The complete genome of WH1 was compared with the other four complete genomes of D. zeae strains (EC1, EC2, MS1, and MS2). Their genome sequences were retrieved from the NCBI genome database. The basic genomic profile features of these four strains were taken from the NCBI GenBank database and the bioinformatics resource center proteome comparison tool of the pathosystems resource integration center (PATRIC) Web Server (Davis et al., 2020). To identify the lists of shared and unique gene clusters, Pan-genome analysis was performed using Pan-Genomes Analysis Web Server . The orthologous gene clusters were identified by the Gene Family (GF) method set to its default parameters.
The soft rot bacteria within Dickeya genus macerate the plant tissue and acquire nutrients from the dead cells. Most pathogenicity determinants related encoding genes, including PCWDEs, type secretion systems (types I-VI), extracellular polysaccharides (such as CPS and EPS), bacterial attachment operons (type IV pili), flagella and chemotaxis, zeamine synthesis, quorum sensing systems, and biofilm formation, have been described in strains EC1, EC2, MS1, or MS2 (Feng et al., 2019;Zhang et al., 2020Zhang et al., , 2022Liu et al., 2022). In this study, we mainly analyzed and compared the similarities, differences, or absence of gene encoding virulence determinants between WH1 and these four D. zeae strains.

Extracellular enzyme activity assays
Pectate lyases (Pel), cellulases (Cel), and proteases (Prt) produced by Dickeya genus were the major causes resulting in foot rot in rice and other plants (Boluk et al., 2021). Therefore, their activities were detected according to the protocols as previously described with a few modifications (Chatterjee et al., 1995;Hu et al., 2022). In brief, the media for enzyme activity assays as follows: (1) Pel assay medium (per liter): 10 g polygalacturonic acid, 10 g yeast extract, 0.38 µmol CaCl 2 , 100 mmol Tris-HCl, 15 g agar, pH8.5; (2) Cel assay medium (per liter): 1.0 g carboxymethyl cellulose, 3.8 g Na 3 PO 4 , 15 g agar, pH7.0; (3) Prt assay medium: LB agar medium containing 1% (v/v) skimmed milk. The media were poured and allowed to solidify and 5 mm diameter wells were punched into the agar and sealed the bottom with worm 1.5% agar solution. 25 µL overnight cultures were applied to each well, and plates were incubated at 28 • C. After 15 h, Pel assay plates were flooded with 1 M HCl; Cel assay plates were flooded with 1% Congo red solution for 1 h and washed twice with 1 M NaCl solution. Next, the treated Pel and Cel assay plates were incubated at room temperature for 24 h and the diameter of haloes around the wells were measured. For Prt assay, the haloes around the wells were measured after 24 h without any treatment. Each treatment was carried out at least three times, and all assays were repeated at least three times.

Flagellum-mediated motility assays
These assays were performed according to the methods as previously described (Boluk et al., 2021). Media for swimming (per liter contains 10 g tryptone, 5 g NaCl, and 3 g agar) and swarming (per liter contains 5 g peptone, 3 g yeast extract, and 4 g agar) assays were prepared. One microliter of overnight bacterial culture was spotted onto the center of a plate (60 mm diameter) containing about 10 mL of each medium. The plates were incubated at 28 • C for 20 h before measurement of the diameters of each bacterial motility zone. Each experiment was repeated at least three times in triplicate.

Assays of exopolysaccharide production
Exopolysaccharide (EPS) was detected in WH1 culture according to a previous study with minor modifications (Hu et al., 2018). In brief, a single colony of WH1 strain was picked and transferred into a 2 mL LB medium and incubated at 28 • C overnight. Afterward, 1 mL of which was applied into 100 mL LB medium and grown with shaking at 150 rpm for 12 h. Cultures were centrifuged at 8,000 rpm for 10 min, and then at 4,000 rpm for 20 min to obtain 80 mL supernatants.
Double volumes of ethanol were added to the supernatants, mixed thoroughly, stored at 4 • C overnight for precipitation, and subjected to centrifuge at 8,000 rpm for 20 min. Finally, supernatants were discarded and pellets were weighed after drying at 55 • C overnight. The experiment was repeated three times in triplicate.

Biofilm formation assays
An overnight bacterial culture was diluted at 1:100 with LB medium; 200 µL of the culture was added into each well of 96-well microtiter plates and statically incubated at 28 • C for 20 h. For qualification of the total biofilm mass, the suspension culture was removed and the plate was washed three times with PBS. After drying for 2 h at 60 • C, the biofilm was stained with 50 µL of 0.1% crystal violet (w/v) for 15 min. The wells were also washed three times with PBS to remove unbound crystal violet dye and dried for 2 h at 60 • C. Biofilms formed in 96well plate were observed. Additionally, the biofilms developed in a 48-well microtiter plate containing one glass coverslip per well were used during observation with a fluorescent microscope, as previously described (Merritt et al., 2005). The plate was statically incubated at 28 • C for 20 h. The AO/EB staining kit (Sango Biotech, Shanghai, China) was used to stain biofilms developed on glass coverslip after being washed three times with PBS to discard planktonic cells. The images of biofilms were observed through the fluorescent microscope (Leica Microsystem, Germany).

Acylated homoserine lactones production assay
Chromobacterium violaceum CV026 is usually used as a biosensor to visualize acylated homoserine lactones (AHLs) with N-acyl side chains from C4 to C8 in length production by Gramnegative bacteria . In this study, one colony of CV026 was inoculated into 2 mL of LB medium and incubated at 28 • C overnight. The CV026 cultures were added to a warm LB agar (1.5%) medium at a ratio of 1:100, and then the mixtures were poured immediately over the surface of LB agar plates prepared in Petri dishes. When the overlaid agar solidified, a well (5 mm diameter) was made in the center of each plate and sealed the bottom with warm agar solution, 25 µL of overnight WH1 culture (10 9 CFU/mL) was added to the well. Meanwhile, 25 µL of overnight Pseudomonas aeruginosa PAO1 culture (10 9 CFU/mL) was performed under the same conditions as the reference. Violacein halo production was observed after being incubated at 28 • C for 24 h. The diameter of the violacein halo was also measured.
AHL production was determined by inoculating 200 µL of an overnight culture of WH1 into 200 mL of LB. After 20 h incubation at 28 • C, the sterile supernatant was polled and extracted with acidified ethyl acetate according to a protocol described previously (Zhou et al., 2017). The AHL produced by WH1 was analyzed using an HPLC system (Shimadzu, Japan) equipped with a C18 column by UV absorbance at 210 nm. The mobile phase A was water, and the mobile phase B was methanol. The following rate was set as 0.8 mL/min. The injection volume was 20 µL. The peak corresponding to 3OC6-HSL was identified according to the retention time of commercial 3OC6-HSL standard (Aladdin, Shanghai, China) with the same HPLC protocol.

Pathogenicity assays
Thirty rice seeds were put into 10 mL of sterile water containing 1 mL of WH1 overnight culture and incubated at room temperature for 5 h. The rice seeds were rinsed three times with sterile water and transferred onto the top of two moistened filter papers placed on a sterilized plate at room temperature. Rice seeds were treated with the same volume of sterile water as a negative control. The rice seeds were then incubated at 33 • C under 8 h dark/16 h light with supplementation of sterile water when necessary. The seed germination rate was measured after 7 days.
To investigate the pathogenicity of strain WH1 on other plants, we selected potato tubers, carrots, radishes, and Chinese cabbage. The plants were washed with deionized water and then surface-disinfected with 70% ethanol, cut into slices, and dried for 30 min in a bio-safety hood. 25 µL of overnight culture was inoculated on the center of the tissue slices. After inoculation, tissue slices were placed on wet sterilized filter papers in a sterile double-layer drain basket and kept in a growth chamber in conditions of 28 • C and 75% relative humidity. The lesions were measured every day for 3 days. Tissue slices inoculated with LB medium were used as a control.

Zeamine detection assay
The zeamine detection assay was performed according to a method previously described with a few modifications (Chen et al., 2016). In brief, plates were prepared by pouring about 10 mL of LB agar medium, overlaid with 5 mL of 1% agar containing 50 µL overnight culture of fresh Escherichia coli DH5α after solidification. 25 µL of overnight WH1 culture in LS5 medium (per liter contains: 5.25 g K 2 HPO 4 , 2.25 g KH 2 PO 4 , 10.0 g sucrose, 3.6 g NH 4 NO 3 , 1.0 g KCl, 0.25 g MgSO 4 , pH7.0) were added to the wells in plates. The bioassay plates were incubated at 28 • C for 24 h before photography. The diameter of the inhibition zone was measured.

Statistical analysis
All assays were carried out at least three times independently unless otherwise stated. The results obtained were summarized in figures as mean ± standard deviations. Statistical significance was evaluated using a two-tailed Student's t test. A p-value < 0.05 was considered significant.

FIGURE 1
Phylogenetic placement of Dickeya zeae strain WH1. Phylogenetic analysis of Dickeya strains based on concatenated sequences of housekeeping genes, aptD, gyrB, infB, and rpoB, 16 Dickeya strains with complete genomes were included in a phylogenetic analysis with a neighbor joining algorithm and bootstrapped at 1,000 replications.
Frontiers in Microbiology 05 frontiersin.org Results WH1 isolated from rice root was classified as Dickeya zeae based on phylogenetic analysis In this study, although the strain WH1 was isolated from healthy rice root in Wuhu City, China (Eastern China), the pure culture was inoculated to the rice root, and typical soft rot symptoms were observed (data not shown), validating its role as rice pathogen. Meanwhile, we isolated one strain, Paenibacillus polymyxa, from the healthy rice root, which can inhibit the growth of strain WH1 (Supplementary Figure 1). Therefore, it is not surprising that WH1 was isolated from healthy rice roots. Additionally, we found that WH1 can inhibit the growth of Candida albicans, while P. aeruginosa PAO1 can inhibit the growth of WH1 (Supplementary Figure 2). To identify the taxonomic status of strain WH1, MLSA was performed based on the complete sequences of four housekeeping genes including aptD, gyrB, inf B, and rpoB. Phylogenetic analysis of concatenated sequences of four genes distinctly separated D. zeae strains from other species. Among these D. zeae strains, EC1, EC2, and WH1 isolated from rice were grouped into a clade differentiated from other D. zeae strains. WH1 was placed at the same branch with EC2 isolated from rice (Figure 1). ANI analysis indicated that the ANI values between WH1 and previously identified strains from rice ranged from 97.34 to 98.28% and were higher than the cut-off value (96%) for species delineation. Meanwhile, the isDDH values between WH1 and D. zeae strains isolated from rice ranged from 76.3 to 83.9%, which were also higher than the cut-off value (70%) ( Table 1).
In summary, these data establish that strain WH1 belongs to the species D. zeae.

Genomic differences between WH1 and other Dickeya zeae strains
The depth of the assemblies was 243.51 and 322.37 for Illumina and Nanopore platforms, respectively. The complete genome of strain WH1 consisted of a single circular chromosome of 4,675,890 base pairs with a 53.59% GC content. General features of the genome are listed in Table 2, and a detailed circular map is shown in Figure 2. In total, 22 predicted rRNA genes and 76 tRNA genes were identified, and 86.92% of nucleotides are predicted to coding proteins. By a combination of coding protein prediction and homology search, 4096 coding DNA sequences (CDSs) with an average length of 979 bp were identified on the genome. Of these, 274 CDSs were annotated by COG function groups and classified into 26 COG categories. Meanwhile, 1,572 CDSs were annotated by KEGG. 5.81% (238) of total ORFs were annotated as hypothetical proteins. Moreover, four representative D. zeae strains, EC1, EC2, MS1, and MS2, were isolated from rice and banana, resulting in severe soft rot disease. Their genome sequences were downloaded from NCBI and used for genomic comparison analysis. Their genome features are also listed in Table 2.
Pan-genomes analysis was used to develop a robust comparative genomics analysis. The results showed that the number of core gene clusters among the five genomes was 3,040. The number of specific clusters in WH1, EC1, EC2, MS1, and MS2 were 284, 284, 134, 197, and 225, respectively (Figure 3). The unique clusters in WH1 contained some important

Plant cell wall-degrading extracellular enzymes
Most of these extracellular enzymes, pectate lyases (Pel), and other pectinase such as pectin methylestrases (Pem) and pectin lyase (Pnl) produced by D. zeae strains have been shown to play a major role in the virulence and tissue maceration (Boluk et al., 2021). The genome WH1 contained a total of 14 genes encoding pectin degradation enzymes (Supplementary Table 3). These pectin degradation genes were highly conserved in various D. zeae strains, except that pelB encoding pectate disaccharide lyase was absent in the MS1 strain and pemB encoding pectinesterase B precursor was absent in all analyzed strains (Supplementary Table 3). Although the pehKNVWX genes encoding polygalacturonases were also conserved in D. zeae strains, parts of which (pehNVW) were absent in all analyzed genomes. In addition, the endoglucanase genes (celY and celH) and beta-glucosides encoding genes (bglABCD and bgxA) are involved in the degradation of cellulose to glucose. Although these genes were highly conserved in the Dickeya genus, bglC was absent in strains MS1 and MS2 isolated from banana, bglD was just present in EC1, celY was just present in WH1 and EC1. In addition, the homologies of bglB and bglC in WH1 were 53 and 52% with EC1, respectively, although both were isolated from rice. The xylanases (xylBAFGHR and xynA) genes are involved in the degradation of xylan and xyloglucan mainly present in plant cell walls (Boluk et al., 2021). All genomes of analyzed D. zeae strains contained these genes. In summary, these pel, peh, and cel genes were conserved in genomes of D. zeae strains (Supplementary Table 3).
The prtABCDEFG gene cluster encoding four proteases and three protease secretion-associated proteins (T1SS) was located on the positive strand of WH1, EC1, MS1, and MS2 chromosomes, but located on the negative strand of the EC2 chromosome. Alignment with other D. zeae strains showed that they were highly conserved in the genome of WH1 (Supplementary Table 3).

Secretion systems
The type II secretion system (T2SS) is used by gramnegative bacteria to translocate extracellular proteins across the outer membrane (Tan et al., 2021). In D. zeae strains, the T2SS is encoded by the out genes to allow the secretion Venn diagram representing the numbers of common or specific gene clusters among analyzed Dickeya zeae strains in PGAP analysis.
of several proteins including most pectinases. Therefore, T2SS plays an important role in the pathogenicity of Dickeya spp. on host plants. WH1 chromosome harbored a T2SS gene cluster (outOMLKJIHGFEDCBS), covering about 13 kb with 14 ORFs. The gene cluster shared high similarity with that of other analyzed D. zeae strains (similarity 99%) at the amino acid level (Supplementary Table 4). Therefore, the out gene cluster was highly conserved among D. zeae strains.
Besides T2SS plays an important role in the pathogenicity of Dickeya spp., T3SS was also confirmed as a pathogenicity factor in many phytopathogenic bacteria . Similarly, the hrp genes in Dickeya spp. have also been reported to play an important role in their pathogenesis and interaction with host plants (Hu et al., 2018). In the genome of WH1, a large hrp gene cluster spanning 27.3 kb was identified, which was composed of 29 genes, while 27 genes in EC1, three genes in EC2, 26 genes in MS1, and 31 genes in MS2 (Supplementary Table 4). The homologies of hrp gene cluster in WH1 were high with EC1, MS1, and MS2, arranged from 80 to 100% at the amino acid level.
The type VI secretion system (T6SS) was identified to take part in bacterial pathogenicity and host colonization Lin et al., 2017). Many studies showed that T6SS functions in various biological processes, such as mediating cooperative or competitive interactions between bacteria and eukaryotes, and bacterial biofilm formation . T6SS is typically encoded by a cluster of 12-20 genes, with a minimum of 13 genes for the production of a functional apparatus (Zhou et al., 2015). In this study, a gene cluster encoding T6SS was present in D. zeae WH1, spanning 44 kb with 37 ORFs. In this gene cluster, in addition to 17 genes described previously with high identity at the amino acid level, there are 20 additional ORFs inserted between vgrG and tssB genes in genome WH1, while 16 genes in EC1, 22 genes in EC2, 21 genes in MS1, and 20 genes in MS2. For these inserted genes, two genes encoding ankyrin repeat proteins were present in WH1 and EC1, while three genes were in EC2, MS1, and MS2. Additionally, two genes encoding hypothetical proteins were inserted between both genes encoding ankyrin repeat proteins in WH1. IS family transposase was present in EC2 and MS1, but absent in WH1, EC1, and MS2. Two genes encoding plasmid stabilization system proteins ParDE were just found in WH1 and EC1 (Supplementary Table 4). In summary, different genes were inserted between vgrG and tssB of the T6SS gene cluster in all analyzed genomes.

Flagellar and chemotaxis genes
Motility and chemotaxis are essential for phytopathogens when searching for favorable sites to enter into the plant apoplast. Therefore, they play important roles in plant infection by Dickeya spp. The flagellar biosynthesis and chemotaxis clusters constitute 20 fli genes, 14 flg genes, and 5 flh genes, involved in the flagella synthesis, four flagellar rotation genes and six chemotaxis associated genes were present and highly conserved in analyzed D. zeae genomes. Interestingly, two fliC encoding flagellin FliC were found in EC1, MS1, and MS2 genomes, but absent in WH1 and EC2 genomes. Moreover, there were 7 additional genes inserted between fliA and fliC genes in the WH1 genome, encoding a class I SAM-dependent methyltransferase (ctg_02571), a glycosylase (ctg_02573), two hypothetical proteins (ctg_02572 and ctg_02575), a WbqC family protein (ctg_02574), a DegT/DnrJ/EryC1/StrS family aminotransferase (ctg_02576), and a methyltransferase regulatory domain-containing protein (ctg_02577), which were not found in other four D. zeae genomes (Supplementary Table 5). The type IV pilus biogenesis encoding system consists of pilW, pilT, pilABC, and pilMNOPQ genes. The pilW and pilT genes were located distant from the type IV pilus biogenesis cluster. They were highly conserved in D. zeae genomes (Supplementary Table 5).

Polysaccharides biosynthesis genes
Besides lipopolysaccharide (LPS) as a component of the cell wall in Gram-negative strains, Dickeya genus can biosynthesize other polysaccharides, such as capsular polysaccharide (CPS) and exopolysaccharide (EPS). CPS can form a discrete capsule around the cell and is intimately associated with the cell surface; EPS is one of the important virulence factors for bacterial phytopathogens and the main toxic factors leading to waterlogging and wilting on plants after infection (Hu et al., 2018). Both are important factors in infection and biofilm formation. In our analysis, both clusters were present in WH1, MS1, and MS2 genomes, but the CPS cluster was absent in EC1, and the EPS cluster was absent in EC2 (Supplementary Table 6).

Quorum sensing systems
Many Gram-negative bacterial pathogens utilize the luxI/luxR quorum sensing system to regulate the expression of virulence genes and biofilm formation (de Kievit, 2009). Generally, luxI encodes a synthase for the production of AHL family quorum sensing signals, and luxR encodes an AHL signal receptor. Upon interaction with the AHL signal, LuxR becomes an active transcription factor and hence modulating the expression of virulence genes. Previous studies showed that N-3-oxohexanoyl-homoserine lactone (3OC6-HSL), N-3-oxo-octanoyl-homoserine lactone (3OC8-HSL), N-hexanoyl-homoserine lactone (C6-HSL) and N-decanoylhomoserine lactone (C10-HSL) were produced by Dickeya genus (Liu et al., 2022). Of which, 3OC6-HSL was produced by D. zeae EC1 encoding by the luxI homologous expI . A Blast search of the WH1 genome revealed only one copy of expI and one well-conserved luxR homolog expR. ExpI produced by WH1 was highly homologous with the other four D. zeae strains (>99%) at the amino acid level and these four strains produce 3OC6-HSL (Supplementary Table 7). Therefore, WH1 would produce 3OC6-HSL, which was confirmed using HPLC (Data shown below description).
Interestingly, a novel QS system vfm gene cluster upstream of the expI/expR was found in the WH1 and other genomes analyzed in this study (63-100% identity at amino acid level). The gene cluster was originally identified in D. dadantii 3937 isolated from the lesion spot of a wilted African violet and was associated with the regulation of virulence factor production and pathogenesis (Potrykus et al., 2018). The vfm gene cluster, WH1, EC1, EC2, and MS2 genomes contained 26 genes including vfmAZBCDEFGHIJXWVUTSRQPOMNLKY, while the MS1 genome contained 27 genes with one gene encoding fatty acid CoA ligase family protein inserted between vfmN and vfmM (Supplementary Table 7). However, the biological significance of various QS systems in WH1 and other Dickeya species needs to be further investigated.

WH1 produces toxins different from zeamines produced by EC1
Dickeya zeae EC1 produces at least two polyamino phytotoxins and antibiotics, zeamines, which were shown to be the major virulence determinants, affecting rice seed germination and inhibiting the growth of E. coli DH5α (Zhou et al., 2015;Lv et al., 2022). However, no zeamine cluster (zmsO-zmsN) was found in the WH1 and EC2 strains from rice and MS1 and MS2 strains from banana (Supplementary Table 8).
Although zeamine biosynthesis genes were absent in the WH1 genome, the strain still inhibited rice seed germination, and the ability to infect other plants was not reduced (data shown in results of pathogenicity assays). Therefore, the above analysis indicated that the zeamines might not be a necessary virulence factor for this new strain from rice. Other virulence factors could be present in strain WH1 and need to be further investigated.
Phenotype assays WH1 genome harbored the genes encoding proteases (prt cluster), pectate lyases (pel cluster), and cellulases (cel5Z, celH, and celY). The strain produced Prts, Pels, and Cels in the plate assays (Figure 4). Moreover, WH1 generated large swarming (diameter greater than 50 mm) and swimming (diameter greater than 40 mm) diameters when incubated at 28 • C for 20 h (Figures 5A,B). Both CPS and EPS clusters were found in the WH1 genome, which plays important roles in biofilm formation and pathogenicity for D. zeae strains. The results showed that WH1 has the ability to produce exopolysaccharides ( Figure 5C). In addition, WH1 produced thick biofilm in a 96-well plate and the biofilms were verified with fluorescent images (Figure 6). The expI/expR quorum sensing system cluster was highly conserved in the five genomes. Therefore, AHL produced by WH1 was qualified with CV026. The results showed that WH1 produced more AHLs than QS model strain P. aeruginosa PAO1 when incubated at the same conditions (Figure 7). Additionally, the results of HPLC demonstrated that the AHL produced by WH1 was 3OC6-HSL (Figure 8). Acylhomoserine lactones (AHLs) production by Dickeya zeae strain WH1 and Pseudomonas aeruginosa PAO1 at the same conditions. (A) Images are presentative of WH1 and PAO1 grown on Chromobacterium violaceum CV026 lawn for 24 h at 28 • C; (B) quantification of AHLs production abilities, ***p < 0.001.

Pathogenicity assays
The results of pathogenicity tests on potato tubers, carrots, radishes, and Chinese cabbage confirmed that the WH1 strain infected all plants with serious maceration symptoms. When infected for 72 h, almost plant tissues were macerated (Figure 9). We also tested the inhibitory activity of WH1 against rice seed germination. The results indicated that WH1 had a strong inhibitory effect on rice seed germination with an inhibitory rate being about 87% ( Figure 10A). Zeamine detection assay demonstrated that WH1 did not inhibit the growth of E. coli DH5α ( Figure 10B). Therefore, the virulence factors inhibiting rice seed germination were not zeamines.

Discussion
Soft rot Pectobacteriacaea (SRP) includes the genera Pectobacterium and Dickeya and are widespread globally, which cause blackleg disease in potatoes and soft rot disease in a wide range of hosts, including vegetables, fruits, ornamental plants, and the most important crops, rice, and maize. Surprisingly, there is no report on Pectobacterium spp. isolated from rice and maize. In contrast, Dickeya spp. was isolated from a wide range of hosts compared to Pectobacterium spp. . Although there are some differences in isolated hosts, SRP shares many virulence factors, including PCWDE, secretion systems, motility, and plant responses to SRP infections (Davidsson et al., 2013;Reverchon et al., 2016). In Dickeya genus, D. zeae strains and their genomes were used for references to study virulence factors and pathogenicity of the genus (Zhang et al., 2014FIGURE 8 3OC6-HSL produced by Dickeya zeae strain WH1 measured by HPLC. Blue indicates 3OC6-HSL present in the supernatant from D. zeae strain WH1; pink indicates commercial 3OC6-HSL standard.  Zhou et al., 2015;Lv et al., 2022). However, there is a relatively small amount of research on the above diseases in rice (Van Gijsegem et al., 2021).
Research to date indicates that D. zeae strains can infect a wide range of host plants, including 24 dicotyledonous plants and 23 monocotyledonous plants . Therefore, both crop rotation and interplants might be impacted once infected by this kind of strain. Currently, most Dickeya strains were isolated from diseased plants. In this study, although WH1 was initially isolated from the healthy rice root surface, the pure culture has serious pathogenicity on rice and other plants. Additionally, a potential biocontrol strain, P. polymyxa, also isolated from the healthy rice root surface, can significantly inhibit the growth of WH1. Importantly, D. zeae are often present in latent infections on many host crops and they can persist over winter in contaminated plant residues and soil. When the optimal conditions of temperature, humidity, and other factors occur, latency is broken and bacteria proliferate and cause decay (Bez et al., 2021). Therefore, the rice and other plants are grown in the field infected by WH1 face the risk of soft Phytotoxins production by Dickeya zeae strain WH1. (A) Inhibitory activity of toxins from D. zeae strain WH1 against rice seeds germination; (B) the WH1 culture did not inhibit the growth of Escherichia coli DH5α.
rot disease once microbial balance leads to a high abundance of the D. zeae strain.
Based on our analysis, we conclude that the novel strain WH1 is diverse and phylogenetically close to D. zeae strain EC2 isolated from rice in Guangdong, China. Given that the 16S rRNA gene sequence is insufficient to assign the taxonomy of the bacterium, the joint phylogenetic tree based on MLSA analysis was considered. We found that WH1, EC1, and EC2 were grouped into a clade, even though WH1 and EC1, and EC2 were isolated from different places (Eastern China and Southern China). Meanwhile, the results of ANI and DDH analyses also confirmed the MLSA analysis. Interestingly, EC1 harbors a zms gene cluster involving zeamines in biosynthesis, which can inhibit the growth of E. coli DH5α and rice seed germination (Hu et al., 2018), was absent in WH1 and EC2 strains isolated from rice and MS1 and MS2 isolated from banana. However, WH1 produced phytotoxin(s) with high inhibitory activity on rice seed germination and no effect on the growth of E. coli DH5α, which further revealed that WH1 produces other toxins instead of zeamines. Additionally, we found WH1 can inhibit the growth of Candida albicans, while P. aeruginosa PAO1 can inhibit the growth of WH1, which was different from EC1 and MS2 (Hu et al., 2018). Therefore, the toxins produced by WH1 would be different from other D. zeae strains, especially different from EC1. The toxins produced by WH1 require further exploration in the future.
Plant cell wall degrading enzymes produced by the Dickeya genus are considered essential virulence factors for host plants and disease development . Our genomic comparative analyses demonstrated that the main encoding genes associated PCWDEs were existed and conserved in WH1 and other D. zeae strains. Besides, the pathogenesis of Dickeya spp. involved in several complex secretion systems (T1SS-T6SS) to translocate a wide range of extracellular enzymes and effector proteins from the periplasm across the outer membrane. The Prts, which are crucial for virulence, are secreted by the T1SS known as the prtDEF operon. The T1SS cluster, mainly transporting metalloproteinase, was present in all analyzed strains. Previous studies have revealed that Dickeya species mainly used the T2SS to translocate extracellular proteins such as pectinases and cellulases. Similar results were obtained in our analyses. Additionally, Hu et al found that the deletion of T3SS substantially diminished MS2's pathogenicity on both dicots and monocots . Meanwhile, they suggested that the loss of T3SS and reduced PCWDE activity together might have contributed to the host specificity and virulence of D. zeae JZL7, one of the monocot-specific strains isolated from Clivia miniata. The effect of T3SS on WH1's pathogenicity will therefore require further investigation in the future. T6SS was recognized as a distinct class of bacterial protein secretion system and identified as a virulence factor in Gram-negative bacteria . However, the biological function of T6SS in the Dickeya genus has not yet been determined. In this study, a gene cluster encoding T6SS was present in all analyzed genomes.
Among Gram-negative bacteria, EPS and CPS are often essential virulence determinants in plant pathogens. EPS is a main component of the bacterial biofilm matrix and is responsible for adhesion to plant surfaces Zhang et al., 2020). Strains EC1 and EC2 isolated from rice produced either EPS or CPS, while banana isolates MS1 and MS2 generally produced significantly more EPS than strains from rice and exhibited different biofilm-associated phenotypes (Hu et al., 2018). A previous study demonstrated that MS1 could form a thick biofilm, while EC1 biofilm appeared in clumps and EC2 biofilm was thinner . WH1 had the entire EPS and CPS gene clusters, and our study found that the strain easily formed thick biofilm in vitro. Therefore, the relationship between the ability of WH1 to form biofilm and its pathogenicity will be required investigation in the future. Flagellar is used for both swimming and swarming motility. Jahn et al. (2008) have proven that the deletion of the fliA gene encoding a sigma factor obstructed bacterial motility, and limited Pel production and bacterial attachment to plant tissues in D. dadantii 3937. Moreover, the type IV pilus was responsible for motility twitching in P. aeruginosa and Dickeya aquatic (Duprey et al., 2019). Therefore, flagellar biosynthesis and chemotaxis proteins and the type IV pilus play important roles in plant bacterial pathogenicity. We found these clusters present in all analyzed D. zeae strains.
Dickeya species have at least two QS systems, expI/expR, and vfm systems. A previous study demonstrated that deletion of expI in D. zeae EC1 abolished AHL production, increased swimming and swarming motility, disabled biofilm formation, attenuated virulence of pathogen on potato tubers, and reduced Frontiers in Microbiology 13 frontiersin.org the inhibitory activity on rice seed germination through decreasing the expression of zms genes (Hussain et al., 2008). The inactivation of expI in D. zeae MS2 affected motility and cell clumping, while it did not affect infections of banana seedlings (Feng et al., 2019). Therefore, the effects of the expI/expR system on the virulence and pathogenicity of D. zeae strains is complex. The genome of WH1 also contains expI/expR clusters with an identity of 99% with other D. zeae strains. Meanwhile, HPLC demonstrated the AHL produced by WH1 was 3OC6-HSL. However, there are different virulence determinants between WH1 and EC1, such as the zeamine biosynthesis gene cluster that is present in EC1 but absent in WH1. Therefore, the effect of expI/expR system on virulence factors production of WH1 remains to further investigate in the future. Vfm-QS signal is uniquely present in Dickeya strains, participating in and dominating the macerating soft rot process of the pathogen in plant tissues. Previous studies demonstrated that the vfm QS system participated in the regulation of PCWDEs production and virulence in the Dickeya genus (Zhang et al., 2014). However, the chemical structure of the Vfm signal has not been characterized yet. In this study, the entire vfm gene cluster was found in the WH1 genome, but its effect on virulence has not yet been investigated through experimental evidence.

Conclusion
This study presents a high-quality complete genome sequence of a novel D. zeae strain WH1 isolated from rice root. Phylogenetic tree analysis and genome-genome nucleotide comparison indicated that WH1 was placed at the same branch with EC2 isolated from rice. The genome sequence of WH1 was annotated and compared with the representative genome sequences of other D. zeae strains, focusing on virulence determinants and potential regulatory mechanisms. Several groups of virulence genes, such as those encoding for PCWDEs, T1SS gene cluster, T2SS gene cluster, T6SS gene cluster, flagellar and chemotaxis gene cluster, and quorum sensing gene cluster were highly or fully conserved in all five genomes isolated from different hosts. Interestingly, the T3SS gene cluster was absent in strain EC2, T4SS gene cluster was absent in WH1, EC2, and MS2. Genes inserted between vgrG and tssB within the T6SS gene cluster were different in all analyzed genomes. Both EPS and CPS gene clusters were present in WH1, MS1, and MS2, while either EPS or CPS was present in EC1 or EC2. Gene cluster encoding zeamines biosynthesis was only present in EC1. Additionally, related phenotypic and pathogenicity analyses were performed for strain WH1. The strain can produce pectate lyases, cellulases, and proteases, easily formed thick biofilm in vitro, and had strong AHL quorum sensing activity. Besides, WH1 can infect plants with serious maceration symptoms and strongly inhibit rice seed germination through other toxins instead of zeamines.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi. nlm.nih.gov/PRJNA858112, CP101400.


s2fieldsofstudy: Biology Medicine

text:
Bilirubin ameliorates murine atherosclerosis through inhibiting cholesterol synthesis and reshaping the immune system

Atherosclerosis is a chronic inflammatory disease caused mainly by lipid accumulation and excessive inflammatory immune response. Although the lipid-lowering and cardioprotective properties of bilirubin, as well as the negative relationship between bilirubin and atherosclerosis, were well documented, it is not yet clear whether bilirubin can attenuate atherosclerosis in vivo. In this study, we investigated the role of bilirubin in improving atherosclerosis. We found that mildly elevated bilirubin significantly reduced the risk factors of atherosclerosis, such as plasma glucose, total cholesterol, and low-density lipoprotein cholesterol, and the formation of atherosclerotic plaques, liver total cholesterol, and cholesterol ester concentration in apolipoprotein E-deficient (ApoE−/−) mice fed a western-type (high fat) diet. It was further found that bilirubin could promote the degradation of 3-Hydroxy-3-Methylglutaryl-CoA Reductase (HMGCR), a rate-limiting enzyme for endogenous cholesterol synthesis. Using mass cytometry-based high dimensional single cell analysis, we observed a decrease of natural killer cells and an increase of dendritic cells and myeloid-derived suppressor cells, which all are closely associated with atherosclerosis risk factors and contribute to the improvement of atherosclerosis, in ApoE−/− mice treated with bilirubin. By in-depth analysis, modulation of multiple spleen or peripheral blood T cell clusters exhibiting either positive or negative correlations with total cholesterol or low-density lipoprotein cholesterol was detected after bilirubin treatment. In this study, we demonstrate that bilirubin serves as a negative regulator of atherosclerosis and reduces atherosclerosis by inhibiting cholesterol synthesis and modulating the immune system. Supplementary Information The online version contains supplementary material available at 10.1186/s12967-021-03207-4.

Introduction
Cardiovascular disease (CVD) has become the leading cause of death and morbidity worldwide [1]. Atherosclerosis is the principal cause of CVDs, cerebral infarctions, and peripheral vascular diseases. Atherosclerosis is defined as a chronic inflammatory disease of the arteries. Plaque decomposition or plaque rupture can block arteries, causing severe acute cardiovascular events such as myocardial infarction, stroke, and sudden death [2][3][4][5].
Bilirubin is a final product of heme catabolism. In microsomes, heme is converted into carbon dioxide, iron, and biliverdin IX (biliverdin) and thereafter biliverdin reductase transforms biliverdin into bilirubin [6]. The physiological concentration of bilirubin is beneficial for oxidative stress-related pathology, a key event in the development of atherosclerosis [7][8][9]. Bilirubin levels have been shown to be inversely correlated with CVD mortality and severity [10][11][12][13][14][15][16]. Numerous studies have shown that bilirubin offers cardioprotective properties in Gilbert syndrome (GS) patients, who suffer sustained unconjugated hyperbilirubinemia. A low prevalence (2%) of ischemic heart disease (IHD) was found in GS patients, as compared to a higher prevalence (12.1%) among the general population [17]. Hyperbilirubinemia in GS subjects or Gunn rats causes a decrease in circulating cholesterol, triacylglycerol, and oxidized low-density lipoprotein (LDL) by increasing the serum anti-oxidative functions [18,19]. Bilirubin has the potential to lower lipid, which is associated with cardiovascular protection in GS, as well as dyslipidemia, one of the leading risk factors for atherosclerosis [20]. Oxidative stress plays a critical role in the pathogenesis of atherosclerosis. Bilirubin, a natural endogenous antioxidant, protects vascular cells from oxygen radical damage, inhibits oxidative modification of LDL and improves lipoprotein composition [21,22]. Bilirubin also inhibits inflammation by downregulating the expression of endothelial cell adhesion molecules, preventing leukocyte adhesion, rolling, and infiltration into vessels [23][24][25]. Additionally, bilirubin reduces neoplastic intima formation by inhibiting smooth muscle cell proliferation and migration. Because of these properties, bilirubin is involved in every aspects of atherosclerosis, inhibiting plaque formation and promoting plaque stabilization [26,27]. These findings greatly demonstrate the importance of bilirubin in the regulation of lipids and improvement of atherosclerosis. To date, the underlying mechanisms of hypocholesterolemia in hyperbilirubinemia have not been fully elucidated. Atherosclerosis is a disease associated with inflammation that relies heavily on the immune system for its development and modulation [28]. New high-throughput single cell technologies, such as mass cytometry and single-cell RNA sequencing, enable a comprehensive analysis of the immune system during atherosclerosis progression through uncovering the diversity of > 15 immune cell populations [29][30][31]. The dissection of immune modulation has greatly contributed to the understanding of atherosclerosis mechanisms and has provided novel insights into its treatment via immune therapy. The hypolipidemic effect of bilirubin led us to examine its effect on the development of atherosclerosis in an apolipoprotein E-deficient (ApoE −/− ) mouse model, as well as the mechanisms involved from the perspective of lipid regulation and immunomodulation using mass cytometry-based single cell analysis. In the present study, we found that bilirubin negatively regulates atherosclerosis through inhibiting cholesterol synthesis and modulating the immune system.

Experimental animals
Animal experiments in our study were approved by the Institutional Committee for the Use and Care of Laboratory Animals of Guangzhou Medical University (2017-014) and performed in accordance with the guidelines from Directive 2010/63/EU of the European Parliament on the protection of animals used for scientific purposes or the NIH guidelines. 8-week-old male ApoE −/− mice (C57BL/6.129P2-APOE/J, Beijing Vital River Laboratory Animal Technology Co., Ltd) were divided into two subgroups. One group (Bilirubin, n = 14) received an intraperitoneal injection of bilirubin (20 mg/kg/d, B4126, Sigma) every other day for 10 weeks. Mice injected with the vehicle served as a control (Control, n = 8). Both groups started the western-type diet supplemented with 21% fat (wt/wt) and 0.15% cholesterol (wt/wt), simultaneously receiving the injection of bilirubin or vehicle. All mice were maintained on these diets for 10 weeks. The body weight, diet, and drinking of mice were monitored for a continuous week. The mice were euthanized with an overdose of pentobarbital sodium (200 mg/Kg) for tissue collection.

Biochemical assays
Peripheral blood was collected into a heparin-coated tube by cardiac puncture and centrifuged to isolate plasma. Total bilirubin (T-BIL) in the plasma was determined by the vanadic acid oxidation method. Liver enzymes, including alanine aminotransferase (ALT) and aspartate aminotransferase (AST), were measured by the ultraviolet lactate and malate dehydrogenase methods. Total cholesterol and triacylglycerol levels were measured by an oxidase method. High-density cholesterol (HDL-C) and low-density cholesterol were measured by a direct method. The total bilirubin assay kit was from FUJIFILM Wako pure Chemical Corporation, Japan. The kits for glucose, liver enzymes, lipids, and lipoproteins were from Kehua Bio-Engineering, China. All the spectrophotometric methods were performed according to the manufacture's instructions using an autoanalyzer.

Histological and morphometric analysis
The heart was perfused with phosphate-buffered saline. Whole aorta from the heart, extending 5-10 mm after bifurcation of the iliac, including the brachiocephalic trunk, left common carotid artery, and subclavian arteries, were dissected free of fat tissues and opened longitudinally before being fixed in formalin solution at room temperature (RT) for 24 h. After staining with Oil Red O (ORO), aortae were pinned on a black wax plate and photographed under stereoscopy at standardized magnification and illumination. The ORO positive area and entire aorta tree area were quantified using Image J software. Aorta lesion areas were expressed as a percentage of the whole aorta tree area. The aortic root and brachiocephalic trunk were serially sectioned into 5 μm cryosections. The first section of the aortic sinus was harvested and collected in pairs on glass slides when all three aortic valves became visible in the lumen of the aorta. Neutral lipids deposition in the lesion was visualized by hematoxylin and eosin (H&E) or ORO staining, respectively. The collagen of the lesion was detected by Masson's Trichrome staining. The liver was dissected and weighed, expressed as a ratio of liver weight to body weight. Liver sections were collected and performed H&E or ORO staining for histological and lipid deposition analysis.

Immunofluorescence staining
The contents of macrophages and smooth muscle cells in aortic sinus lesions were determined by staining with anti-Galectin 3(ab53082) and anti-α-SMA (ab5694) and detecting with Alexa-conjugated antibodies (ab150073 and ab150080). Nuclei were counterstained with DAPI. All antibodies were purchased from Abcam.

Total cholesterol, cholesterol ester, and free cholesterol assay
Hepatic total cholesterol and free cholesterol concentration were measured using a tissue total cholesterol assay kit (E1015; Applygen Technologies, Inc.) and free cholesterol assay kit (E1016; Applygen Technologies, Inc.) following the manufacturer's instructions. Cholesteryl ester was quantified by subtracting the free cholesterol values from the total cholesterol value.

Western blot
Cells or tissues were collected and homogenized in icecold RIPA buffer. Lysates were maintained constant agitation for 2 h at 4 °C. Gently aspirate the supernatant to a prechilled tube after the centrifugation at 13,200 rpm for 10 min at 4 °C. The protein level of lysates was determined by the BCA protein assay. An aliquot of proteins was separated by SDS-PAGE and then transferred to the immunoblot PVDF membrane. The immunoblots were incubated with primary antibodies at 4 °C overnight, followed by HRP-conjugated secondary antibodies (1:5000) at RT for 1 h. The signals of targeted proteins were visualized by employing the ECL detection system (Amersham). Primary antibodies used for immunoblots were rabbit polyclonal antibodies against HMGCR (ab174830, diluted 1:1000) and mouse monoclonal antibody against GAPDH (ab8245, diluted 1:1000). All antibodies were purchased from Abcam.

Co-immunoprecipitation (Co-IP)
Cells were collected and treated with IP lysis buffer containing protease inhibitor cocktails at 4 °C for 30 min. Cell lysates were centrifugated at 13,200 rpm for 10 min at 4 °C to remove cell aggregates. 50 μl of protein extract were dissolved in a DTT-containing loading buffer as an input sample. The remaining protein extract was adjusted to a concentration of 2 μg/μl of Co-IP dilution buffer. 500 μl of protein extract were immunoprecipitated with magnetic beads conjugated with indicated antibodies at RT for 1-2 h in a roller. The beads were washed with prechilled Co-IP washing buffer three times and then heated in 2 × loading buffer for 10 min at 70 °C. After centrifugation at 12,000 rpm for 10 min at 4 °C, the beads were collected by inserting the protein samples in a magnetic stack for 10 min at RT. The well-prepared protein samples were analyzed by western blot. The immunoblots were probed with rabbit polyclonal antibody against HMGCR (ab174830, diluted 1:1000), rabbit polyclonal antibody against ubiquitin (ab179434, diluted 1:5000), and mouse monoclonal antibody against GAPDH (ab8245, diluted 1:1000). All antibodies were purchased from Abcam.

Protein turn-over assay
To evaluate protein half-life, we treated LO2 and HepG2 cells with or without bilirubin (12 μM) in the presence of cycloheximide at a final concentration of 50 µg/ml to block protein synthesis for the indicated time (0, 6, 12, and 24 h). Protein samples were subjected to western blot to analyze the expression of HMGCR.

Peripheral blood and spleen cell preparation
The peripheral blood was collected from the heart after the injection of pentobarbital sodium. After mice were sacrificed, the spleen was isolated and spleen cells suspension was separated by gently crushing the spleen. Cells from peripheral blood or spleen were fixed with Fix I buffer (Fluidigm) for 10 min at RT. After washing with cold PBS, red blood cell lysis buffer was used to remove erythrocytes from peripheral blood or spleen cells.

Cell immunostaining
Peripheral blood or spleen cells were barcoded separately with premade combinations of six different palladium isotopes using the Cell-ID 20-Plex Pd Barcoding Kit (Fluidigm) for 30 min at RT. The samples were washed twice with cell staining buffer and then pooled together. The pooled samples were incubated with anti-CD16/32 (FcR III/II, Biolegend, CA, USA) for 10 min at RT to lower nonspecific antibody binding. These cells were washed twice with cell staining buffer and stained with a cocktail of 26 metal isotope-conjugated antibodies (Additional file 1: Table S1) for 30 min at RT. Stained cells were washed and stained with Intercalator-Ir (Fluidigm) at 4 °C overnight. After washing with cell staining buffer and ddH 2 O, cells were resuspended in ultrapure water supplemented with EQ Four Element Calibration Beads (Fluidigm) and analyzed using Helios mass cytometer (Fluidigm).

Mass cytometry data processing and analysis
The resulting flow cytometry standard (FCS) files generated by mass cytometry were normalized, randomized and debarcoded using the CyTOF software (Fluidigm). The debarcoded files were uploaded to cytobank.cn, an online platform for cytometry-based single cell analysis. The single cell data without debris and doublets were used for high-dimensional analyses, such as viSNE and FlowSOM.

Statistical analysis
The Shapiro-Wilk test was used to determine whether the data was normally distributed. A parametric test (Student's t test) was used to determine the statistical significance between 2 groups if the data is normally distributed, otherwise, a nonparametric test (Mann-Whitney test) was used. Correlation analyses were performed using Person's correlation coefficient. Error bars represent the mean ± standard deviation (SD). p < 0.05 was regarded as statistically significant.

Bilirubin lowers the risk of atherogenesis in ApoE −/− mice
To investigate whether bilirubin can protect atherogenesis in vivo, we used a murine model in which ApoE −/− mice were fed on a western-type (high fat) diet. Mice received an intraperitoneal injection of bilirubin or vehicle every other day for 10 weeks. There was no significant difference between the average body weight of the two groups of mice throughout the 10 weeks of intervention. (Fig. 1A). In comparison with a vehicle control group, exogenous bilirubin administered resulted in a fourfold increase in total plasma bilirubin ( Fig. 1B). Treatment with bilirubin reduced plasma glucose, total cholesterol, and LDL cholesterol levels in mice. (Fig. 1C and D), whereas both groups had similar levels of triglycerides and high-density lipoprotein (HDL) cholesterol. Furthermore, multiple correlations were found between glucose, total cholesterol, LDL, HDL, triglyceride, and total bilirubin. (Fig. 1E). For example, total bilirubin levels were negatively and significantly correlated with the levels of glucose and LDL levels in the peripheral blood. The LDL level was positively and significantly correlated with the total cholesterol level. Mice with higher levels of glucose also had higher levels of total cholesterol and LDL. There were significant and positive correlations between the concentration of triglyceride and total cholesterol or HDL (Fig. 1E). Considering the long-term exposure, we also investigated the bilirubin toxicity to the liver through histological and biochemical analysis. Bilirubin-treated mice displayed a downward trend in a lower hepatic enzyme activity (AST and ALT) than control mice (Additional file 1: Fig. S1B). But both groups showed no change in liver weight/body weight. (Additional file 1: Fig. S1A). These findings illustrate that a mildly elevated bilirubin reduced the risks of developing atherosclerosis without causing liver toxicity, suggesting a protective role of bilirubin in atherogenesis.

Bilirubin reduces atherosclerosis in ApoE −/− mice
In vivo, bilirubin demonstrated a risk-lowering effect, thus we investigated whether it affects atherosclerotic lesion formation. Morphometric assays were used to measure the aortae and the aortic sinus lesion formation after 10 weeks of intervention. In both groups of mice, lesions were present throughout the aorta, but lesions were reduced in mice treated with bilirubin, as determined by en face ORO staining. (Fig. 2A). The bilirubin treatment significantly decreased the size and lipid load of the aortic sinus plaque and the aortic sinus lesion area. (Fig. 2B). Mice treated with bilirubin also showed reduced brachiocephalic arteries (Fig. 2C). Additionally, bilirubin treatment promoted collagen deposition ( Fig. 2B) and α-SMA ( +) smooth muscle cell accumulation without affecting the distribution of Galectin-3( +) macrophages (Fig. 2D). The liver plays a crucial role in lipid metabolism. In the liver, cholesterol is derived from dietary, de novo synthesis, and transport from extrahepatic tissues. In mice fed a western diet for a long time, excessive lipid accumulates in the liver, increasing the risk of atherosclerosis and liver steatosis. H&E and ORO staining showed that bilirubin-treated mice had less hepatic steatosis than vehicle-treated mice (Fig. 2E). Next, The liver was then analyzed using an enzymatic method to determine the content of total cholesterol, cholesterol ester, and free cholesterol. The bilirubin treatment decreased liver total cholesterol and cholesterol ester concentration, but did not affect free cholesterol (Fig. 2F). As a result, bilirubin significantly inhibits the formation of atherosclerotic plaques and favors the phenotype of stable plaque.

Bilirubin promotes HMGCR protein degradation
The liver is the primary organ that synthesizes cholesterol and around 20-25% of total daily cholesterol production occurs here. Since HMGCR is the key ratelimiting enzyme for endogenous cholesterol synthesis, we examined HMGCR expression in vivo and in vitro. HMGCR expression in the liver was significantly reduced in ApoE −/− mice treated with bilirubin compared to vehicle-treated mice ( Fig. 3A and Additional file 1: S2A). In some cases, two bands of HMGCR have been detected in our and other studies [32,33]. The upper one is the HMGCR band and the lower one is the non-specific band. The expression of HMGCR in the liver was also decreased in BALB/c mice treated with bilirubin for three consecutive days (Additional file 1: Fig. S2B). In two hepatic cell lines, bilirubin did not decrease the mRNA expression of HMGCR in vitro (Fig. 3B), but decreased the protein level of HMGCR in a time-and dose-dependent manner ( Fig. 3C and D), suggesting that bilirubin may affect HMGCR expression at the translational or post-translational level. After treating two hepatic cell lines with cycloheximide (CHX), we found that bilirubin treatment decreased protein synthesis (Fig. 3E). Considering that HMGCR is degraded through ubiquitination at lysine 89 and 248 [34], we further evaluated the effect of bilirubin on the ubiquitination of HMGCR. Bilirubin increased HMGCR ubiquitination in HepG2 hepatic cells treated with MG132, a proteasome inhibitor (Fig. 3F). Furthermore, bortezomib, a potent protease inhibitor, reversed the decrease in HMGCR levels induced by bilirubin in two hepatic cell lines (Fig. 3G). These results indicate that bilirubin reduces HMGCR level mainly by promoting its ubiquitination and degradation, consequently reducing endogenous cholesterol synthesis.

Bilirubin affects immune system in the spleen
Since bilirubin significantly improves atherosclerosis which is closely associated with the immune response, it is likely that bilirubin also has an effect on the immune system. Thus, we used a high-throughput single cell technology, mass cytometry [31,35], to simultaneously measure the expression of 26 markers on immune cells from the spleen or peripheral blood at the single cell level. ViSNE analysis of CD45 + cells from the spleen was performed to visualize these high-dimensional mass cytometry data in two dimensions. According to the expression of 13 surface markers on the cells displayed on the viSNE map ( Fig. 4A), 18 major immune cell populations were identified and gated (Fig. 4B). In each gated population, the expression of 13 markers was the same as the phenotypes of the indicated types of immune cells (Fig. 4C). Next, we compared the proportions of all these cell populations and found that after bilirubin treatment, the percentages of CD11b low dendritic cells (CD11b low DCs) significantly increased, while the proportion CD11b high DCs markedly decreased (Fig. 4D). In addition, the percentage of CD4 T cells in the spleen T cells was significantly reduced, while the proportion of CD8 T cells increased (Fig. 4D), suggesting that bilirubin induced a change in spleen T cell composition. Interestingly, the proportion of CD4 T cells among spleen T cells was strongly, significantly, and positively related to the LDL and (total cholesterol) TCHO concentration in the blood, while CD8T was strongly, significantly, and negatively related to the LDL and TCHO concentration in blood bilirubin (Fig. 4E), indicating that spleen CD4 T cells are positively associated with atherosclerosis. The proportion of CD11b high DC that decreased by bilirubin showed a strong and negatively relationship with total bilirubin (Fig. 4E). Moreover, the proportion of CD8 T cells in spleen T cells, as well as the percentages of CD11b low DC, were significantly and negatively associated with the LDL and TCHO concentrations in the blood (Fig. 4E), implicating that these immune cells may favor the improvement of atherosclerosis.

Bilirubin alters the composition of T cells in the spleen
T cells can influence the process of atherosclerosis by regulating their anti-inflammatory and pro-inflammatory functions and their functions are mainly determined by the direction of their differentiation. T cells are highly heterogeneous and contain various subtypes of differentiated T cells. Thus, we used 12 T cell-related markers to systematically analyze the spleen T cell. After viSNE analysis, heterogeneity in markers was assessed at the single-cell level. From the viSNE map of T cells, we also observed a huge heterogeneity of the T-cell compartments (Additional file 1: Fig. S3A). FlowSOM, an unsupervised cluster analysis for visualizing and interpreting cytometry data, was used on these T cells to map T cell phenotypes exhaustively and to further define any exclusive or significantly altered T cell clusters. 30 T cell clusters containing cells with similar phenotypes were identified and visualized on the viSNE map (Fig. 5A). The expression profiles of all the T cell clusters were visualized in a heatmap (Fig. 5B) and this approach led to the identification of 12 CD4 + clusters, one double-positive cluster, eight CD8 + clusters, five double-negative phenotypes, one natural killer (NK) T cluster and three TCRβ − clusters (Fig. 5B). Among these 30 T cell clusters, proportions of T5 (CD44 + CD62L − CXCR3 + CD4T), which mainly contains effector memory (em) Th1 cells, and T16 (Ly6C + CD44 + CD62L + CD8T), which mainly contains central memory Ly6C + CD8 T (Tcm) cells, were significantly increased after bilirubin treatment (Fig. 5C). In addition, these two clusters, together with T23 which was increased with a trend close to significance, had a mutual, significant, and negative relationship with the concentrations of LDL and TCHO in the blood ( Fig. 5D and E). Furthermore, T6 and T13 proportions decreased with a trend close to significance after bilirubin treatment (Fig. 5C). Bilirubin did not affect the proportion of T1 (CD3 + CD4 + CD25 + CD127 − Ly6C + ) and T2 (CD3 + CD4 + CD25 + CD127 − Ly6C − ) clusters which mainly contain regulatory T (Treg) cells (Additional file 1: Fig. S3B). These data indicate that bilirubin also can regulate the composition of spleen T cells which are related to the improvement of atherosclerosis.

Bilirubin regulates myeloid-derived suppressor cells (MDSCs), DCs and NK cells in the peripheral blood
After bilirubin administration, the immune cells in the peripheral blood were also comprehensively analyzed using mass cytometry. viSNE was also used to identify the cell populations (Fig. 6A). 19 immune cell populations were obtained and the expression of their phenotypic markers was confirmed ( Fig. 6B and C). The percentages of these 19 cell populations were then compared. Proportions of DC, Ly6C + NK, and all NK cells were significantly reduced in the peripheral blood of mice treated with bilirubin (Fig. 6D). The percentages of polymorphonuclear (PMN)-MDSC, neutrophils, and all MDSC significantly increased (Fig. 6D) and were strongly, significantly, and positively correlated with bilirubin concentration (Fig. 6E and F), implicating that bilirubin treatment may directly increase these cells. In addition, the proportion of all MDSC was significantly and negatively correlated with the concentrations of LDL and TCHO in the blood (Additional file 1: Fig. S4A). Considering that higher levels of LDL and TCHO are closely associated with atherogenesis, MDSCs increase by bilirubin may contribute to the improvement of atherogenesis though regulating the levels of LDL and TCHO. Moreover, proportions of Ly6C + NK and all NK cells showed a significant negative correlation with blood bilirubin level (Fig. 6G), suggesting that the reduction of these NK cells may be directly caused by bilirubin administration. The proportion of Ly6C + NK cells correlated positively and significantly correlated with the level of LDL in the blood (Additional file 1: Fig. S4B). Furthermore, B cells were negatively correlated with the concentration of bilirubin and LDL, and positively associated with TCHO concentration (Fig. 6F and Additional file 1: Fig. S4C). These data propose that bilirubin affects the MDSC and NK cells which are closely associated with the improved atherosclerotic blood lipid profiles.

Bilirubin mildly affects T cell composition in the peripheral blood
Peripheral blood T cell composition after bilirubin treatment was also analyzed using viSNE (Additional file 1: Fig. S5A). 29 T cell clusters, including 13 CD4 + clusters, eight CD8 + clusters, four NKT clusters, two double-negative phenotypes, and two TCRβ − clusters were identified ( Fig. 7A and B). Treg cell clusters, T1 and T2, were not changed by bilirubin treatment (Additional file 1: Fig. S5B). T7 (CD44 + CCR4 + CD4 T), T9 (CD44 + CCR4 + CCR6 + CD4 T), T13 (CD44 − CD62L + CD4 T), and T17 (Ly6C + CD44 − CD62L + CD8 T) were significantly decreased, whereas T20 (Ly6C − CD44 − CD62L − CD8 T) was increased after bilirubin treatment (Fig. 7C). Additionally, proportions of two TCRβ − clusters, T28 and T29, were also changed by bilirubin with a trend close to significance (Additional file 1: Fig. S5C). The percentages of T9 and T28 showed a negative and significant correlation with the total bilirubin level, whereas T20 and T29 showed a positive correlation with total bilirubin (Fig. 7D and E and Additional file 1: S5C). There was a positive and significant correlation between the proportions of T7, T9, T13, and T17 and the amount of LDL in the blood. (Fig. 7F). Except for T17, the other clusters were also positively correlated with TCHO concentration (Additional file 1: Fig. S5D). These results suggest that bilirubin can induce changes of T cell subclusters and that these changed T cells also have a strong correlation with the improvement of atherosclerosis.

Discussion
Accumulating epidemiological evidence has demonstrated the negative relationship between bilirubin and CVD. Bilirubin, a potential predictor and protective marker of CVD, might become a new therapeutic target in the clinical setting. However, little is known about the direct cardioprotective effects of bilirubin. In the present studies, exogenous bilirubin inhibited lesion formation and lowered circulating cholesterol, LDL cholesterol, glucose concentrations, and hepatic lipid accumulation in ApoE −/− mice receiving western-type diet. We confirmed that bilirubin prompted HMGCR degradation via the proteasome pathway, suggesting that the cardioprotective effects of bilirubin may also be governed by its lipidlowering properties. Furthermore, bilirubin significantly modulated the immune system, which resulted in a cardioprotective immune profile. These findings highlight the negative regulatory role of bilirubin in atherosclerosis while suggesting bilirubin as a potential treatment option for atherosclerosis. Dyslipidemia and diabetes are common risk factors for atherosclerosis. A higher level of circulating LDL cholesterol or total cholesterol is independently and positively correlated with the prevalence of atherosclerotic peripheral vessel diseases [36]. In our study, there was a significant reduction of total cholesterol (− 17%), LDL-C (− 28%) and blood glucose (− 26%) in mice treated with bilirubin, confirming the lipid-lowering properties of bilirubin. Studies in animals (Gunn rats) or humans (GS subjects) of hyperbilirubinemia have consistently shown that bilirubin protected against IHD and atherosclerosis in dyslipidemia or diabetes high-risk population [17,19,37,38]. In light of this circumstantial evidence, along with our finding showing that bilirubin attenuates atherosclerotic lesion formation and promotes a stable atherosclerotic plaque, further confirm that bilirubin is a negative regulator in atherosclerosis.
Bilirubin protects against atherosclerosis due to its anti-oxidative and anti-inflammatory properties, as well as the regulation of leukocyte migration mediated by endothelial adhesive molecules [7]. In atherosclerosis, foam cells are formed by the uptake of oxidatively modified lipids as oxLDL by macrophages [39]. Bilirubin has been shown to inhibit lipoprotein oxidation and therefore protects endothelial cells from being damaged by oxidative stress [22]. McNamara et al. showed that the serum remnant-like particle cholesterol (RLP-c), an independent risk factor for CVD, was higher in the male patients in the high bilirubin group than those in the low bilirubin group [40]. It is showed that bilirubin regulates cholesterol and lipid metabolism via activating the aryl hydrocarbon receptor (AhR) signaling pathway, which has been shown to protect against cardiogenesis and oxidative stress [41,42]. Notably, bilirubin, as a novel PPARα agonist, suppressed lipid accumulation in 3T3-L1 adipocytes and PPARα KO mice [43]. According to a recent study, increased unconjugated bilirubin promotes the clearance of excess cholesterol out of the body by promoting transintestinal cholesterol secretion [18]. These previous studies have shown that bilirubin can regulate LDL via multiple pathways, which is also in accordance with our findings. The liver plays a critical role in maintaining cholesterol homeostasis by regulating multiple mechanisms, including exogenous cholesterol uptake via LDLR and cholesterol biosynthesis via HMGCR activity, esterification for storage, and reverse cholesterol transport (RCT) [44]. A disruption in any one of these pathways can lead to dyslipidemia. Newly synthesized cholesterol assembles into very low-density lipoprotein and triglycerides, which is then excreted into the bloodstream or stored in lipid droplets. Ectopic lipid accumulation and the deposition of lipid within non-adipose tissues (e.g., muscle and liver) contribute to the pathogenesis of metabolic diseases, such as diabetes and atherosclerosis [45,46]. Hepatic steatosis, an indicator of atherosclerosis, can be caused by excessive lipid deposits in the liver. We found that ApoE −/− mice treated with bilirubin displayed a decrease in lipid droplets accumulation in the liver. In addition, bilirubin also inhibited cholesterol and cholesterol esters production, suggesting that bilirubin attenuates hepatic steatosis and atherosclerosis by interrupting cholesterol biosynthesis. Our findings are consistent with a previous study, which demonstrated mild hyperbilirubinemia conditions significantly reduced the intracellular lipid accumulation in C2C12 skeletal mouse muscle and HepG2 human liver cells [47].
HMGCR, a key rate-limiting enzyme, mediates liver de novo biosynthesis of cholesterol. Statins can significantly suppress cholesterol biosynthesis through inhibiting HMGCR activity and thus are widely used in treating hypercholesterolemia. Two negative feedback pathways are available to regulate the expression of target genes involved in cholesterol metabolism: the processing of sterol regulatory element-binding proteins 2 (SREBP2), a crucial transcriptional regulator of cholesterol biosynthesis, and sterol-induced HMGCR degradation [48]. Here, bilirubin did not reduce the mRNA expression of HMGCR, but significantly decreased its protein level, suggesting the involvement of a post-translational mechanism. Song et al. reported that the mTORC1-USP20-HMGCR axis, a post-transcriptional mechanism, participated in the turnover of HMGCR and that the inhibition of USP20 might be beneficial for a variety of metabolic diseases, such as diabetes, hyperlipidemia, and CVD [48]. Our findings demonstrated that bilirubin promotes the ubiquitination and degradation of HMGCR, thus contributing to a reduction of cholesterol biosynthesis and the improvement of atherosclerosis.
Atherosclerosis is considered a chronic inflammatory disease with autoimmune responses. Both the innate and adaptive immune responses have been implicated in atherosclerosis progression. Here, we found that bilirubin displays atheroprotective properties in ApoE −/− mice and acts as a regulatory role on both innate and adaptive immune cells. DCs present antigens to T cells and thus provide an important bridge between innate and adaptive immunity. DCs involves atherosclerosis through multiple mechanisms, inducing lipid uptake, antigen presentation, and cytokine production [49,50]. There are three DC subsets, including conventional DC (cDC) 1, cDC2, and plasmacytoid DC [51]. cDC1, which express low levels of CD11b, play atheroprotective function through inducing Treg cells [52]. We showed that bilirubin increases spleen CD11b low DCs, which correlates negatively with LDL and TCHO levels. Thus, bilirubin likely exerts atheroprotective effects by regulating the functions of these DCs.
CD4 + T cells are important mediators of the adaptive immune response and their subtypes, Th1 cells, have pro-atherogenic functions mainly through secreting proinflammatory cytokines, such as IFN-γ, IL-2, IL-3, and tumor necrosis factor [53,54]. In the present study, CD4 T cells were significantly reduced by bilirubin treatment and this decrease is positively associated with LDL and (See figure on next page.) Fig. 6 Bilirubin reshapes the immune system in the peripheral blood. ApoE −/− mice fed with a western-type diet were intraperitoneally injected with (bilirubin, n = 14) or without (control, n = 8) bilirubin. The animals were sacrificed after ten weeks. The peripheral white blood (PWB) cells were collected and simultaneously stained with 26 metal isotope-labeled antibodies. The expression of 26 surface markers was then determined by mass cytometry. TCHO, suggesting the involvement of these CD4 T cells in atheroprotective role of bilirubin. Moreover, after indepth single cell analysis, a decrease of T6, which mainly are Th1 Tem and accounts for above 15% of spleen T cells, contributes to the overall reduction of CD4 T cells.
Considering the pro-atherogenic functions of Th1 cells, decrease of effector memory Th1 cell may be one of the causes of bilirubin-induced atherosclerosis improvement. NK cells, a critical component of the innate immune system, have been demonstrated as an atherogenic factor. NK cells accumulate in the atherosclerotic plaques, and contribute to the expansion of necrotic cores and atherosclerotic lesion development through producing perforin and granzyme B [55]. Antibodies-mediated depletion of NK cells in ApoE −/− mice greatly attenuated atherosclerosis [55]. Another study has also proved the pro-atherogenic role of NK cells by showing that depletion of NK function decreases atherosclerosis in Ldlr −/− mice [56]. Here, bilirubin-induced decrease of NK cells, especially Ly6C + NK cells, was observed. Due to the proatherogenic functions of NK cells, the reduction of NK cells after bilirubin treatment may benefit the remission of atherosclerosis.
MDSCs comprise immature myeloid cells that consist mainly of progenitors and precursors of granulocytes, macrophages, and dendritic cells [57]. MDSCs primarily mediate immunosuppressive function through inhibiting T cell activation and proliferation, promoting the development of regulatory T cells, blocking antigen recognition, and suppressing NK cell cytotoxicity [57,58]. Based on MDSC phenotypic and morphological features, two major subsets, named polymorphonuclear (PMN) and monocytic (M)-MDSC, were distinguished. In mice, the expression of CD11b, Ly6G, and Ly6C were introduced to identify these two MDSC subpopulations: PMN-MDSC (CD11b + Ly6G + Ly6C low ) and monocytic (M)-MDSC (CD11b + Ly6G − Ly6C high ) [59]. MDSCs have been extensively investigated in tumor immunology, and accumulating evidence has demonstrated their involvement in obesity and atherosclerosis [60,61]. By suppressing T cells, MDSCs are able to reduce atherosclerotic lesion development in Ldlr −/− mice [61]. In our study, bilirubin significantly attenuated the atherosclerosis and caused a significant increase in all MDSC, especially PMN-MDSCs. Moreover, MDSCs are also significantly and negatively correlated with LDL and TCHO, two important markers of atherosclerosis. Thus, bilirubin-induced MDSC increase may also improve atherosclerosis.
Overall, our data clearly show the lipid-lowering and cardioprotective properties of bilirubin in atherosclerosis and that the detailed mechanisms involve inhibition of cholesterol synthesis through promoting HMGCR degradation and regulation of the immune systems.
Notably, mass cytometry-based single cell analysis was introduced to comprehensively dissect the regulatory effects of bilirubin on different immune cell lineages and T cell subsets. Further, these comprehensive data suggest that bilirubin negatively regulates atherosclerosis not only by affecting cholesterol synthesis but also by reshaping the immune system. In addition to the methodological novelty of our work, our results directly propose bilirubin as a novel therapeutic strategy for improving atherosclerosis treatment.
Additional file 1: Figure S1. (A) ApoE −/− mice fed with a westerntype diet were intraperitoneally injected with or without bilirubin. After 10 weeks, mice were scarified, and the body and liver weight were measured. The relative liver weight was expressed as liver weight to body weight ratio. (B) Concentrations of liver enzymes were determined by biochemical analysis. Error bars represent mean ± SD. n (control) =8, n (bilirubin) =14. *P < 0.05. Differences were determined by an unpaired t-test. Figure S2. (A) Quantification of hepatic HMGCR protein expression in ApoE −/− mice treated with bilirubin (n=3) or vehicle (control, n=3). (B) BALB/c mice were received intraperitoneal injection of bilirubin (20 mg/ Kg/day, n=4) or vehicle (n=4) once a day for 3 days and sacrificed at day 4. The expression levels of HMGCR in the livers were determined by western blot. Bar plots showing the relative HMGCR protein expression in the livers (right panel). *P < 0.05. Figure S3. After determining the expression of 27 surface markers on the spleen cells using mass cytometry, the T cells expressing CD3 were gated and analyzed using viSNE. viSNE was used to visualize the distribution of the spleen T cells from both groups. Cells on the viSNE map were colored by normalized expression of indicated surface markers. (B) Bar plots showing the frequencies of T1 or T2 in spleen T cells obtained from ApoE −/− mice treated with or without bilirubin. Error bars represent mean ± SD. n (control) =8, n (bilirubin) =14. Figure S4.  Figure S5. After determining the expression of 27 surface markers on the PWB cells using mass cytometry, the T cells expressing CD3 were gated and analyzed using viSNE. (A) viSNE was used to visualize the distribution of the peripheral blood T cells from both groups. Cells on the viSNE map were colored by normalized expression of indicated surface markers. (B) Bar plots showing the frequencies of T1 or T2 in peripheral blood T cells obtained from ApoE −/− mice treated with or without bilirubin. (C) Bar plots showing the frequencies of T28 or T29 in peripheral blood T cells obtained from ApoE −/− mice treated with or without bilirubin (left panel). Dot plots showing the Pearson correlation coefficients for relationships between the concentrations of total bilirubin in peripheral blood and frequencies of T28 or T29 in peripheral blood T cells. (D) Dot plots (n=22) showing the Pearson correlation coefficients for relationships between the concentrations of TCHO in peripheral blood and the frequencies of T7, T9, T13 in peripheral blood T cells. Error bars represent mean ± SD. n (control) =8, n (bilirubin) =14. Differences were determined by an unpaired t-test. Correlations were determined by a Pearson test. Table S1. Mass cytometry antibody reagents.