A single-nucleotide mutation within the TBX3 enhancer increased body size in Chinese horses

Summary Chinese ponies are endemic to the mountainous areas of southwestern China and were first reported in the archaeological record at the Royal Tomb of Zhongshan King, Mancheng, dated to approximately ∼2,100 YBP.1 Previous work has started uncovering the genetic basis of size variation in western ponies and horses, revealing a limited number of loci, including HMGA2,2LCORL/NCAPG,3ZFAT, and LASP1.4,5 Whether the same genetic pathways also drive the small body size of Chinese ponies, which show striking anatomical differences to Shetland ponies,6 remains unclear.2,7 To test this, we combined whole-genome sequences of 187 horses across China. Statistical analyses revealed top association between genetic variation at the T-box transcription factor 3 (TBX3) and the body size. Fine-scale analysis across an extended population of 189 ponies and 574 horses narrowed down the association to one A/G SNP at an enhancer region upstream of the TBX3 (ECA8:20,644,555, p = 2.34e−39). Luciferase assays confirmed the single-nucleotide G mutation upregulating TBX3 expression, and enhancer-knockout mice exhibited shorter limbs than wild-type littermates (p < 0.01). Re-analysis of ancient DNA data showed that the G allele, which is most frequent in modern horses, first occurred some ∼2,300 years ago and rose in frequency since. This supports selection for larger size in Asia from approximately the beginning of the Chinese Empire. Overall, this study characterized the causal regulatory mutation underlying small body size in Chinese ponies and revealed size as one of the main selection targets of past Chinese breeders.

In brief Liu et al. discover a regulatory mutation in TBX3 driving size variation in Chinese horses. Reporter assays and mice models confirm that it can increase TBX3 transcription and the limb length. They find this mutation first occurred 2,300 years ago and rose in frequency since, revealing size as one of the main selection targets of past Chinese breeders.

Genome sequencing and phylogenetic analysis
We sequenced the genomes of 48 ponies and nine horses at an average depth of coverage of $11.13 and combined them with 130 genomes previously characterized 7,8 (Data S1A; STAR Methods). Our dataset includes 52 Chinese ponies from four breeds, 100 Chinese native horses from 12 breeds, 23 imported Falabella ponies, seven Thoroughbred horses, and five Przewalski's horses. Stringent quality filtering identified a total of 20,428,318 SNPs, 65.1% of which were within introns and 20.2% within intergenic regions (Data S1B; STAR Methods). Principal component analysis (PCA) and phylogenetic tree reconstruction revealed that all Chinese ponies clustered together in a single, monophyletic group (hereafter referred to as Southern ponies, SC; Figures S1A and S1B). Two other main clusters consisting of Northern horses (NC) and Qinghai-Tibetan horses (QT) were identified, in line with previous analyses. 8 Pony breeds showed genome-wide heterozygosity levels generally higher than those found in horse breeds (Data S1C; STAR Methods). Their genomes were also characterized by shorter runs of homozygosity fragments (ROHs) ( Figure S1D), possibly reflecting relaxed selective pressures in Chinese ponies relative to horses.

Selection signatures among horses and ponies
We scanned the genomes of 152 Chinese native ponies and horses for signatures of positive selection that may underpin their different statures. To achieve this, we calculated three complementary statistics along the horse reference genome using 100-kb-long sliding windows (step size = 15 kb). The first statistic was the pairwise genetic differentiation F ST index for identifying those genomic regions maximizing allelic frequency differences between ponies and horses. The second statistic measured relative differences in nucleotide diversity between the two groups (qp(pony/horse)), with high or (low) qp values indicative of positive selection in Chinese horses or ponies, respectively. Selection candidates were also identified through the calculation of crosspopulation extended haplotype homozygosity (XP-EHH). In order to limit false positive identifications, we considered the top 1% windows from all three scans, which provided a total of 153 windows encompassing 4.66 Mb and 69 protein-coding genes (Figure 1A; Data S2A-S2D; STAR Methods). These selection candidates showed limited genetic diversity in horses, but not in ponies, suggesting positive selection for horse-related phenotypes ( Figure 1A). Most of the protein-coding genes were significantly enriched for functional Gene Ontology categories related to the anterior/posterior pattern specification (p = 7.77eÀ7, following g:SCS ''Set Counts and Sizes'' correction 9 ), the embryonic skeletal system development (p = 3.07eÀ6), embryonic skeletal system morphogenesis (p = 1.27eÀ5), and pattern specification process (adjusted p = 3.59eÀ5) (Data S2E). Importantly, all the four outlier genomic regions providing the strongest signatures of selection contained TBX3 (ECA8:20,640,001-20,740,  000 bp, ECA8:20,655,001-20,755,000 bp, ECA8:20,625,001-20,725,000 bp, and ECA8:20,610,001-20,710,000 bp) ( Figure 1A; Data S2D) and were further refined into an 80-kb region (ECA8: 20,640,000-20,720,000 bp) repeating calculations in 20-kb sliding windows ( Figure 1B; Data S2F-S2H).
TBX3 is a T-box gene involved in the development of the anterior/posterior axis, 10 which was previously reported to undergo selection in Debao ponies. 2 This indicated that TBX3 could play an important role in driving body size differences between Chinese ponies and horses. Interestingly, the strongest selection signature identified in Falabella ponies, which represent the smallest miniature horse breed on the planet, 11 was located around another gene, HMGA2 ( Figure S2A), which was previously reported to be associated with short stature in Shetland ponies. 12 Including other western pony breeds to the selection analysis confirmed previously reported selection candidates, such as FGF12, DLX3, and PDK2, 13,14 but not TBX3 (Data S2I-S2K). This suggested a potentially different genetic basis for size-related phenotypes in Chinese and western ponies (Data S2J and S2K).
In addition to showing selection signatures in Debao ponies, 2 TBX3 was previously described to cause asymmetric deposition of skin pigments characteristic of Dun phenotypes when present in the ancestral, non-deleted allelic state. 15 In order to check whether the selection signature identified above could in fact target Dun phenotypes, we screened the presence of not deleted (Dun and non-dun1) and deleted (non-dun2) alleles in all 187 genomes present in our panel ( Figure S2B; Data S2L). We found Dun alleles in five Przewalski's horses (PrZ) and one Chinese Hequ horse, but not in any of the remaining 99 Chinese horses, 52 Chinese ponies, 23 Falabella ponies, and 7 Thoroughbred horses. This ruled out Dun phenotypes driving the selection signatures identified.   Report Genome-wide association study We next carried out a genome-wide association analysis for the wither height of all 152 Chinese samples (STAR Methods). This analysis accounted for the underlying population structure by adjustment of above mentioned PCA covariates and genetic relatedness among all pairs of individuals (STAR Methods). It revealed two statistically significant regions (p < 6eÀ8), encompassing 82 SNPs ( Figures 1C and 1D; Data S2M). The first region spanned the TBX3/TBX5 locus (ECA8:20,566,781-20,865,109), which also provided the top-selection candidate. The second region comprised the HMGA2 locus ECA6:82,576,739-82, 701,867. This result indicated that TBX3, TBX5, and/or HMGA2 may provide the genetic basis for the short stature of Chinese ponies.
In order to further identify the main underlying genetic driver, we investigated patterns of differential placental gene expression in four Chinese YiLi horses and three NiQi ponies (Figure S2C). We focused on the placenta as this tissue is essential for fetal growth, 16 and the placental surface structure is known to be associated with the neonatal body size. 17 Overall, we identified a total of 620 differentially expressed genes in the transcriptomes of YiLi horses and NiQi ponies (Data S2N). Interestingly, TBX3 was highly expressed in horse placenta. In contrast, none of the other two candidate genes (TBX5 and HMGA2) showed substantial expression levels (fragments per kilobase of exon per million fragments mapped, FPKM < 0.3). TBX3 placental expression was also 4.3-fold higher in YiLi horses than NiQi ponies (p = 0.04; Figure 1E), indicating positive association of TBX3 expression levels and adult wither height. This finding was in line with the absence of missense mutations in the TBX3 region examined, and the lower nucleotide diversity and negative Tajima's D statistics found in this region in horses ( Figure 1B). Altogether, this suggested that the genetic variation present within the TBX3 region could regulate gene expression and potentially drive size differences in Chinese horses and ponies. We thus decided to focus the following experimental validation efforts on the TBX3 locus.
Identifying TBX3 causative variant As the TBX3 region showing selection signatures and associated with wither height encompassed approximately $80 kb, we carried out a number of analyses to narrow down the list of possible causative variants (STAR Methods). We first explored the patterns of linkage disequilibrium (LD) in the region using a haplotype block analysis (Data S2O), which revealed nine blocks (H1-H9, haplotype p < 1eÀ5 and #SNPs > 1) containing 16 candidate SNPs ( Figure 1F). Of these, only the H2 block, which contained two SNPs (TBX3-EN1, ECA8:20,644,525 and TBX3-EN2, ECA8:20,644,555), overlapped with a 243-bp region showing extremely elevated conservation GERP scores in 90 mammals (score = 504.5, p = 1.51eÀ186). These two SNPs also overlapped an enhancer region previously identified in large-scale mouse embryonic studies 18 ( Figure 1F). Together, this suggested that the TBX3-EN1 and TBX3-EN2 SNPs may play a strong role in the regulation of TBX3 expression ( Figure 1F; Data S2P). We next measured the allelic frequency at the TBX3-EN1 and TBX3-EN2 SNPs in a comprehensive panel of 189 Chinese ponies and 574 Chinese horses using kompetitive allele specific PCR (KASP) assays (Data S3A-S3C; STAR Methods). Chisquare tests confirmed strong statistical association with the horse wither height, which was stronger for TBX3-EN2 (p = 2.34eÀ39) than TBX3-EN1 (p = 8.18eÀ29) ( Figure 2). Our association model suggested that the variation at TBX3-EN2 (TBX3-EN1) could explain as much as $20.3% ($15.1%) of the height variation ($10 cm) in this panel ( Figure 2A).
Interestingly, the A allele at TBX3-EN2 was most common in seven pony breeds (frequency = 65.6%, N = 189) and the Przewalski's horses (frequency = 100%, N = 13). The G allele, however, reached a higher frequency in 18 horse breeds (frequency = 81.2%, N = 561). The frequency of the G allele was most limited (21.1%) in Debao ponies (wither height < 100 cm) but the highest in YiLi horses (85.7%, wither height > 150 cm) and Northern horse  . This confirmed a model in which the A variant could represent the ancestral allele at TBX3-EN2 and the derived G allele increased in frequency among Chinese horses in response to selection for higher wither height. The distribution of the genetic variation at TBX3-EN1 in our panel suggested a possibly similar scenario for the T allele.
We further reconstructed the frequency trajectories for the TBX3-EN1 and TBX3-EN2 derived alleles during the last $6,000 years, leveraging an extensive ancient genome time series previously released for horses. 19 We found that most ancient samples carried the A allele at TBX3-EN2 until $2,300 YBP, a time after which the G allele increased in frequency steadily until the present time ( Figures 3A and 3B; STAR Methods). Importantly, post-mortem DNA decay can lead to the accumulation of G/A substitutions in ancient DNA data. 20 The allelic trajectory reconstructed here was, however, not correlated to postmortem DNA damage levels, as measured through mapDam-age2 deamination parameters 21 (Figures 3B and S3). Therefore, our analyses supported ongoing selection for the G allele at TBX3-EN2 in relation with increased wither height from the late first millennium BCE. The T allele present at the TBX3-EN1 also strongly associated with horse wither height and was randomly distributed at the beginning but followed a similar trajectory with the increase of the TBX3-EN2 G frequency, in line with the selection acting on linked genetic variants ( Figure 3B). Moreover, the genotyping profiles of 245 ancient horses 19 within the 10-kb region flanking the TBX3 enhancer containing TBX3-EN1 and TBX3-EN2 were closer to those found in Chinese ponies (N = 52) than in Chinese horses (N = 100) ( Figure 3C). Modern Chinese horses were characterized by limited sequence variation, in line with the positive selection acting in the region.

Functional validation of the TBX3 causative variant
The analyses based on ancient and modern genomic data supported a model in which the selection for the T mutation at TBX3-EN1 and G mutation at TBX3-EN2 may have driven increased sizes in Chinese horses over the last 2,300 years. Reporter assays in HEK293T cells, however, indicated no significant difference in expression levels between the vector carrying both derived alleles at TBX3-EN1 and TBX3-EN2 (T and G), and that carrying the derived allele at TBX3-EN2 only (C and G; Figure 4A). This contrasted with the significant drop detected for vectors carrying the ancestral A allele at TBX3-EN2, regardless of the allele associated at TBX3-EN1 (T and A or C and A). Together, this indicated that the G mutation at TBX3-EN2 alone was sufficient to increase TBX3 expression levels and, thus, likely represented the causative allele. We then speculated that the G allele at TBX3-EN2 could affect the TBX3 enhancer region and TBX3 expression and may have driven larger wither height in Chinese horses during the last $2,300 YBP. Interestingly, the horses examined for placental transcriptome expression supported such a model, as all four YiLi horses were homozygous G at TBX3-EN2 while all three NiQiang ponies were homozygous A at TBX3-EN2 ( Figure 1E).
In order to further validate the functional consequences of TBX3 mutations, we constructed knockout (KO) mice for the TBX3 enhancer using the CRISPR/Cas9 technology (STAR Methods), in which a core region showing high sequence similarity among vertebrates and including the TBX3-EN2 SNP was deleted (Data S4A; Figure 4B). Such TBX3 enhancer À/À mice showed a 6-fold reduction in TBX3 forelimb transcription relative to the WT (p = 0.004). Similar reduction in TBX3 transcription levels was also detected in embryos and 2-week leg bones (p < 0.05) (Figures 4C and S4A). No differences in TBX3 transcription were detected in muscle, skin, or heart tissues. Our results thus supported that the 35-bp region acts as a limb-specific enhancer for TBX3 expression, in line with dual-luciferase reporter experiments in the HEK293T cell lines ( Figure 4A).
Finally, we used micro-CT and X-ray imaging techniques on 11-and 2-week mice to measure the length of all bone structures involved in limb formation (i.e., the humerus, radius, metacarpal, and fingers for forelimbs and the femur, tibia, metatarsal, and toes for hindlimbs) ( Figure S4B). The 11-week mice included the littermate of six males (WT:KO = 3:3) and six females (WT:KO = 3:3), while 2-week mice include the littermate of ten males (WT:KO = 3:7) and eight females (WT:KO = 5:3). KO and WT mice showed statistically significant differences in the length of all measured limb bones, with size reduction in KO mice ranging from 1.06% to 21.42% (p < 0.05). This finding was consistent across sexes and the two developmental stages investigated (Data S4B-S4D; Figures 4D, 4E, and S4C). Interestingly, the size reduction was more pronounced in forelimbs than hindlimbs, especially in 2-week female mice, reaching almost 20% for the third and fourth metacarpal digits (Figures 4E and  S4C). This was in agreement with previous studies reporting forelimb malformation in TBX3 mice mutants. 22,23 In this study, we used whole-genome sequence data to identify candidate regions underlying the striking differences in wither height found among Chinese horse and pony breeds. One such region encompassed the TBX3 gene and showed strong signatures of positive selection, indicating possible functional relevance. We further restricted and functionally validated genotype-phenotype association for one single A/G nucleotide variant located within the TBX3 enhancer (ECA8:20,644,555), with the G allele predominantly found in horses with increased TBX3 transcription in cell reporter assays, the placenta, and during limb development. Mice models deleted for the 35-bp region carrying the mutation within the TBX3 enhancer also showed reduced limb sizes at 2-and 11-week developmental stages. Altogether, our work demonstrated that TBX3 expression played a crucial role in driving body size variation in Chinese horses and ponies. Interestingly, the A allele appeared ubiquitous among (E) Metacarpal length of TBX3 enhancer in KO À/À and WT mice at the age of 2 weeks; the x axis 1-5 means the five digits (refer to D) and sex separately.

ll
OPEN ACCESS ancient horses predating the Iron Age in the most extensive ancient genome time series available for horses. 19 The G allele was first detected among Berel' Pazyryk horses some $2,300 YBP in the Altay region 24 and steadily increased to presentday frequencies thereafter. This supported a scenario in which past breeders started to select horses of larger sizes concurrently with the beginning of the Chinese Empire during the Qin dynasty. 25,26 Importantly, all Przewalski's horses investigated carried the A ancestral allele, which was also predominant among Chinese ponies. Przewalski's horses and Chinese pony breeds were, thus, not subject to selection for size during their recent evolutionary history. As the G allele, most frequent in horses, is also dominant in western ponies and miniature horses, including Falabellas, our work unveiled different mechanisms underlying size shifts in the history of horse breeding. It also offered a new genetic marker to the horse breeding industry for size improvement and/or pony conservation.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:

DECLARATION OF INTERESTS
The authors have a patent related to this work. Report RESOURCE AVAILABILITY

Lead contact
Further information and requests for resources, material and reagents should be addressed and will be fulfilled by the lead contact, Lin Jiang (jianglin@caas.cn).

Materials availability
All plasmids and mouse lines generated in this study are available upon request.

Data and code availability
The NCBI accession number for the genomic data reported in this paper is PRJNA780752. All other previously published genomic data used in this study are available in the key resources table. This paper does not report original code. Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

EXPERIMENTAL MODEL AND SUBJECT DETAILS
Candidate SNP genotyping in a larger population We used KASP (Kompetitive Allele-Specific PCR) 43 as a cost-effective approach to genotype a panel of 763 horses and ponies originating from 26 Chinese populations at the for two TBX3 enhancer positions (TBX3-EN1 and TBX3-EN2; Data S3A and S3B). The allele-specific primers were designed by carrying the standard FAM and HEX tails and with the targeted SNP at the 3 0 end (Data S3C). PCR cycling was performed using the following protocol: 1) hot start, 95 C, 15 min; 2) ten touchdown cycles, 95 C, 20 s; 3) touchdown at 65 C initially and decreasing by À1 C per cycle for 25 s; and 4) annealing, 95 C, 10 s, 57 C, 60 s, 30 cycles. The genotyping of KASP markers was conducted at the China Golden Marker (Beijing) Biotech.

Dual-luciferase expression assays
Dual-luciferase assays were carried out in HEK293T (Human embryonic kidney 293 cells). HEK293T cell line was purchased from Beijing university medical college hospital. HEK293T cells were cultivated in Dulbecco's Modified Eagle Medium DMEM basal medium (Glibco, USA), supplemented with 10% heat-inactivated fetal bovine serum and penicillin (0.2 U/mL)/streptomycin (0.2 mg/mL)/ L-glutamine (0.2 mg/mL) (GIBCO). The 63-bp TBX3 enhancer region carrying both TBX3-EN1 and TBX3-EN2 mutations was isolated by PCR with the primers listed in Data S3D, following purification and digested with KpnI and NheI restriction enzymes before being inserted into PGL4.23-Basic vectors (EasyGeno, QIAGEN). A total of four vectors displaying different combinations of mutations of TBX3-EN1 and TBX3-EN2 were constructed and tested with reporter assays for regulatory expression effects. These included: the horse derived mutation (T and G), the TBX3-EN2 single mutation (C and G), TBX3-EN1 single mutation (T and A), and the double TBX3-EN1 and TBX3-EN2 ancestral allele (C and A), where the first letter provides the TBX3-EN1 allele, and the second letter the TBX3-EN2 allele. The TBX3 enhancer region was constructed into dual luciferase reporter vector pGL4.74 and transfected into HEK293T cell by Lipofectamine3000 (ThermoFisher Scientific, USA), following the procedure provided by the manufacturer.


Linear Six-Carbon Sugar Alcohols Induce Lysis of Microcystis aeruginosa NIES-298 Cells

Cyanobacterial blooms are a global concern due to their adverse effects on water quality and human health. Therefore, we examined the effects of various compounds on Microcystis aeruginosa growth. We found that Microcystis aeruginosa NIES-298 cells were lysed rapidly by linear six-carbon sugar alcohols including mannitol, galactitol, iditol, fucitol, and sorbitol, but not by other sugar alcohols. Microscopic observations revealed that mannitol treatment induced crumpled inner membrane, an increase in periplasmic space, uneven cell surface with outer membrane vesicles, disruption of membrane structures, release of intracellular matter including chlorophylls, and eventual cell lysis in strain NIES-298, which differed from the previously proposed cell death modes. Mannitol metabolism, antioxidant-mediated protection of mannitol-induced cell lysis by, and caspase-3 induction in strain NIES-298 were not observed, suggesting that mannitol may not cause organic matter accumulation, oxidative stress, and programmed cell death in M. aeruginosa. No significant transcriptional expression was induced in strain NIES-298 by mannitol treatment, indicating that cell lysis is not induced through transcriptional responses. Mannitol-induced cell lysis may be specific to strain NIES-298 and target a specific component of strain NIES-298. This study will provide a basis for controlling M. aeruginosa growth specifically by non-toxic substances.

INTRODUCTION
Excessive proliferation of harmful cyanobacteria, referred to as cyanobacterial blooms or algal blooms, which mainly occurs with high nutrients (nitrogen and phosphorus), high temperatures, and high residence time, causes serious water quality deterioration such as scum formation, unpleasant odor and taste, hypoxia, and mass fish death in aquatic ecosystems (Codd, 2000;Jankowiak et al., 2019). Further, many harmful cyanobacteria produce toxic compounds such as microcystin and nodularin that can be lethal to animals and humans (Falconer, 1996;van Apeldoorn et al., 2007). Therefore, cyanobacterial blooms are a serious global concern, and their prevention and control are of increasing importance for protecting freshwater sources as well as animal and human health Zhao et al., 2019). Various cyanobacteria including Microcystis, Anabaena, Aphanizomenon, Lyngbya, Cylindrospermopsis, Nodularia, Oscillatoria, Trichodesmium, and Planktothrix are known to be associated with cyanobacterial blooms (Kruskopf and Du Plessis, 2006;Westberry and Siegel, 2006;Briand et al., 2008;Ploug, 2008;Sharp et al., 2009;Al-Tebrineh et al., 2012). Among these, Microcystis aeruginosa is a predominant and infamous cyanobacterium in most cyanobacterial blooms and produces microcystin, which is known to cause liver toxicity and neurotoxicity in humans (Ozawa et al., 2005;Hu et al., 2016).
One of the ways to prevent cyanobacterial bloom occurrence is to reduce the nutrient contents, the main cause of cyanobacterial blooms, in the water body (Xu et al., 2015). Various physical, chemical, and biological approaches have been used to control or suppress existing cyanobacterial blooms (Lawton et al., 1998;Wu et al., 2012;Fan et al., 2013;Visser et al., 2016;Wu et al., 2016;Yu et al., 2017). Physical methods such as mechanical removal, mixing, and ultrasound are difficult to use for large water bodies and consume large amounts of energy (Lawton et al., 1998;Wu et al., 2012;Visser et al., 2016). Chemical agents such as copper algaecides, ferric salt flocculants, clay flocculants, and glyphosate have been relatively widely used for controlling cyanobacterial blooms, but can be harmful to other living things and result in secondary pollution (Fan et al., 2013;Wu et al., 2016;Yu et al., 2017). As alternatives to physicochemical approaches, biological methods using antagonistic organisms or microorganism-or plantderived natural compounds specifically targeting cyanobacteria, particularly M. aeruginosa, have attracted much attraction (Van Wichelen et al., 2016;Lee et al., 2020). For instance, Lee et al. (2020) reported that amentoflavone derived from the plant, Selaginella tamariscina, showed a selective killing effect on M. aeruginosa.
During the investigation of interactions between M. aeruginosa and symbiotic heterotrophic bacteria, a mutation in a carbon metabolism-related-genes of a symbiotic bacterium showed a significant effect on M. aeruginosa growth. We hypothesized that organic carbons might have different effects on M. aeruginosa growth and we tested the effects of certain organic compounds on the growth of M. aeruginosa NIES-298. Unexpectedly, we found that mannitol, considered a beneficial compound with antioxidant or cell-protective effects in organisms, induced rapid cell lysis in M. aeruginosa NIES-298 specifically (Smirnoff and Cumbes, 1989;Dianawati et al., 2012). Therefore, in this study we investigated the mannitol-induced cell lysis in M. aeruginosa using various tests and microscopic observations. These results can provide some insights in M. aeruginosa control. ·s −1 with a 12 h/12 h light/dark cycle and were maintained by transferring 2% (v/v) of cultures to fresh BG-11 medium every 3 weeks. All experiments used the exponentially growing cells, approximately 15 days after their transfer according to the growth curve. All treatments were performed at 3 h of the light cycle to ensure high metabolic activity (Straub et al., 2011).

Effects of Various Organic Compounds on Microcystis aeruginosa Growth
To assess the effects of various organic compounds (glucose and fructose as basic carbon sources; pyruvate, gluconate, and mannitol) and yeast extract as a known death-inducer of Synechococcus (Christie-Oleza et al., 2017) on M. aeruginosa growth, each organic compound was added to 30 ml of strain NIES-298 culture in a 60 ml culture bottle with a filter cap allowing aeration (SPL Lifescience, 70325) to a final concentration of 50 µM [0.1% (w/v) for yeast extract] and the culture bottles were incubated at 25 • C under 40 µmol photons m −2 ·s −1 with a 12 h/12 h light/dark cycle, and the growth of M. aeruginosa was monitored by measuring optical density (OD, 680 nm) using a microplate reader (Synergy TM HTX, Biotek, United States). Because mannitol had a clear growth inhibition effect on strain NIES-298, the effects of sugar alcohols (glycerol, erythritol, inositol, iditol, fucitol, galactitol, sorbitol, mannitol, and maltitol) depending on their chemical structures as well as of different mannitol concentrations (0-100 µM) on M. aeruginosa growth were also assessed using the same procedure.
Because M. aeruginosa dwells with other bacteria in the natural freshwater condition, it was tested if mannitol could induce cell lysis of M. aeruginosa in the presence of other bacteria. Pseudomonas sp. MAE1-K is a symbiotic bacteria isolated from xenic M. aeruginosa KW culture. MAE1-K was inoculated into M. aeruginosa NIES-298 culture to be 10 7 cells/ml to axenic. The co-cultures of M. aeruginosa NIES-298 and Pseudomonas sp. MAE1-K were treated with 50 µM mannitol and the cell lysis of M. aeruginosa NIES-298 was assessed.
The effect of light intensity on the mannitol-induced cell lysis was assessed by culturing M. aeruginosa NIES-298 treated with 50 µM mannitol under continuous dark, 40, 80, and 160 µmol photons m −2 ·s −1 with a 12 h/12 h light/dark cycle at 25 • C. Chlorophyll a concentration was also measured using an acetone extraction method, for the growth assessment (Inskeep and Bloom, 1985). For measuring the chlorophyll a and organic matter released from cells, 1 ml of M. aeruginosa culture was centrifuged for 10 min at 13,300 × g and the supernatant was filtered using a 0.45 µm PVDF syringe filter (Millipore, United States). The concentrations of chlorophyll a and organic matter in the filtrates were determined via the acetone extraction and absorbance at 254 nm (A 254 ) (Weishaar et al., 2003), respectively. The effects of mannitol on the growth of other cyanobacterial strains (strains PCC 7806, NIBR 18, N452, KCTC AG20470, and FBCC 010003) were tested with 50 µM mannitol as described above and assessed by measuring OD (680 nm).

Measurement of Extracellular Mannitol Concentration
The mannitol concentrations in the supernatants of M. aeruginosa cultures at 0, 2, 12, and 24 h after 50 µM mannitol treatment were measured using an Agilent Infinity 1290 UHPLC equipped with an Agilent Infinity Lab Poroshell 120 HILIC-Z column (2.1 mm ×100 m, 2.7 µm) at 35 • C coupled with an Agilent 6550 QTOF mass spectrometer (LC-QTOF-MS). Raffinose (2 µM) was used as the internal standard to normalize the peak intensity of each sample. Water containing 0.3% ammonium acetate (A) and acetonitrile (B) were used as mobile phases at a flow rate of 0.4 ml/min with the following gradient: 0-1 min 87% B, 1-5 min 50% B, and 5-7 min 87% B. The mass spectrometry was operated with electrospray ionization under the following operation parameters: polarity, negative; gas temperature, 250 • C; nebulizer, 35 psig; and MS range, 20-1,000 m/z.

Measurement of Caspase Activity and Oxidative Stress
Because caspase induction has been considered a representative programmed cell death response due to oxidative stresses (Bhattacharjee and Mishra, 2020), the caspase activity of M. aeruginosa NIES-298 cells was assessed using a caspase-3 activity kit (Biovision, United States), according to the manufacturer's instructions. Briefly, M. aeruginosa cells were treated with 25 and 50 µM mannitol and 1 mM H 2 O 2 (as a positive control) and were incubated under the same culture conditions described above. Cells were lysed using a lysis buffer and incubated with the caspase-3 substrate, DEVD-AFC (7-amino-4-trifluoromethyl coumarin), at 37 • C for 1 h. Fluorescence was measured at 404 nm excitation and 505 nm emission and was normalized to the OD values at 680 nm. To assess oxidative stress generation by mannitol treatment, 30 ml of NIES-298 cell cultures in 60-ml air-filtered culture bottles were treated with superoxide dismutase (SOD, 50 U·ml −1 ), catalase (50 U·ml −1 ), and glutathione (1 mM), along with mannitol (50 µM), and their growth was assessed by measuring OD (680 nm).

Observation of Morphological Changes Using Microscopy
Morphological changes in M. aeruginosa cells by mannitol treatment were observed using differential interference contrast microscopy (DICM; Carl Zeiss Axio Scope.A1, Germany), scanning electron microscopy (SEM; Zeiss Supra 55VP, Germany), and transmission electronic microscopy (TEM; Talos L120C, Czech). For electron microscopic analyses, cells were primarily fixed overnight at 4 • C using 2.5% (w/v) glutaraldehyde, rinsed three times with 0.05 M sodium cacodylate buffer (pH 7.4), and post-fixed with 1% (w/v) osmium tetroxide. The fixed cells were dehydrated in an increasing ethanol series (30 min each in 30, 50, 70, 80, 90, and 100% ethanol). For TEM analysis, the fixed cells were embedded in Spurr's resin, sectioned, stained, and examined by TEM (Spurr, 1969). For SEM analysis, the fixed cells were dried with a critical point dryer using liquid carbon dioxide, sputter-coated with gold, and then examined with SEM (Jeon et al., 2003).

Transcriptomic Analysis
To investigate mannitol-induced lysis of M. aeruginosa NIES-298 cells at the transcriptional level, total RNA was extracted from M. aeruginosa NIES-298 cells treated with 25 and 50 µM mannitol for 10 min using a RNeasy Mini kit (Qiagen, Germany). rRNA was removed from the total RNA and the resulting mRNA was sequenced on an Illumina Hiseq 4000 platform at Macrogen (Korea). Low-quality sequencing reads were removed using Sickle (version 1.33) 1 and high-quality mRNA reads were mapped onto the complete genome of M. aeruginosa NIES-298 (CP046058-9; sequenced by our study group) using Burrows-Wheeler Aligner (BWA) software package (version 0.7.17-r1188) (Li and Durbin, 2010). Gene expression levels were calculated as reads per kilobase of gene per million mapped reads (RPKM). Transcriptomic data was deposited to NCBI Bioproject with the accession number PRJNA691993.

Statistical Analysis
All tests were performed in triplicate. The results were statistically analyzed using unpaired Student's t-test available in R software (version 4.1.0) and were expressed as means ± standard deviations. Statistical significances of the results were calculated compared to the control group at the same incubation time and were indicated with p-values < 0.05 ( * ), < 0.01 ( * * ), or < 0.001 ( * * * ).

Effects of Linear Six-Carbon Sugar Alcohols on the Growth of Microcystis aeruginosa NIES-298
We assessed the effects of various organic compounds on the growth of axenic M. aeruginosa NIES-298 (hereinafter M. aeruginosa) and found that their effects were quite different depending on the organic compounds used (Figure 1). M. aeruginosa cultures supplemented with glucose and pyruvate displayed growth profiles similar to the control cultures (no addition of organic compounds), whereas the growth of cultures supplemented with fructose, gluconate, or yeast extract was shown to be significantly inhibited after 7 days. The growth inhibition of M. aeruginosa by yeast extract was partially consistent with a previous report showing that yeast extract accelerated the death of Synechococcus sp. WH78909 cells (Christie-Oleza et al., 2017).
It has been reported that mannitol is not toxic to many other cyanobacteria, marine and freshwater microalgae (e.g., Synechococcus, Emiliania huxleyi, and Monoraphidium griffithii, ·s −1 with a 12 h/12 h light/dark cycle. Data are indicated as mean values ± standard deviations. Statistical significances depending on organic compounds were calculated compared to the control group (no treatment) and their significant differences were indicated with different colors. * * p < 0.01; * * * p < 0.001. respectively); in contrast, some eukaryotic microalgae can produce mannitol and use it as a storage compound (Obata et al., 2013;Jacobsen and Frigaard, 2014;Yee, 2015). However, in this study, mannitol induced a remarkable decrease in cell density soon after its addition to the cultures. As mannitol is a sugar alcohol with six carbons, we further investigated whether other sugar alcohols could induce M. aeruginosa cell death (Figure 2). Glycerol (3 carbons), erythritol (4 carbons), and inositol (6 carbons, but with a ring structure) showed no growth inhibitory effect on M. aeruginosa. On the contrary, all linear six-carbon sugar alcohols including galactitol, iditol, fucitol, and sorbitol clearly induced M. aeruginosa cell death, similar to mannitol. Maltitol, a 12-carbon sugar alcohol known as 4-O-α-glucopyranosyl-D-sorbitol, had a growth inhibitory effect on M. aeruginosa though its effect was weaker than those of other linear six-carbon sugar alcohols (maltitol). These results suggest that linear six-carbon sugar alcohol moieties (structures) have a detrimental effect on the growth of M. aeruginosa NIES-298. Hereinafter, we performed further studies using mannitol as a representative linear six-carbon sugar alcohol.
To assess the threshold mannitol concentration inhibiting the growth of M. aeruginosa, different concentrations of mannitol were used to treat M. aeruginosa cultures (Figure 3). M. aeruginosa growth was very weakly affected by 10 µM mannitol, whereas concentration-dependent growth inhibition of M. aeruginosa was observed with more than 25 µM mannitol ( Figure 3A), however, the green color of the cultures was maintained for a relatively long time after mannitol treatment (inset of Figure 3B), suggesting that unlike stress-induced death of M. aeruginosa, mannitol does not induce chlorophyll bleaching (Park et al., 2010;Alexova et al., 2011;Yang et al., 2017). The effect of light intensity on the cell lysis of M. aeruginosa NIES-298 by mannitol was also evaluated and the mannitol-induced cell lysis  Figure 1), suggesting that the mannitol-induced cell lysis may proceed through light-independent reactions.

Possible Mechanisms Involved in Mannitol-Induced Microcystis aeruginosa Cell Death
The possible mechanisms by which linear six-carbon sugar alcohols can induce cell death of M. aeruginosa were examined. The first possibility is osmotic stress because cyanobacterial cells may be susceptible to osmotic stress (Shapiguzov et al., 2005). However, because M. aeruginosa was reported to be tolerant even to 7 g/L NaCl (approximately 120 mM), 25 µM mannitol is unlikely to cause osmotic stress enough to induce cell death (Tonk et al., 2007). The second possibility for cell death induction is the accumulation of reactive oxygen species (Zhou et al., 2018). However, treatments of SOD, catalase, and glutathione failed to rescue the mannitol-induced cell death in M. aeruginosa (Supplementary Figure 2), suggesting that oxidative stress may not be the relevant underlying cause. Further, mannitol is known as an antioxidant or compatible solute that protect cells from oxidative and osmotic stresses rather than as an agent inducing oxidative or osmotic stress (Smirnoff and Cumbes, 1989;Sand et al., 2013). The third possibility is the accumulation of organic matter through mannitol metabolism because buildup of organic matter in cultures was reported to be potentially toxic to axenic Synechococcus sp. WH7803 (Christie-Oleza et al., 2017). However, mannitol analysis showed that the mannitol concentration in M. aeruginosa cultures remained unchanged (Supplementary Figure 3). Genome analysis of M. aeruginosa NIES-298 (NZ_CP046058-9) also further indicated that no membrane transport systems and metabolic genes (e.g., mannitol dehydrogenase) were identified for linear six-carbon sugar alcohols. These results suggest that linear six-carbon sugar alcohols may not cause a build-up of organic matter via their metabolism, because M. aeruginosa cannot metabolize and uptake linear six-carbon sugar alcohols.
Recently, programmed cell death has been proposed as an important cell death mechanism in cyanobacteria including Anabaena, Trichodesmium, and Microcystis (Ning et al., 2002;Berman-Frank et al., 2004;Ding et al., 2012). Programmed cell death is initiated by the induction of caspase, which is resulted in various morphological and biochemical responses (Saraste and Pulkki, 2000;Bhattacharjee and Mishra, 2020). Therefore, the caspase activity of mannitol-treated cells was measured; H 2 O 2 , which is known to induce caspase activity by programmed cell death (Zhou et al., 2018), was used as the positive control. Treatment of M. aeruginosa cultures with 1 mM H 2 O 2 significantly induced caspase-3 (-like) activity at approximately fourfold compared to the untreated cells after 48 h (Figure 4). However, no caspase-3 activity was induced by the treatments with 25 and 50 µM mannitol, suggesting that programmed cell death is not the cause of mannitol-induced cell death in M. aeruginosa. Taken together, these results indicate that the cell death mechanism underlying death induction in M. aeruginosa by mannitol may not be explained using the existing cell death mechanisms. Therefore, we investigated the morphological changes induced by mannitol treatment in M. aeruginosa cells using various microscopic observations to further scrutinize the mannitol-induced cell death process.

Morphological Changes Induced by Mannitol Treatment in Microcystis aeruginosa Cells
To investigate the mannitol-induced cell death process, the morphological changes induced by mannitol treatment in The caspase activities were expressed as arbitrary unit (AU). Data are indicated as mean values ± standard deviations. Statistical significances depending on the treatments were calculated compared to the control group (no treatment) and their significant differences were indicated with p-values < 0.001 ( * * * ).
Frontiers in Microbiology | www.frontiersin.org M. aeruginosa cells were observed by DICM, SEM, and TEM ( Figure 5). The DICM images indicated that M. aeruginosa cells without mannitol treatment were in a healthy state showing green color as well as round and regular shapes (Figure 5A). TEM and SEM images also showed that mannitol-untreated cells had round and regular shapes without any cell wall damages ( Figures 5F,L). In particular, mannitol-untreated cells were shown to have intact and normal cell walls, cell membranes, thylakoid membranes, and storage granules in a constant state in the TEM images ( Figure 5F). However, upon mannitol treatment, the cell walls burst rapidly in 10 min ( Figure 5B) and intracellular matter including chlorophyll was released out of cells, while maintaining green color (Figures 5C,D); eventually, the green color was bleached (Figure 5E), indicating the death of M. aeruginosa cells.
The TEM images displayed more marked morphological changes in M. aeruginosa cells upon mannitol treatment. After mannitol treatment, the inner and outer membranes of cell walls were shown to be crumpled with uneven cell shapes, and the periplasmic spaces between inner and outer membranes were increased upon the separation of inner and outer membranes (Figures 5G,H). The thylakoid membranes and granules inside cells became unclear, suggesting the disruption of intracellular structures. Interestingly, many outer membrane vesicles (OMVs) were identified from the outer membranes of cells treated with mannitol (no OMV was detected in mannitol-untreated cells), and eventually, the outer membranes of the cell walls disappeared (Figures 5H-J), probably leading to cell wall disruption and intracellular matter release ( Figure 5K). SEM images taken after 2 h of mannitol treatment showed that the mannitol-treated cells were shrunken, probably due to cell wall disruption and intracellular matter release, and that there were some pores on the cell walls probably leading to the intracellular matter release (Figure 5M). These results suggest that mannitol may induce the increase in periplasmic space and formation of OMV, leading to cell wall disruption and intracellular matter release. H 2 O 2 treatment is also reported to cause M. aeruginosa cell shrinkage, but not to cause the perforation of cell walls (Zhou et al., 2020), suggesting that mannitol and H 2 O 2 induce M. aeruginosa cell death via different processes.
Microscopic observations of M. aeruginosa cells suggest that the disruption of cell walls leading to the release of intracellular matter may occur very rapidly by mannitol treatment. Therefore, the contents of organic matter and chlorophyll a released from M. aeruginosa cells were measured after mannitol treatment (Figure 6). The concentrations of organic matter and chlorophyll a increased quickly and in a concentration-dependent manner by mannitol treatment and these increases were completed in 60 min, suggesting that mannitol-induced cell lysis progresses rapidly to cell wall disruption in M. aeruginosa. Morphological changes and cell membrane disruption in M. aeruginosa have also been reported to be induced by allelopathic biochemicals such as ferulic acid, prodigiosin, and amentoflavone (Matthijs et al., 2016;Wang et al., 2016;Yang et al., 2017;Lee et al., 2020). However, mannitol-induced cell lysis occurred within a short time (minutes to hours), whereas cell lysis by many other allelopathic biochemicals generally occurs over several days.
Based on eukaryotic cell death mechanisms, Zheng et al. (2013) proposed four death modes for cyanobacterial cells, apoptotic-, autophagic-, autolytic, and necrosis-like deaths, to explain the cell death of a cyanobacterial symbiont of Azolla microphylla (Zheng et al., 2013). The cyanobacterial cell death modes were mainly classified based on the structural changes in cell walls, cell membranes, peptidoglycan layers, and intracellular contents, as well as the order of cell death progress (Bhattacharjee and Mishra, 2020). Yang et al. (2017) reported that the death of M. aeruginosa FACHB-1752 cell by prodigiosin proceeded by a combination of necrotic-and apoptotic-like process (Yang et al., 2017). However, our microscopic observations showed that mannitol-induced lysis of M. aeruginosa cells may not be classified into the above-mentioned four cyanobacterial cell death modes and may proceed through a yet unidentified novel cell lysis mode.
Our microscopic observations showed that OMVs that probably led to cell disruption were formed by mannitol treatment. OMV formation, which has been mainly studied in gram-negative bacteria, has been reported from various cyanobacteria, including filamentous freshwater cyanobacteria such as Anabaena (Oliveira et al., 2015) and Cylindrospermopsis (Zarantonello et al., 2018), and marine cyanobacteria such as Prochlorococcus, Synechococcus (Biller et al., 2014), and Synechocystis (Pardo et al., 2015). Generally, OMV formation is known to increase under stress conditions (Oliveira et al., 2016;Zarantonello et al., 2018). It is also reported that Prochlorococcus and Synechococcus may use OMVs as a decoy for a protection from phage infection (Biller et al., 2014) and that E. coli may use OMV to excrete misfolded proteins (McBroom and Kuehn, 2007). To the best of our knowledge, this is the first report regarding OMV formation in M. aeruginosa. It may be difficult to describe the role of OMV in M. aeruginosa; however, our microscopic observations indicate that the mannitolinduced lysis of M. aeruginosa cells may be associated with OMV formation.

Transcriptomic Analysis and Strain Specificity of Mannitol-Induced Cell Death
Because mannitol-induced lysis of M. aeruginosa cells differs significantly from the previously proposed cell death mechanisms in M. aeruginosa, a comparative transcriptomic analysis was performed to determine the mannitol-induced mechanism of M. aeruginosa cell lysis. Only eight and ten genes exhibited different expression with more than twofold changes in RPKM values in 25 µM and 50 µM mannitol-treated cells compared to the control (untreated) cells (Supplementary Table 1). In particular, only four and eight genes were upregulated with more than twofold changes in 25 µM and 50 µM mannitol-treated cells, respectively. Proteins such as hypothetical protein, restriction endonuclease, and 3-oxoacyl-ACP synthase that may not be associated with general cell death processes were differentially expressed, whereas SOD genes (D3800_00390 and D3800_14880) were not differentially expressed. These results contrasted markedly with the previous transcriptomic analyses indicating that hundreds of genes related to oxidative stress defense, electron transfer, photosystem, carbon fixation, central metabolic pathways, and microcystin biosynthesis were differentially expressed in M. aeruginosa under stress conditions such as oxidative stress and toxic chemicals (Harke et al., 2017;Jiang et al., 2020;Yang et al., 2020;Zhou et al., 2020). Four putative caspase genes (D3800_09480, D3800_10800, D3800_10885, and D3800_18795) were identified from the genome of M. aeruginosa NIES-298 (Klemenčič et al., 2015) and these were not differentially expressed under the mannitol-treated conditions (Supplementary Table 2), which was consistent with our result that caspase activity was not induced by mannitol treatment (Figure 4). These transcriptomic results suggest that mannitol-induced lysis of M. aeruginosa NIES-298 cells may not proceed through a cell-death cascade associated with the transcriptional responses proposed previously (Zhou et al., 2020).
Because M. aeruginosa co-exists with other microorganisms in the natural environment, the effect of mannitol on M. aeruginosa was tested in the presence of other species. Pseudomonas sp. MAE1-K (CP023641) isolated from a xenic M. aeruginosa culture were inoculated into the axenic M. aeruginosa NIES-298 culture, followed by mannitol treatment. The results showed that although the lysis effect was slightly weakened, the mannitol-induced cell lysis of M. aeruginosa NIES-298 occurred rapidly (Supplementary Figure 4), suggesting that mannitol can induce the cell lysis of M. aeruginosa even in the presence of other microbial cells. Mannitol treatment was also applied to various other cyanobacterial strains, including M. aeruginosa PCC 7806, M. aeruginosa NIBR 18, M. aeruginosa N452, Synechococcus sp. KCTC AG20470, and Anabaena sp. FBCC010003; however, mannitol did not induce cell lysis or death in any of the other tested strains (Supplementary Figure 5), suggesting that mannitolinduced cell lysis is specific to M. aeruginosa NIES-298 and that linear six-carbon sugar alcohols may target specific components, encoded by strain NIES-298, probably through a physicochemical interaction.

SUMMARY
In this study, we found that linear six-carbon sugar alcohols including mannitol, sorbitol, galactitol, iditol, and fucitol induced cell lysis (death) in M. aeruginosa NIES-298. The mannitolinduced lysis of M. aeruginosa cells occurred rapidly and concentration-dependently and was specific to strain NIES-298. Further, mannitol treatment induced an increase in the periplasmic space of M. aeruginosa NIES-298 through separation of its inner and outer membranes and disappearance of the outer membranes probably through OMVs leading to cell wall destruction and release of intracellular matter including chlorophylls. These cell lysis (death) processes were quite different from the previously proposed cell death modes for M. aeruginosa. Further, we demonstrated that linear six-carbon sugar alcohols may induce the cell lysis of M. aeruginosa NIES-298 through structure-based physicochemical interaction with specific components encoded by strain NIES-298, and not through a cell-death cascade process associated with transcriptional responses.
The high specificity of linear six-carbon sugar alcohols including mannitol to the cell lysis of strain NIES-298 will limit the use as a control agent of algal blooms because algal blooms are not occurred by a single cyanobacterial strain. However, the specificity of linear six-carbon sugar alcohols to strain NIES-298 may be very good merit to develop an environment-friendly agent that controls M. aeruginosa specifically without toxicity to other organisms. Therefore, to apply this study to aquatic environments, further research to investigate the relevant cell lysis mechanisms in more detail is necessary. Even though few previous research succeeded mutagenesis of M. aeruginosa, genetic manipulation, random mutagenesis, or long-term laboratory evolution under the sublethal concentration of mannitol may be able to identify the mechanism of mannitol-induced cell lysis. It will provide a basis for developing novel environment-friendly agents or strategies for the specific control of M. aeruginosa, key cyanobacteria that cause algal blooms in freshwater.

DATA AVAILABILITY STATEMENT
RNA-seq data is publicly available from NCBI BioProject under the accession number of PRJNA691993.

AUTHOR CONTRIBUTIONS
JJ, YS, SJ, and CJ designed the research. JJ conducted experiments and performed the main data analysis. YS conducted the main experiments. JB performed LC/MS analysis. SJ and HP conducted data analysis. CJ supervised the study and obtained funding. YS, JJ, and CJ wrote the manuscript. All authors have read and approved the final version of the manuscript.

FUNDING
This work was supported by grants from the National Institute of Biological Resources (NIBR202129201), Republic of Korea.