A 2-bp frameshift deletion at GhDR,which encodes a B-BOX protein that co-segregates with the dwarf-red phenotype in Gossypium hirsutum L.

WANG Xue-feng,SHAO Dong-nan,LlANG Qian,FENG Xiao-kang,ZHU Qian-hao,YANG Yong-lin,LlU Feng,ZHANG Xin-yu,Ll Yan-jun,SUN Jie,XUE Fei

1 Key Laboratory of Oasis Eco-Agriculture,College of Agriculture,Shihezi University,Shihezi 832000,P.R.China

2 CSIRO Agriculture and Food,Canberra,ACT 2601,Australia

3 Cotton Research Institute,Shihezi Academy of Agriculture Sciences,Shihezi 832000,P.R.China

Abstract Plant architecture and leaf color are important factors influencing cotton fiber yield. In this study,based on genetic analysis,stem paraffin sectioning,and phytohormone treatments,we showed that the dwarf-red (DR) cotton mutant is a gibberellin-sensitive mutant caused by a mutation in a single dominant locus,designated GhDR. Using bulked segregant analysis (BSA) and genotyping by target sequencing (GBTS) approaches,we located the causative mutation to a～197-kb genetic interval on chromosome A09 containing 25 annotated genes. Based on gene annotation and expression changes between the mutant and normal plants,GH_A09G2280 was considered to be the best candidate gene responsible for the dwarf and red mutant phenotypes. A 2-nucleotide deletion was found in the coding region of GhDR/GH_A09G2280 in the DR mutant,which caused a frameshift and truncation of GhDR. GhDR is a homolog of Arabidopsis AtBBX24,and encodes a B-box zinc finger protein. The frameshift deletion eliminated the C-terminal nuclear localization domain and the VP domain of GhDR,and altered its subcellular localization. A comparative transcriptome analysis demonstrated downregulation of the key genes involved in gibberellin biosynthesis and the signaling transduction network,as well as upregulation of the genes related to gibberellin degradation and the anthocyanin biosynthetic pathway in the DR mutant.The results of this study revealed the potential molecular basis by which plant architecture and anthocyanin accumulation are regulated in cotton.

Keywords: cotton,BBX,dwarf,anthocyanin,gibberellin

1.lntroduction

Plant architecture is the morphological basis of crop yield.An ideal plant architecture has the potential advantages of increasing cotton planting density,improving photosynthetic efficiency,adapting to mechanized production and increasing yield. In commercial cotton farming,the architecture of cotton plants is modulated through the exogenous application of plant growth regulators (PGRs) (Zhang Jet al.2020). Although multiple applications of PGRs can effectively control plant architecture,chemical sprays can cause environmental pollution and increase production costs,so alternative growth control methods are needed. Breeding cotton varieties with ideal architecture is the most efficient and sustainable approach for controlling cotton architecture(Zhang Xet al.2020).

Plant height is one of the most important traits of plant architecture. In addition to environmental factors,such as photoperiod,temperature,water and fertilizer,plant height development is also regulated by genetic composition and plant hormones such as gibberellic acid(GA),brassinosteroids (BRs) and auxin (IAA). The first Green Revolution in agriculture was characterized by reduced plant height and semi-dwarf breeding,which was achieved by introducing the mutant gene “sd-1” defect in gibberellin synthesis into rice (Spielmeyeret al.2002),and by introducing the gain-of-function mutant gene “Rht-1”involved in gibberellin signal transduction into wheat(Pearceet al.2011;Van De Velde Ket al.2021). As a result,the new crops have a semi-dwarf plant stature and are more resistant to lodging,which makes high density planting possible and has brought the world’s food production to a new level.

Recently,several individual genes involved in the regulation of plant architecture have been characterized in cotton.GhPAG1andGhPAS1regulate cotton plant development and architecture by affecting the levels of endogenous BRs and mediating BR signaling (Yanget al.2014;Wuet al.2021).Gh_A06G1386encodes a cytochrome P450 family protein,ent-kaurenoic acid oxidase (KAO),which is one of the key enzymes involved in gibberellin synthesis (Zhaoet al.2017).GhGA2OX8(Ghir_A05G017390) encodes gibberellin 2-oxidase,which is also involved in the biosynthesis of gibberellin (Li Let al.2021).GhDREB1B,anEXTRDWARFgene,encodes a dehydration response element binding (DREB) transcription factor,and the overexpression ofGhDREB1Bincreases the level of GA2ox oxidase and reduces the accumulation of bioactive GA (Jiet al.2021). Auxin efflux carrier protein (PIN) was reported to play a rate-limiting role in catalyzing the transport of IAA and it participates in the plant height developmental silencing ofGhPIN3(Gh_D01G1471),leading to a taller plant height in cotton(Maet al.2019). Down-regulating the expression of the strigolactones (SLs) signaling geneGhMAX2in cotton results in reduced plant height and shortened internodes(Penget al.2022). Most of the cotton dwarf genes reported thus far are related to GA,BR and IAA. The available genetic resources for cotton plant architecture breeding are still very limited,and the genetic resources related to important plant architecture traits need to be further explored.

Several cotton mutants showing red coloration in the whole plant or certain organs due to altered anthocyanin accumulation have been reported,including red plant(R1) (Liet al.2019),sub-red plant (Rs) (Lianget al.2020),red petal spot (R2) (Abidet al.2022),and dwarf red mutant (Rd) (McMichael 1942).R1(GhPAP1D) andRs(GhPAP1A) are a pair of homoeologs encoding the R2R3-MYB transcription factors,and the overexpression of either gene results in increased anthocyanin accumulation,leading to red coloration (Liet al.2019;Lianget al.2020). In the promoter regions ofGhPAP1DandGhPAP1A,228-and 50-bp tandem repeats were identified,respectively. The additional repeats contain the binding site for the transcription factors that enhance the expression ofGhPAP1DandGhPAP1Aand the accumulation of anthocyanins.R2controls the appearance of the red spot at the base of each petal (Abidet al.2022). The genetic basis of the Rd mutant has not yet been uncovered. Of these four anthocyanin-related cotton mutants,R1,R2 and RS show normal plant height,while Rd is a dwarf.

Previous studies have found certain relationships between plant height and pigment metabolism. For example,the high-tillering and dwarf phenotype of rice was caused by a single nucleotide transition in theHTD12gene that functions in carotenoid biosynthesis and regulates plant height by affecting the level of phytohormone strigolactone (Zhouet al.2021). Three other photomorphogenesis genes involved in the regulation of both anthocyanidins and plant height areCRY(Cryptochrome) (Chatterjeeet al.2006;Sharmaet al.2014),ELONGATEDHYPOCOTYL5(HY5)(Gangappa and Botto 2016),andBBX21(Xuet al.2018).Both photosynthesis and plant architecture are important factors for determining the fiber yield in cotton,but little is known about the relationship between photosynthesis and plant architecture from a genetic perspective.

In this study,we characterized a natural cotton mutant with an extremely dwarfing plant height and red leaves,herein designated dwarf and red (DR). This mutant is a good candidate for investigating the relationship between photosynthesis and plant architecture,as its dwarfism is closely associated with red coloration. We showed that the mutant phenotype of theDRmutant is controlled by a single dominant locus,and identified the candidate gene underlying the mutant phenotype using mapbased cloning. We further characterized the mutant by a combination of microscopic observations of stem sections,a comparison of anthocyanin accumulation between the mutant and normal plants,and comparative transcriptome analysis,to reveal the molecular mechanism underpinning the dwarf and red phenotype.

2.Materials and methods

2.1.Plant materials and mapping population construction

TheDRmutant was kindly provided by the National Cotton Germplasm Medium-Term Bank in China. An F2segregating population,derived from a cross between the normal phenotype Upland cotton (Gossypiumhirsutum)cultivar Xinluzao 74 (X74) and theDRmutant,was developed for the mapping ofGhDR,the gene responsible for theDRmutant phenotype. The F2population was grown in the field (with 528 and 972 individuals in 2020 and 2022,respectively) and the greenhouse (with 316 individuals) to investigate the segregation ratio of the mutant phenotype and to map the corresponding genetic locus.DRand normal F2plants grown in the greenhouse were also used in the transcriptome sequencing experiment. Recombinants selected from the F2population with crossovers between the flanking markers ofGhDRwere fixed in the F3generation (with 204 individuals) to verify the mapping results and to further fine map theGhDRlocus. The phenotypic segregation ratio of the three F2populations was analyzed by aχ2goodness-of-fit test in SPSS.

2.2.Phenotypic and histological observations

The plant height was measured with a ruler. The standard for plant height measurement is the distance from the cotyledon node to the growing point. For stem section histological observations,tender stems of X74 andDRwere fixed in 10% FM fixative solution (ethyl alcohol 90 mL,formalin 5 mL,and glacial acetic acid 5 mL). The fixed stems were dehydrated,embedded in paraffin,and sectioned with a microtome. The sections were then rehydrated for histological analysis. The anthocyanin contents were analyzed by high performance liquid chromatography-mass spectrometry (HPLC-MS) using the extracts obtained from different tissues according to a previously published solid-liquid extraction (SLE) method(Barneset al.2009).

2.3.Treatment of phytohormones

Cotton seedlings were treated with brassinolide (BR)(16 mg L-1),gibberellic acid (GA3) (25 mg L-1),and auxin (IAA) (17.5 mg L-1) at the one true leaf stage. The cotton seedlings used in the treatments were grown in an artificial climate chamber with a 16 h light/8 h dark photoperiod at a temperature of (28±1)°C. Experiments were carried out with three independent biological replicates per treatment (each replicate with eight plants).

2.4.Detection of plant endogenous gibberellins

The leaves and young stems of X74 orDRwere mixed and stored at -80°C until use. The cryopreserved plant material samples were ground into powder with a grinder(MM 400,Retsch) (30 Hz,1 min). For each sample,50 mg of freshly ground powder was taken and 10 μL of 100 ng mL-1internal standard mixed working solution was added.Then solution was extracted with 500 μL of acetonitrile solution (acetonitrile/water at 90:10,v/v) by vortexing for 15 min,then centrifuged at 4°C and 12 000 r min-1for 10 min to collect the supernatant. Aliquots of 10 μL TEA and 10 μL BPTAB were added to the supernatant.After 1 h of reaction at 90°C,the solution was blown dry with nitrogen. The dried material was reconstituted with 100 μL of acetonitrile/water (90:10,v/v) solution,passed through a 0.22-μm filter,and placed in a sample vial for LC-MS/MS analysis according to a published protocol(Denget al.2017).

2.5.Genomic DNA and RNA extraction

Young leaves from the two parents and F2individuals were collected and stored at -80°C. Young leaf tissues were ground into a fine powder in liquid nitrogen using a hybrid grinding machine. Genomic DNA was then extracted using the CTAB method and stored at -20°C.Total RNA was extracted using the EASY-spin Plus Plant RNA Kit (Aidlab,China) according to the manufacturer’s instructions. The RNA integrity,purity and concentration were determined using an ND2000. Suitable RNA samples (i.e.,those with A260/A280ratio=1.8-2.0,A260/A230ratio ＞1.5,and RNA integrity number ＞8) were stored at-80°C for subsequent qRT-PCR analysis.

2.6.BSA-seq,genotype analysis and fine mapping

Two DNA pools prepared from either dwarf/red or normal F2plants were used in the bulked segregant analysis(BSA) sequencing experiment. The dwarf/red pool contained equal amounts of DNA from 30 dwarf and red individuals and the normal pool contained DNA from 30 normal individuals. The DNA was also extracted from the two parents,X74 andDR. The four DNA samples(two parents and two pools) were processed by the NEBNext?Ultra? II DNA Library Prep Kit for Illumina?(GENOSEQ,Wuhan,China) to construct sequencing libraries,which were sequenced using Illumina HiSeq TM PE150 with a sequencing coverage of 30×. The sequenced reads were aligned to the reference genome TM-1_ZJU_v2.1 (ftp://ftp.bioinfo.wsu.edu/species/Gossypium_ hirsutum/ZJU_G.hirsutum_AD1genome_v2.1) using the maximal exact matches(MEM) algorithm of BWA (version 0.7.15-r1140) to obtain the alignment results in SAM format,which was then converted to BAM format using samtools (version 1.3.1). After sorting the reads with SortSam in the Picard tool (version 1.91),the resulting BAM file was used for analyzing the coverage depth and variant calling. The single nucleotide polymorphisms (SNPs)were used to calculate the Δ(SNP-index) value of each mutation site using the QTLseqr (R package) Software(Takunoet al.2013),and the distribution of the Δ(SNPindex) within a 2-Mb sliding window was plotted across the whole genome. The genomic region with a Δ(SNPindex) value over the threshold line was considered as the candidate region ofGhDR.

To narrow down the interval of theGhDRlocus identified based on the BSA experiment,we used the genotyping by target sequencing (GBTS) platform (Guoet al.2019)to obtain the genotyping results of the 528 F2individuals among the different SNP/InDel sites in the target region.A total of 47 primer pairs (Appendix A) were designed to validate the SNP/InDel variants in the candidate region.These primers covered the SNP sites and their flanking sequences. The BSA-seq and genotype analysis were performed by GENOSEQ Co.,Ltd.,Wuhan,China.

2.7.Transcriptome analysis and qRT-PCR

Transcriptome sequencing was conducted for theDRand normal segregants of the F2population using young stem samples. Three biological replicates were used for both dwarf/red plants (namedDR1,DR2,andDR3) and normal plants (named WT1,WT2,and WT3),and each replicate included five plants with uniform phenotypes. Total RNAs were extracted and reverse transcribed into cDNAs to construct the cDNA libraries. The RNA-seq libraries were sequenced with the Illumina HiSeq platform. The resulting data were filtered to obtain clean data,which were then aligned to the reference genome (TM-1_ZJU_v2.1) using the BWA Software (version 0.7.15-r1140) in order to obtain mapped data. Fragments per kilobase of transcript per million fragments mapped (FPKM) values were used to measure the levels of transcripts or gene expression.

For qRT-PCR,total RNA was reverse transcribed into cDNA using an Easy Script One-Step gDNA Removal and cDNA Synthesis Super Mix Kit (TransGen Biotech,Beijing,China). qRT-PCR was carried out on an ABI7500 real-time fluorescent quantitative PCR instrument (Shanghai,China) in a 10-μL volume containing 200 ng of cDNA,10 pmol L-1of each primer,and 6 μL of AceQ qPCR SYBR Green Master Mix (TransGen Biotech,China) according to the manufacturer’s protocol. The qPCR conditions were as follows: primary denaturation at 95°C for 20 s,followed by 40 amplification cycles of 3 s at 95°C,and 30 s at 60°C.Melting curve analysis was performed to ensure the absence of primer-dimer formation and primer specificity.The primers used in the qRT-PCR are listed in Appendix B. Three biological replicates with independently isolated RNAs were used for each sample,and each RT reaction was loaded in triplicate. The relative gene expression levels were calculated using the 2-ΔΔCtmethod (Livak and Schmittgen 2001).

2.8.Gene cloning and multiple sequence alignment

The primers used in gene and promoter amplification were designed using Primer3 (https://www.primer3plus.com/) and are shown in Appendix B. The cDNA and genomic DNA from X74 andDRwere used to amplify the full-length coding sequences ofGhDRand its promoter,respectively. The PCR products were purified using a product purification kit and sequenced by Sangon Biotech(Shanghai,China). The resulting sequences were aligned with the DNAMAN Software.

2.9.Subcellular localization of GhDR

TheGhDRcDNA was amplified fromDRby primers GH_A09G2280/216-GFP-F/R (Appendix B),and cloned into 35S::GFP to generate the 35S::GH_A09G2280/216-GFP fusion construct. Similarly,the cDNA ofGH_A09G2280/238was amplified from X74 by primers GH_A09G2280/238-GFP-F/R and cloned into 35S::GFP to generate the 35S::GH_A09G2280/238-GFP fusion construct. The vector of 35S::GFP (as the control),35S::GH_A09G2280/216-GFP and 35S::GH_A09G2280/238-GFP were individually introduced into tobacco leaves,which were incubated for 2 days in a low-light environment. GFP signals and chlorophyll autofluorescence in the protoplasts were examined under a confocal laser-scanning microscope (Nikon C2-ER,Shanghai,China) at excitation wavelengths of 488 and 640 nm,respectively.

3.Results

3.1.Phenotypic characterization of the DR mutant

TheDRmutant showed reddish coloration after seed germination. Compared with the normal cotton plants of X74,the stature of the mutant was much shorter from the seedling stage to the mature stage due to a less elongated stem and fewer internodes (Fig.1-A-C). In the end,the plant height of X74 was 86.1 cm and the plant height of theDRmutant was 16.2 cm,so the plant height of theDRmutant was only about one-quarter that of the X74 plant (Fig.1-D). In order to determine the reason for the dwarfing,cytological analysis was performed on the stems of X74 andDR. The results of paraffin sections showed that there was no difference in the cell sizes in the transverse sections of X74 andDR(Fig.2-D and E),but the cell size ofDRin the longitudinal section was significantly smaller than that of X74 (Fig.2-B and C).

The stems and leaves ofDRshowed a red phenotype,largely due to high contents of anthocyanins in all theDRtissues analyzed,including hypocotyls,cotyledons and true leaves (Appendix C).

3.2.GA could partially rescue the dwarf phenotype of the DR mutant,but BR and lAA could not

Fig. 1 Phenotypic comparison between the dwarf-red (DR) mutant and Xinluzao 74 (X74). Comparison of the plant heights of X74(left) and DR (right) at the seedling stage (A),bud stage (B),and mature stage (C). D,plant heights of X74 and DR at the mature stage. The numbers on top of the bars are the internode numbers. ＊＊,Student’s t-test,P＜0.01.

Fig. 2 Comparison of stem paraffin sections of Xinluzao 74 (X74) and dwarf-red (DR). A,stems of X74 (left) and DR (right) at the bud stage. Red frames indicate the positions from which the stem samples were collected and used for sectioning. B,longitudinal section of X74 stem. C,longitudinal section of DR stem. D,cross section of X74 stem. E,cross section of DR stem. The areas of the black boxes in B-E are 40 000 μm2. Scale bars are 100 μm.

Fig. 3 Responses of the dwarf-red (DR) mutant to gibberellic acid (GA3) treatment. A,comparison of DR seedlings in response to H2O or GA3 treatments. B,plant height of the DR mutant following H2O (left) or GA3 (right) treatments at the maturation stage. C,measurements of plant height of the DR mutant following H2O (left) or GA3 (right) treatments at the maturation stage. The numbers on the bar graph represent the number of internodes. ＊＊,Student’s t-test,P＜0.01.

Fig. 4 Map-based cloning of GhDR. The GhDR locus was narrowed down to a 197-kb interval that includes 25 genes. GH_A09G2280 is the candidate gene of GhDR. A 2-bp deletion at the 1 245-1 246th positions of GhDR,resulting in truncation at the 3′ terminus of the DR protein (216 amino acids instead of 238 amino acids due to a premature stop codon). X74,Xinluzao 74;DR,dwarf-red.

Fig. 5 The expression levels of GhDR in leaves (A) and stems(B) of X74 and the dwarf-red (DR) mutant as determined by qRT-PCR. X74,Xinluzao 74. Bars mean SD (n=3). ＊＊,Student’s t-test,P＜0.01.

Fig. 6 The expression levels of the expressed genes in the mapping interval. Dwarf-red 1 (DR1) to DR3 are three replicates of the dwarf-red plants from the F2 population;WT1 to WT3 are three replicates of the normal plants from the F2 population.

Fig. 7 Expression profiles of the genes of the gibberellin biosynthesis pathway in dwarf and red (DR) and X74. The square boxes represents the intermediate products of gibberellin metabolism,and those with reduced content are highlighted in green. The oval boxes represent the differentially expressed genes,of which those with high expression are represented in red,and those with low expression are represented in blue. The three biological replicates of DR and WT are shown in the heatmap from left to right.

Fig. 8 The anthocyanin biosynthesis pathway and the expression profiles of the genes differentially expressed in dwarf and red(DR) and X74. The square boxes represent the intermediate products of anthocyanin metabolism,the red oval boxes represent the highly expressed genes,and the dotted boxes represent the three processes of anthocyanin synthesis. The three biological replicates of DR and WT are shown in the heatmap from left to right.

Fig. 9 Schematic diagram of the formation of the dwarf red phenotype under the regulation of BBX. A model showing the BBX24/25 proteins in the regulation of photomorphogenesis,affecting anthocyanin accumulation and cell elongation arrest.Lines with arrows indicate positive regulation,lines with flat ends indicate negative regulation,and dotted boxes indicate failure.

Most dwarf mutants are caused by impaired phytohormone biosynthesis and/or signaling,so we treated theDRmutant with different plant hormones and found that theDRmutant responded obviously to GA3while it was not sensitive to the treatments with either BR or IAA. GA3treatment significantly promoted stem elongation in theDRmutant (Fig.3),although the final plant height was not as high as that of X74. The plant height of theDRmutant remained unchanged after treatment with either BR or IAA (Appendix D). Previously,the hypocotyl elongation of the BR and GA deficient or insensitive mutants was shown to be inhibited under dark conditions (Alabadiet al.2004). In order to confirm whetherDRis a hormone deficient or insensitive mutant,we dark-treated X74 andDR. After 5 days,there was no significant difference in the lengths of their hypocotyls (Appendix E). We tested the contents of different endogenous gibberellins inDRand found that,except for GA5 and GA8,all other eight gibberellins had a lower levels in theDRmutant than in X74 (Appendix F). Among them,the contents of GA1 and GA4,the biologically active gibberellins,were significantly reduced in theDRmutant,in line with the result of a GA rescue experiment. These results demonstrate that theDRmutant is a gibberellin-sensitive mutant due to deficiencies in active GAs,and the dwarfism is caused by a reduced level of endogenous GA content.

3.3.Genetic analysis of the DR mutant

To investigate the genetic basis of theDRmutant,we generated an F2population by crossingDRto X74. We planted three F2populations in the field or greenhouse and found that the red phenotype and the dwarf phenotype showed a phenomenon of linked inheritance. Based on the phenotype of the F2population,we performed a correlation analysis using leaf color and the mean value of plant height corresponding to each leaf color. There was a strong correlation between leaf color and plant height (R2=0.9647 andR2=0.9999 in 2020 and 2022,respectively) within the 95% confidence interval (Appendix G),suggesting that the red and dwarf phenotype of theDRmutant could be controlled by either the same gene or two different genes that are tightly linked. In either case,we could use the red phenotype to investigate the inheritance mode of plant height.

We phenotyped the F2plants at the mature stage by classifying them into three groups: red,intermediate and green. The segregation ratio of the plant color in the 528 F2plants did not meet the ratio of 1:2:1,with fewer red plants than expected,likely resulting from the death of some of the red and dwarf plants because of their compromised photosynthetic capacity. We thus grew two new F2populations derived from the same cross,one grown in the field (with 972 individuals) and the other grown in the greenhouse (with 316 individuals)with controlled conditions,and recorded the plant color phenotype at the seedling stage. The segregation in both F2populations fitted the 1:2:1 ratio (χ21:2:1=0.025 and 0.925 in field and greenhouse,respectively,df=2,P＞0.05),indicating that the dwarf and red phenotype of theDRmutant is conferred by a single dominant locus (Table 1).

Table 1 Segregation ratios for the dwarf-red (DR) mutant and wild-type categories in the F2 populations

Based on visual observation,we noticed that the red color intensity of the F1plants derived fromDR×X74 was weaker than that ofDR,and the height of F1was intermediate between X74 andDR. This was confirmed by measuring plant height and anthocyanin contents in different tissues fromDRmutants,X74 and their F1,suggesting that the dwarf and red phenotype is controlled by an incompletely dominant locus (Appendix C).

3.4.Fine-mapping the GhDR locus

We performed BSA-seq to map and identify the genomic region underlying the dwarf and red mutant traits. Two DNA pools were made using the F2individuals showing the normal (tall and green) or mutant (dwarf and red)phenotype. Each pool contained the DNA from 30 individuals and was sequenced to a 30× depth. At the same time,the two parents were also sequenced with the same sequencing depth.

The variants between the two pools were called,including 897 951 SNPs and 105 184 InDels. Based on the distribution of Δ(SNP-index),we could locate the candidateGhDRlocus to the interval between 75.89 and 82.43 Mb on chromosome A09 (Appendix H).

Based on the results of BSA-seq,we developed 18 markers in the interval containing theGhDRlocus and genotyped the F2population. An analysis of markerphenotype association allowed us to anchor theGhDRlocus in the region defined by markers M11 and M12,with a genetic distance of～0.5 cM (Fig.4). There are 36 annotated genes in that interval (based on TM-1_ZJU_v2.1).To shorten the physical mapping interval,we developed 22 additional markers between M11 and M12 and used them to genotype 40 individuals showing crossover in the interval.To further fine map theGhDRlocus,the recombinants were selected from the F2population with crossovers between the flanking markers ofGhDR,and a total of 204 individuals were fixed in the F3generation,all of which were genotyped by seven markers (M11,M8457,M8917,M5491,M1085,M5196,and M7212) (Fig.4).

Finally,theGhDRlocus was delimited to a 197-kb region flanked by M8457 and M5491 (Fig.4). This region contains 25 genes (based on TM-1_ZJU_v2.1). Among those 25 genes,the only coding sequence related variant between X74 and theDRmutant (a 2-bp deletion inDR) was found in the third exon ofGH_A09G2280(Appendix I). An InDel marker (M8917) was developed based on this variant,and it was found to co-segregate with the mutant trait. Based on these results,we consideredGH_A09G2280to be the best candidate ofGhDR,andGH_A09G2280encodes a protein with 238 amino acids. In theDRmutant,the 2-bp deletion at the 1 245th and 1 246th nucleotide positions caused a frameshift and induced a pre-mature stop codon inGhDR/GH_A09G2280,leading to 3′ terminal truncation ofGhDR. The length of the mutatedGhDR/GH_A09G2280 protein is 216 amino acids (Fig.4).

3.5.Functional annotation and expression profile of GhDR

As further evidence thatGH_A09G2280is the best candidate ofGhDR,of the 25 genes in the mapping interval,GH_A09G2280was annotated as a homolog of Arabidopsis BBX24,a gene reported to be associated with hypocotyl elongation and anthocyanin accumulation(Appendix J).GH_A09G2280is relatively highly expressed in all tissues analyzed,particularly in stem and leaf (Appendices K and L). Based on the qRT-RCR results,the expression levels ofGhDR/GH_A09G2280in both stem and leaf were significantly higher in theDRmutant than in X74 (Figs.5 and 6).

3.6.Cloning of the candidate gene

To determine whether there are any other sequence variations in theGhDRgene between X74 and theDRmutant,we amplified the promoters (2 000 bp upstream of ATG) and genomic sequences ofGhDRfrom bothDRand X74. Sequence alignment showed that the 2-bp deletion in the third exon ofGhDRis the only difference between the two accessions (Appendix M).

GhDR/GH_A09G2280is a homolog of ArabidopsisAT1G06040,which encodes a B-box zinc finger protein,BBX24. We aligned the Arabidopsis and cotton B-box protein sequences of AT1G06040/BBX24,AT1G75540/BBX21,GH_D05G1922/BBX21,GH_A05G1184/BBX21,GH_A09G2280/BBX24 and its homoeolog GH_D09G2215,and the mutatedGhDR(the truncated GH_A09G2280 from theDRmutant). This alignment showed that while the two B-boxes are highly conserved among these proteins,the nuclear localization domain and the VP domain are missing inGhDRdue to its 3′ terminal truncation (Appendix N).

3.7.Subcellular localization of GhDR

Given the absence of the nuclear localization domain in the mutatedGhDR(Appendix N),we predicted the subcellular localization ofGhDRusing the online tool Cell-PLoc 2.0 (http://www.csbio.sjtu.edu.cn/bioinf/Cell-PLoc-2/). The results indicated that while GH_A09G2280/238 was predicted to be localized only in the nucleus,GhDR(GH_A09G2280/216) was predicted to be localized in both nucleus and cytoplasm,which was confirmed by subcellular localization experiments using GFP fusion proteins (Appendix O).

3.8.Transcriptome analysis reveals changes in the biosynthesis of GA and anthocyanins in the DR mutant

To investigate the molecular mechanism involved in the formation of theDRphenotype,we performed RNA-seq using RNAs from the stems of X74 andDRand carried out a network analysis. Of the 25 annotated genes located at the mapping interval,16 were expressed and four showed ＞1.5 fold up-regulation in theDRmutant,includingGhDR/GH_A09G2280;however,of these four genes,GhDR/GH_A09G2280was the one with the highest expression level and the other three had relatively low expression levels (Fig.6;Appendix P).

Several gibberellin biosynthesis pathway genes,such as those encoding ent-kaurene synthase (KS,GH_A05G1635),gibberellin 13-beta-dioxygenase (GA13ox/CYP714C2,GH_A07G2012),and gibberellin 3-betadioxygenase (GA3ox,GH_A13G1807),were downregulated in theDRmutant (Fig.7;Appendix Q). In addition,the two key genes (GA2ox,GH_A13G1786andGH_D01G0359) that encode gibberellin 2-betadioxygenase,which catalyzes the inactivation of several biologically active gibberellins,had higher expression levels in theDRmutant. On the other hand,GH_A06G0641andGH_D06G0611,which were identified as DELLA4,were highly up-regulated (1.66 and 2.94-fold changes,respectively) in theDRmutant,and this result was further confirmed by qRT-PCR (Appendix R).

Given the increased accumulation of anthocyanins in theDRmutant,we compared the changes in expression of the genes involved in the phenylpropanoid pathway.Genes encoding phenylalanine ammonia-lyase (PAL) and 4-coumarate CoA ligase (4CL) were up-regulated inDRleaves. Several major structural genes of the anthocyanin biosynthetic pathway were highly expressed inDR,such as the two CHS genes,GH_A09G0001andGH_A09G0002,that were up-regulated by 60.7-and 25.7-fold inDR,respectively. The expression of some genes annotated as F3′H,F3′5′H,DFR and UFGT also showed higher levels inDRwhen compared with X74 (Fig.8;Appendix S). Some genes encoding transporter proteins of the MATE efflux family were shown to be DEGs with the same expression trends as the anthocyanin synthesis structural genes,implying the involvement of MATEs in the vesicular transport of anthocyanins. Interestingly,in addition to the previously reported GSTF12 (GH_A07G0814andGH_D07G0816) involved in anthocyanin accumulation in cotton,several other GST genes were found to be differentially expressed between X74 andDR.Furthermore,four R2R3 MYB genes (GH_A07G0850,GH_A07G0851,GH_D07G0852andGH_D07G0853) and two bHLH genes (GH_D05G1649andGH_D08G2399)were highly up-regulated inDR(Appendix T).

4.Discussion

4.1.The function of the BBX lV sub-family genes

The BBX family is divided into five sub-families (Gangappa and Botto 2014). In Upland cotton,this family has a total of 123 members (Appendices U and V),andGH_A09G2280is a homologous gene of the BBX IV memberAtBBX24. The BBX IV subfamily includes 35 cotton genes and nineArabidopsisgenes (AtBBX18toAtBBX26).AtBBX21can directly bind to the promoter ofHY5through its second B-box to control the expression ofHY5,a key gene related to photomorphogenesis,and consequently the genes regulated byHY5in order to promote photomorphogenesis (Xuet al.2018).AtBBX21can also post-transcriptionally promote the level ofHY5,whileAtBBX24interferes with the binding ofHY5to the promoter of anthocyanin biosynthesis genes by heterodimerizing withHY5. The opposing functions of these two B-box genes have been shown to be specified by their C-terminal sequences,as AtBBX24 contains a typical VP domain at the C-terminus,whereas AtBBX21 does not (Yadukrishnanet al.2018).

In the red pear mutant,a 14-nucleotide deletion was found in the coding region ofPpBBX24. This deletion caused a coding frameshift and the formation of a premature stop codon,leading to loss of the NLS and VP domains inPpBBX24(Ouet al.2020). In this study,the mutation (a 2-bp deletion causing the loss of both NLS and VP domains) in theDRmutant and its consequence are very similar to those of the red pear mutant. While the mutant plants ofArabidopsisBBX24have a normal height,theAtBBX21overexpression plants are dwarf and have significantly increased accumulation of anthocyanins(Jobet al.2018;Xuet al.2018;Burschet al.2021). In this study,while the green cotton plants with the unmutatedGhDRdevelop normally,theDRmutant with the mutatedGhDRshowed enhanced expression of genes related to gibberellin and anthocyanin metabolism(Figs.7 and 8) and dwarf stature. Therefore,just like the importance of the C-terminus in specifying the functions of AtBBX21 and AtBBX24,the C-terminus with the NLS and VP domain seems to be critical for the functionality ofGhDR.

4.2.BBX and anthocyanidins

Transcriptome analysis indicated that enhanced expression ofGH_A09G2280inDRactivated the expression of many genes involved in multiple pathways,including the structural genes encoding anthocyanin biosynthetic enzymes (PAL,C4L,F3′H,F3′5′H,DFR and UFGT),two transport genes (GhGSTF12),and four regulatory genes encoding MYB113 transcriptional factors(TFs). Among them,GH_A07G0850/GhPAP1AandGH_D07G0852/GhPAP1Dhave been reported to positively regulate the biosynthesis of anthocyanins in Upland cotton (Liet al.2019;Lianget al.2020;Shaoet al.2022).Similarly,the red leaf trait observed in theDRmutant was probably caused by the activation of the anthocyanin pathway. The difference is that the upregulation ofGhPAP1AandGhPAP1Dis due to the presence of additional TF binding sites in their own promoters,whereasDRactivates MYB113 to positively regulate the anthocyanin signaling pathway by up-regulating the expression ofHY5(GH_D08G2693).

There are indeed examples of BBX family genes that induce anthocyanin accumulation. For example,SlBBX20was reported to promote the transcription of the anthocyanin biosynthesis geneSlDFRby binding to the G-box1 in its promoter (Luoet al.2021). The overexpression ofBBX23in poplar activated the expression of MYB TFs and structural genes in the flavonoid pathway,thereby promoting the accumulation of proanthocyanidins and anthocyanins (Li Cet al.2021).CRISPR/Cas9-mediated knockout ofBBX23resulted in the opposite outcome. Furthermore,the phenotype induced byBBX23overexpression was enhanced by exposure to high light conditions. BBX23 was capable of binding directly to the promoters of proanthocyanidinand anthocyanin-specific genes,and its interaction with HY5 enhances the activating activity (Li Cet al.2021). Such examples are numerous,PpBBX16is the homologous gene ofAtBBX22,and it cooperates withPpHY5and positively regulates light-induced anthocyanin accumulation by activating MYB10 in red pear. The transcript levels ofPpMYB10andPpHY5in red samples were significantly higher than those in green samples,whereas the results for the normal-typePpBBX24gene were the opposite (Zouet al.2017). In apples,MdBBX1/20/22promote anthocyanin biosynthesis in response to ultraviolet-B (UV-B) radiation or low temperature (Baiet al.2014;Anet al.2020).

4.3.BBX and plant height

BBX24 is a negative regulator of the red-light,far-red light and blue light-mediated inhibition of hypocotyl elongation(Indorfet al.2007). The interaction between the BBX24 transcription regulator and the DELLA protein regulates the activity of PIF4,a key transcription factor in the shadow signaling network. Both BBX24 and DELLA are negative regulators of the GA signaling pathway (Croccoet al.2015;Burschet al.2020),providing evidence that HY5 requires the B-box containing proteins (BBX20,BBX21 and BBX22) for transcriptional regulation. BBX21 binds to the promoter ofHY5in order to promoteHY5transcription,thereby leading to the accumulation of more HY5 protein (Zhaoet al.2020). HY5 has been shown to target many GA metabolism genes (Leeet al.2007). In peas,the HY5 ortholog LONG1 negatively regulates GA levels by inducing the expression ofGA2ox2,a key gene for GA catabolism (Welleret al.2009). In addition,the GA2ox2 enzyme inactivates most active GA,resulting in increased DELLA activity and reducedPIF4transcription activity (Fenget al.2008;de Lucaset al.2008). In this study,the genes encoding GA2ox2 (GH_A13G1786,GH_D01G0359) and DELLA (GH_A06G0641,GH_D06G0611)were significantly expressed in theDRmutant. The typical dwarf phenotype ofDRwas restored by spraying with gibberellin,but not with IAA or BR.

4.4.Hypothesis for the formation of the DR phenotype

The induced expression ofGhDRin theDRmutant prompted us to silenceGhDRin the mutant by virusinduced gene silencing (VIGS) to rescue the mutant phenotypes,i.e.,changing dwarf and red back to normal and green. Unfortunately,the expected phenotypes were not observed despite several attempts. One possibility is the presence of the Dt-subgenome homoeolog (GH_D09G2215) ofGhDR/GH_A09G2280. However,given that the RNAi constructs used in VIGS are expected to target both homoeologs due to their highly homologous sequences,this possibility is less likely. Instead,the main reason for the failure to observe the expected phenotypes might be that the expression ofGhDRcould not be completely silenced by VIGS to induce total loss-of-function and the threshold for the dwarf/red mutant phenotype expression is low;in other words,the presence of a very low level of the mutated transcript/protein ofGhDRis enough to trigger and/or maintain the mutant phenotypes.

We showed that the deletion mutation inGhDRcosegregates with the dwarf and red phenotype. While we were unable to reverse the dwarf and red phenotype back to normal and green phenotype by VIGS,we speculate that the 2-bp deletion inGhDRis the causative mutation responsible for the dwarf and red mutant phenotype based on the results of this study and previous studies onGhDRhomologs in other plant species. The enhanced expression ofHY5caused by the elevated level ofGhDRdue to the mutation in its 3′ terminus by an unidentified mechanism might be associated with the enhanced biosynthesis and accumulation of anthocyanins,leading to the red mutant phenotype. On the other hand,the reduced degradation of DELLA protein caused by the mutation inGhDRmight have repressed GA biosynthesis,leading to the dwarf mutant phenotype.

These positive and negative regulatory roles ofGhDRseem to be achievedviaHY5 (Fig.9),which has been demonstrated to play an essential role in photomorphogenesis and the accumulation of anthocyanins (Jobet al.2018;Xu 2020). HY5 reportedly binds to transcription factors such as MYB (Shinet al.2013),bHLH (Liet al.2020),and WD40 (Wanget al.2020)to promote the accumulation of anthocyanins. Consistently in this study,the expression levels of HY5 and its potential targets,such as MYB,bHLH,and WD40,were significantly increased inDR(Appendix T),suggesting a role for the interaction betweenGhDRand HY5 in the formation of the red phenotype. In addition,the DELLA protein can reportedly prevent TCP (Daviereet al.2014),ARFs (Ohet al.2014),PIFs (de Lucaset al.2008),and BZR1 (Liet al.2012;Liuet al.2018) from functioning by sequestering them. In this study,we found that the expression levels of the genes encoding DELLA proteins were significantly up-regulated inDR,leading to repression of the transcription factors such as TCP,ARFs,PIFs,and BZR1 that are required for the metabolic pathway of gibberellin biosynthesis and development of the normal plant. This result suggests thatGhDRmay be indirectly involved in the regulation of GA biosynthesis including DELLA and its TF targets,ultimately affecting plant architecture. Taken together,these results imply thatGhDRsimultaneously regulates plant architecture and anthocyanin accumulation by regulating a swath of TFs that are direct and/or indirect targets of HY5 and DELLA (Fig.9).

This study provides new genetic resources for further use in molecular breeding methods to reasonably regulate and moderately express candidate genes related to plant architecture. However,the biological functions ofGhDRin photomorphogenesis still need to be confirmed,and its downstream molecular mechanisms need to be further explored.

5.Conclusion

In this study,we characterized a natural cotton mutant with extremely dwarfing plant height and red leaves,and the mutant phenotype was controlled by a single dominant locus. A 2-nucleotide deletion was found in the coding region of the candidate geneGhDR/GH_A09G2280,resulting in a frameshift and truncation ofGhDR. A comparative transcriptome analysis demonstrated thatGhDRsimultaneously regulates plant architecture and anthocyanin accumulation by regulating a broad range of TFs. This study provides the new genetic resources for further exploiting molecular breeding related to plant architecture.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (32160477 and 31960412),the International S&T Cooperation Projects of BINGTUAN,China (2021BC001),and the Young and Middle-aged Leaders in Scientific and Technological Innovation Foundation of Shihezi,China (2021RC02 and 2020CB010).

Declaration of competing interest

The authors declare that they have no conflict of interest.

Appendicesassociated with this paper are available on https://doi.org/10.1016/j.jia.2022.10.007