Mutations in the miRNA165/166 binding site of the HB2 gene result in pleiotropic effects on morphological traits in wheat

Dengji Jing,Lei Hu,Chozhong Zhng,Hongn Li,Zheng Wng,Jin Li,Guiping Wng,Rui Song,To Shen,Hongyu Li,Shengsheng Bi,Ynn Liu,Jin Wng,Ho Li,Jorge Duovsky,Shisheng Chen,*

a Peking University Institute of Advanced Agricultural Sciences,Weifang 261000,Shandong,China

b Department of Plant Sciences,University of California,Davis,CA 95616,USA

c State Key Laboratory of Crop Stress Adaptation and Improvement,College of Agriculture,Henan University,Kaifeng 475004,Henan,China

Keywords:Wheat miRNA165/166 Curled leaves Paired spikelets Homeodomain-leucine zipper transcription factor

ABSTRACT Leaf,spike,stem,and root morphologies are key factors that determine crop growth,development,and productivity.Multiple genes that control these morphological traits have been identified in Arabidopsis,rice,maize,and other plant species.However,little is known about the genomic regions and genes associated with morphological traits in wheat.Here,we identified the ethyl methanesulfonate-derived mutant wheat line M133 that displays multiple morphological changes that include upward-curled leaves,paired spikelets,dwarfism,and delayed heading.Using bulked segregant RNA sequencing(BSR-seq)and a high-resolution genetic map,we identified TraesCS1D02G155200(HBD2)as a potential candidate gene.HB-D2 encodes a class III homeodomain-leucine zipper(HD-ZIP III)transcription factor,and the mutation was located in the miRNA165/166 complementary site,resulting in a resistant allele designated rHb-D2.The relative expression of rHb2 in the mutant plants was significantly higher(P＜0.01)than in plants homozygous for the WT allele.Independent resistant mutations that disrupt the miRNA165/166 complementary sites in the A-(rHb-A2)and B-genome(rHb-B2)homoeologs showed similar phenotypic alterations,but the relative intensity of the effects was different.Transgenic plants expressing rHb-D2 gene driven by the maize UBIQUITIN(UBI)promoter showed similar phenotypes to the rHb-D2 mutant.These results confirmed that HB-D2 is the causal gene responsible for the mutant phenotypes.Finally,a survey of 1397 wheat accessions showed that the complementary sites for miRNA165/166 in all three HB2 homoeologs are highly conserved.Our results suggest that HB2 plays an important role in regulating growth and development in wheat.

1.Introduction

Bread or common wheat(Triticum aestivum;AABBDD genomes)is one of the most important cereal crops in terms of both cultivated area and total production(FAOSTAT,https://www.fao.org/-faostat/en/),and it plays a key role as a staple food for much of the world’s human population.Common wheat has undergone two recent allopolyploidization events,and as a result,the inheritance of morphological traits in wheat tends to be complex and quantitative[1].In addition,the large(approximately 17 Gb)and complex nature of the hexaploid wheat genome has limited indepth genetic studies of morphological traits.

Leaf,spike,stem,and root morphologies(henceforth referred to as plant architecture)are agriculturally important traits that have been extensively studied in Arabidopsis(Arabidopsis thaliana),rice(Oryza sativa),maize(Zea mays)and other plant species[2–4]because of their potential large effects on development and productivity.Scientists have expended considerable effort to identify the key genes that control multiple morphological changes,and thus,some genes related to plant architecture have been identified.In Arabidopsis,the REVOLUTA(REV)gene is necessary for apical meristem development and for limiting cell divisions in the stems and leaves[2].In rice,plants with a mutation in the miRNA165/166 target sequence of the transcription factor lateral florets 1(lf1)show three-floret spikelet and rolled-leaf phenotypes[3,5].Defects in the genes narrow leaf and dwarf 1(ND1),curled and dwarf leaf 1(OsCD1),and narrow and rolled leaf 1(OsCslD4)result in dwarfed plants with narrow and curled leaves[6–8].In maize,rolled leaf1(rld1)encodes a homeodomain leucine zipper class III(HD-ZIP III)protein that controls the upward rolling of the leaf blade[4].

Previous studies have shown that HD-ZIP III transcription factors play important roles in the embryo[9],spike[3],root[10],and shoot and leaf development[11].It is well documented that HD-ZIP III members are regulated by miRNA165/166 and are required for vascular development,shoot apical meristem(SAM)establishment,and polarity formation in lateral organs[12,13].The HD-Zip III family in Arabidopsis contains five genes:REVOLUTA(REV),PHABULOSA(PHB),CORONA(CNA),PHAVOLUTA(PHV),and

ATHB8[14].These genes have redundant or antagonistic roles in Arabidopsis development[11].The HD-Zip III family in rice also consists of five members,including OsHox10,OsHox9,OsHox33,OsHox32,and OsHox29(also named OsHB1 to OsHB5)[15].In contrast,there are few studies in wheat on the identification of genomic regions/genes involved in determining plant architecture.

In the present study,we identified and characterized a mutant line,M133,derived from an EMS-mutagenized population of the wheat cultivar‘Ningchun 4’[16].M133 displayed multiple morphological changes,such as upward-curled leaves,paired spikelets,dwarfism,and delayed heading.A single genetic locus,temporarily designated as Abnormal Plant Architecture-1(APA1),was found to be responsible for the mutant phenotypes.The objectives of this study were to:(1)characterize and genetically map APA1;(2)identify potential candidate genes in the sequenced wheat genome;and(3)validate candidate genes using independent EMS mutants and transgenic plants.

2.Materials and methods

2.1.Plant materials and mapping population

An EMS-mutagenized population of the hexaploid Chinese wheat variety‘Ningchun 4’(NC4)was previously developed to clone the male sterility gene MS1[16].The initial screen of this EMS mutant population was carried out in the 2009–2010 growing season.We identified an M2mutant line(M133)displaying upward-curled leaves,paired spikelets,dwarfism,delayed heading,and decreased fertility.Subsequently,this mutant(M3generation)was grown in a growth chamber environment(25 °C day/23 °C night;18 h light/6 h dark photoperiod),and it was stably inherited through two generations of self-pollination.To genetically characterize this mutant line,M5-generation individuals were crossed with the wild-type NC4,generating an F2population comprised of 168 individuals.To eliminate the influence of background mutations,we generated another BC2F2population consisting of 500 individuals.Phenotypic data from the F1,F2,F2:3,and BC2F2populations was used to map the causal gene.The F2:3families were sown in one-meter-long rows(～25 seeds per family)in the field at Peking University Institute of Advanced Agricultural Sciences,Weifang,Shandong province,China(36°26′04.0′′N,119°26′42.6′′E).

2.2.Bulked segregant analysis-RNA-Seq(BSR-seq)

Based on the phenotypic information,15 homozygous wildtype and 15 homozygous mutant F2:3families were grown at 25 °C day/23 °C night with an 18 h light/6 h dark photoperiod for 14 days,and pooled tissues were collected to construct the wild-type bulk(W-bulk)and mutant bulk(M-bulk)for RNA extraction.Total RNA of the parents and bulks was extracted using the Spectrum Plant Total RNA Kit(MilliporeSigma,St.Louis,MO,USA).RNA-seq was performed on an Illumina HiSeq 4000 platform at Beijing Novogene Bioinformatics Technology Co.,Ltd.(Beijing,China).Raw sequencing reads were trimmed using the software Trimmomatic v 0.32[17]to remove low-quality reads and adaptors.Trimmed reads were aligned to the published reference genome of the hexaploid wheat‘Chinese Spring’(CS)RefSeq v1.1[18]using the software STAR version 2.5.0c[19].Variant calling was carried out using the HaplotypeCaller tool from GATK v3.2–2[20].The single nucleotide polymorphism(SNP)-index andΔ(SNP-index)[21,22]were calculated to identify candidate genomic regions.

2.3.Marker development

SNPs associated with the APA1 locus were identified by BSR-seq analysis and were selected for the development of Kompetitive Allele-Specific Polymerase Chain Reaction(KASP)or Cleaved Amplified Polymorphic Sequence(CAPS)markers.The methods and procedures for developing KASP and CAPS markers have been reported previously[23,24].The sequences flanking the target SNPs were obtained from the‘Chinese Spring’reference genome[18].Genome-specific primers were developed using the software Primer3 web version 4.1.0(https://primer3.ut.ee/)to amplify gene regions containing the target SNPs.NEBcutter v3.0.15(https://nc3.neb.com/NEBcutter/)was used to determine the restriction enzymes that selectively cut the sequences containing the targeted SNPs.To obtain more SNPs between NC4 and the mutant line M133,we aligned the NC4 and M133 sequences from the RNAseq data,identified polymorphisms,and generated markers within the candidate gene region.

2.4.Phylogenetic analysis

After searching the sequenced genomes of rice(Oryza sativa,https://www.ricedata.cn/gene/),maize(Zea mays,https://maizegdb.org/),barley(Hordeum vulgare,https://plants.ensembl.org/Hordeum_vulgare/Info/Index),Brachypodium(https://phytozomenext.jgi.doe.gov/),Aegilops tauschii(https://plants.ensembl.org/Aegilops_tauschii/Info/Index),and the hexaploid wheat variety‘Chinese Spring’(https://wheat-urgi.versailles.inra.fr/Seq-Repository/BLAST)using the sequences of HD-Zip III proteins from Arabidopsis as queries,we obtained 43 protein sequences for phylogenetic analysis.All protein sequences were aligned with Muscle as implemented in the software package MEGA version 7[25],and phylogenetic trees were then produced using the Neighbor-Joining method.Finally,the software Interactive Tree Of Life(iTOL)v.5 was used to visualize the phylogenetic tree(https://itol.embl.de/).

2.5.Quantitative reverse transcription PCR(qPCR)analysis

Plants were grown in growth chambers with the following temperature and photoperiod:25°C day/23°C night,18 h light and 6 h dark.Leaves and roots were collected from BC2F3sister lines at the two-leaf stage.Young spikes(～10 mm)and stems were sampled at later stages.Total RNA was isolated individually from different wheat tissues(leaf,stem,root,and spike)using the Spectrum Plant Total RNA Kit(MilliporeSigma).RNA samples were purified using the Direct-zol RNA MiniPrep Plus kit(Zymo Research,Irvine,CA,USA).To eliminate DNA contamination,we performed the oncolumn RNase-Free DNase I digestion following the manufacturer’s instructions.After examining the integrity of the RNA samples,the concentrations were determined using a NanoDrop OneC spectrophotometer(Thermo Fisher Scientific,Waltham,MA,USA).

First-strand complementary DNA(cDNA)was synthesized from 1 μg of total RNA using ProtoScript II First Strand cDNA Synthesis Kit(New England Biolabs,Hitchin,UK).qPCR reactions were performed on an ABI QuantStudio 5 Real-Time PCR System(Applied Biosystems,Foster City,CA,USA)with PowerUp SYBR Green Master Mix(Thermo Fisher Scientific).The ACTIN gene was used as the endogenous control for normalization of gene expression[26].Transcript levels were presented as fold-ACTIN levels using the 2-ΔCTmethod as described previously[27].

2.6.EMS mutants

The sequenced EMS mutagenized populations of the hexaploid cultivar‘Cadenza’[28]are available on line(https://www.wheattilling.com/).Another sequenced EMS mutagenized population of the Chinese hexaploid cultivar‘Jimai 22’(https://jm.ytjebc.com/)was also used in this study.We screened these two databases for mutations using BLASTN with the sequences of TraesCS1D02G155200,

TraesCS1A02G157500,and TraesCS1B02G173900 as search queries.The plants from the selected‘Cadenza’mutant lines were grown in a greenhouse at the University of California,Davis(CA,USA).

2.7.Energy of miRNA165/166-target interaction

The predicted energy of the interactions between miRNA165/166 and the target sites were determined using the RNAhybrid on line tool[29].For this analysis,we used the miRNA165/166 mature sequence and as targets the wild type HB2 and the mutant alleles M133,Cadenza1761 and Cadenza0269,including 12 nt at each side,which are conserved in the three homoeologs.

2.8.Wheat transformation

We amplified the coding region of the candidate gene from the mutant line M133 using PrimeSTAR Max DNA Polymerase(TaKaRa,Tokyo,Japan),and generated the DNA construct for plant transformation in the binary vector pCAMBIA1300.The primer pair HBCD647 was used to amplify the rHb-D2 coding region(Table S1).The amplified PCR product was gel purified using the NucleoSpin Gel and PCR Clean-up kit(Macherey-Nagel,Germany)and recombined into the linearized pCAMBIA1300 vector using the In-Fusion HD Cloning Kit(Clontech,Mountain View,CA,USA)following the manufacturer’s instructions.The construct pCAMBIA1300-rHb-D2 contained the maize UBIQUITIN(UBI)promoter fused to the amplified rHb-D2 coding region(2517 bp),hereafter referred to as UBI:-rHb-D2.The resulting plasmid was transformed into the hexaploid wheat variety‘Fielder’using Agrobacterium tumefaciens strain EHA105-mediated transformation at the Peking University Institute of Advanced Agricultural Sciences transformation facility.The primer pair HB-JCX developed from the rHb-D2 sequence was used to confirm the presence of the transgene in potential transgenic plants.

2.9.Allelic variation

The published reference genomes of diploid[30,31],tetraploid[32,33]and hexaploid[34,35]wheat varieties were used to detect allelic variations in the HB2 candidate gene.In addition,we also used exome-capture data from three wheat panels to study variations in the HB candidate gene.The three wheat panels consist of 811(after imputation and filtering)[36],512[37],and 59 wheat accessions[38],respectively.

3.Results

3.1.Characterization of the EMS-derived wheat mutant line M133

Based on the phenotyping data,M5plants of the mutant line M133(apa1)exhibited multiple phenotypic changes including upward-curled leaves,paired spikelets,and dwarfism,when compared to the wild-type(WT)control plants(Fig.1).To study this mutant line further,we performed a genetic analysis using F1,F2,and F2:3populations from the cross NC4×M133.We first evaluated the rolled-leaf phenotype,which is a prominent feature of the mutant plants(Fig.1B,D).We observed that the F1plants derived from crossing the parental lines showed an intermediate phenotype(Figs.1D,S1).Among the 168 F2plants and their corresponding F2:3families,42 families showed the homozygous rolledleaf phenotype,88 were segregating,and 38 were homozygous for the normal WT leaf phenotype,which fits the 1:2:1 segregation ratio expected for a single genetic locus(χ2=0.57,P=0.75).

In the same F2and F2:3populations,we also investigated plant height and paired spikelets,which also segregated and showed obvious phenotypic differences(Fig.1C,E).The data showed that the rolled-leaf phenotype is completely linked with the plant height and paired spikelet phenotypes,indicating that they may be controlled by the same gene.To reduce potential variability derived from background mutations,we crossed M133 mutant plants with WT NC4 and then backcrossed two times to NC4.The BC2F1plants were self-pollinated and a BC2F2population consisting of 500 individuals was generated.In this population,we investigated multiple agronomic traits,including rolled leaves,plant height,and paired spikelets,and found that these traits were significantly correlated with each other(P＜0.001).In addition,we observed that the fertility of the homozygous mutant plants in both F2and BC2F2populations was greatly reduced(P＜0.001)when plants were grown at 25 °C day/23 °C night with an 18 h light/6 h dark photoperiod.

3.2.Validation of the phenotypic changes using BC2F3 sister lines

We selected BC2F3sister lines homozygous for the WT and mutant alleles to further determine the phenotypic changes.At maturity,BC2F3plants homozygous for the mutant allele showed a significant reduction in leaf width(49.6%)and length(31.8%)relative to those homozygous for the WT allele(Fig.2A,B).The main stems in the BC2F3mutant plants were significantly shorter than in the WT sister lines(P＜0.0001;Fig.2C).The reduced stem length was associated with decreases in the length of the third,fourth,and fifth internode by 28.5%,52.4%,and 75.1%,respectively,while the lengths of the first and second internode were not significantly different(Fig.2D).For spike morphological traits,in addition to the BC2F3plants homozygous for the WT and mutant alleles,we measured heterozygous BC2F3plants because they also showed paired spikelets(Fig.1E).We observed that the plants carrying the mutant alleles displayed abnormal spikes(Fig.1E)where the awns and spikes were trapped in the rolled leaves.There was no significant difference in the average number of primary spikelets per spike among the three allelic classes(Fig.2E).However,we found that the heterozygous and homozygous mutant BC2F3plants had an average of 8.1 and 11.5 secondary spikelets(paired spikelets),which were not observed in the WT sister line(P＜0.0001;Fig.2F).There was no significant difference in spike length between the WT and heterozygous plants,but spikes from plants homozygous for the mutant allele were significantly shorter(P＜0.0001;Fig.2G).Finally,we examined 6-day-old wheat seedlings(Fig.1F)and found a decrease in the number of seminal roots(6.8%,P＜0.05;Fig.2H)and length of the primary roots(12.9%,P＜0.05;Fig.2I)in the BC2F3plants homozygous for the mutant allele relative to those homozygous for the WT allele.

Fig.1.Phenotypic changes in M133 mutant plants.(A)WT and M133 plants grown in a controlled environment.(B)Close-up views of the main tillers.(C)A close-up of the internodes.Dashed white lines represent the positions of the nodes.(D)The rolled-leaf phenotype in BC2F3 plants homozygous for Apa1 and apa1 and heterozygous plants.(E)Spike phenotypes(red arrows indicate paired spikelets).(F)Root phenotypes(6-day-old seedlings).Plants were grown at 25°C day/23°C night with an 18 h light/6 h dark photoperiod.WT,wild type;Apa1,BC2F3 plants homozygous for the WT allele;apa1,BC2F3 plants homozygous for the mutant allele;Het.,heterozygous(Apa1/apa1)BC2F3 plants.

When BC2F3plants carrying the WT and mutant alleles were grown at 25 °C day/23 °C night with an 18 h light/6 h dark photoperiod,we observed that BC2F3plants carrying the mutant allele flowered 2.3 days later than those carrying the WT allele(P＜0.01,Fig.S2).To further confirm this result,we performed another experiment at a lower temperature with a shorter photoperiod(10 °C/12 h light and 12 h dark).Statistical analyses showed that BC2F3plants homozygous for the mutant allele headed 8.8 days later(P＜0.0001)than those carrying the WT allele under these conditions(Fig.S3),a larger delay than in the previous 25°C day/23°C night long-day experiment.Using the 168 F2:3families described above,we confirmed that heading time was linked to the other phenotypic changes.In addition,we found that the intensity of the morphological defects in the stem and leaf organs(rolled leaf and dwarfism)in the BC2F3mutant plants was weaker when plants were grown at 12 h light/10 °C than when grown at 18 h light/25 °C day/23 °C night(Fig.S3A,B).This result suggests that some of the phenotypic defects in mutant plants may be sensitive to photoperiod or/and temperature.Taken together,these multiple linked mutant phenotypes indicate that the M133 mutated gene has multiple effects on wheat growth and development.

3.3.Genetic mapping of a locus on chromosome 1D by BSR-seq

To map the causal gene,we carried out bulked segregant RNAseq(BSR-seq)on an F2:3mapping population derived from the cross of NC4×M133.RNA sequencing reactions generated,on average,50.1 million 150-base raw read pairs per sample(NCBI SRA database accession number PRJNA826248;NC4,M133,Wbulk,and M-bulk).After trimming the reads and filtering them for quality and adapter contamination,～99%of the reads remained and were aligned to the‘Chinese Spring’wheat reference genome RefSeq v1.1.Approximately 82% of the filtered read pairs mapped uniquely to the reference genome.Subsequent SNP calling identified 4233 high confidence EMS-induced SNPs between the two parents.Based on SNP index analysis,we detected 58 SNPs that were significantly associated with the phenotype.Of these 58 significant SNPs,22(37.9%)were located on chromosome 1D(Fig.3)and were mainly distributed in a region from 234.7 to 324.4 Mb(CS RefSeq v1.1,Table S2)suggesting that the APA1 locus was likely located on chromosome 1D.

Fig.2.Effect of the apa1 mutation on leaf,spike,root,and stem development in wheat.(A)Leaf width;(B)Leaf length;(C)Stem length;(D)Internode lengths;(E)Spikelet number per spike(primary spikelet);(F)Number of paired spikelets(secondary spikelets);(G)Spike length;(H)Root number;and(I)Primary root length.Bars are averages and error bars are standard errors of the mean(SEM).Numbers(n)inside the bars indicate the number of plants examined.The asterisks indicate the statistical significance of differences(Dunnett test;ns,not significant;*,P＜0.05;**,P＜0.01;and***,P＜0.001).Apa1,BC2F3 plants homozygous for the WT allele;apa1,BC2F3 plants homozygous for the mutant allele;Het.,heterozygous BC2F3 plants.

3.4.Isolation of the causal gene

To confirm the mapping location,we selected six SNPs(C234740074T,C238067979T,C240594830T,C259211449T,C299584512T,and C365718293T)surrounding the region of interest on chromosome 1D,and developed six D-genome specific CAPS or KASP markers based on these SNPs(Table S1).Using these markers,we genotyped the 168 F2plants that had been previously phenotyped and mapped the causal gene to a locus 0.6 cM proximal to the CAPS marker loci JY2374,JY2380,and JY2431(Figs.4A,B,S4).To define better the position of APA1,we then developed two new molecular markers using two SNPs(G201643118A and G217637192A)located proximal to marker locus JY2374(C234740074T)and added them to the genetic map(Fig.4B).Based on this population,APA1 was tightly linked to marker locus HB-FR2,and it mapped to a 1.2 cM region between the marker loci 2018FAB and JY2374(Fig.4B).

To delimit the position of the APA1 locus more precisely,we screened another 500 BC2F2plants using the flanking markers 2018FAB and JY2374.We identified 12 plants with informative recombination events in these 500 plants and obtained four recombinants between the same flanking markers used in the previous 168 plants.Using these critical recombinants and three new markers developed in this candidate region,APA1 was fine mapped to a 0.14 cM interval flanked by marker loci KP2121 and KP2150(Figs.4C,S4).The 0.14 cM candidate region defines a 3.0 Mb region(215.7–218.7 Mb)in the CS RefSeq v1.1 genome that contains only 11 annotated high-confidence genes(Table S3).Based on BLASTN searches in the wheat expVIP database(https://www.wheatexpression.com/),we found that nine candidate genes were expressed in various wheat tissues(Table S3).These expressed genes included one gene,TraesCS1D02G155200,which was annotated as encoding a class III HOMEOBOX LEUCINE-ZIPPER(HDZIP III)protein,which is of particular interest because HD-ZIP-type proteins have been reported to be involved in the regulation of leaf,spike,and stem development in the other plants[3,11,12].

Fig.3.Distribution of candidate SNPs associated with the mutant phenotype.A total of 22 significant SNPs were located on chromosome 1D.

Fig.4.Map-based cloning of the causal gene.(A)Colinear region on chromosome 1D from the‘Chinese Spring’reference genome(RefSeq v1.1).Arrows indicate genes.(B)Genetic mapping of APA1 using 336 segregating gametes.(C)High-density map of APA1 using 1336 segregating gametes and prediction of the candidate gene.(D)Gene structure of the candidate gene TraesCS1D02G155200(HB-D2).Black rectangles and lines represent exons and introns,respectively;the mutated site in the mutant line M133 is shown in red.(E)Alignment of miRNA166 with the target sites(cDNA coordinates)in NC4 and M133.(F)Sequencing chromatograms showing the G575A polymorphism between the NC4 and M133 alleles.The mutated nucleotide is indicated by a red arrow.

Within the candidate chromosomal region,there are only two mutations(G216792764A and G217637192A)between the EMS induced M133 and the parental NC4 line based on RNA sequencing data.Both SNPs co-segregated with the mutant phenotypes(Table S1).The first SNP,G216792764A,is unlikely to be the causal mutation because it is located within an uncharacterized gene that was not annotated in the CS reference genome,while the second SNP,G217637192A,is located within the candidate gene TraesCS1D02G155200.Using the primer pair HB-FR2(Table S1),we performed PCR to amplify the region including the mutation and confirmed the presence of a single nucleotide transition from G to A(G575A,cDNA coordinates)in the mutant line M133,which results in an amino acid change from glycine to glutamic acid at position 192 from the starting methionine(G192E;BLOSUM62 score=-2;Fig.4D,F).Further analysis revealed that the mutated site is located in a putative miRNA165/166 target sequence(Fig.4E),which agrees with the partially dominant nature of the mutation.This semi-dominant resistant allele will be designated hereafter as rHb-D2. These results indicate that TraesCS1D02G155200 is the best candidate for the causal gene among the genes present within the candidate gene region.

3.5.TraesCS1D02G155200 gene structure,phylogenetic analysis,and expression pattern

In NC4,the candidate gene TraesCS1D02G155200 spans 8,463 bp from the start to the stop codon with a complete coding sequence of 2517 bp(Fig.4D).This gene is composed of 18 exons interrupted by 17 introns,and it encodes a predicted protein of 838 amino acids.A Neighbor-Joining phylogenetic analysis was carried out to determine the evolutionary relationships between TraesCS1D02G155200 and 42 other HD-Zip III proteins from wheat and six other plant species(Fig.S5).Based on the phylogenetic analysis,the TraesCS1D02G155200 protein was grouped in a clade with OsHB2 from Oryza sativa with a high degree(92.6%)of sequence identity(Fig.S5).TraesCS1D02G155200 is an orthologue of the rice gene OsHB2 and a paralog of lateral florets 1(lf1)and rolled leaf1(rld1)in maize[3,4],and is designated here as HB2(HB-D2).The A-(HB-A2)and B-genome(HB-B2)homoeologs are very similar to HB-D2 with 98.9% and 99.1% sequence identity,respectively.

In general,miRNAs guide the cleavage of target mRNAs by binding to complementary sites[39].Mutations that disrupt the miRNA complementary sites are known to influence miRNA-mediated cleavage in vivo[3,4,40,41].We analyzed HB-D2 transcript levels relative to ACTIN in BC2F3plants homozygous for the WT and mutant alleles by qPCR.The expression data for HB-D2 in leaves,roots,stems,and spikes showed that its transcripts accumulated in all tissues at significantly higher levels in the plants carrying the resistant rHb-D2 allele compared to the Hb-D2 WT allele(P＜0.01;Fig.5).This result suggests that,due to the disruption of the miRNA165/166 complementary site,the mRNA from the resistant rHb-D2 allele is not degraded by miRNA165/166 as efficiently as the WT allele.Transcript levels of HB-A2 and HB-B2 were also evaluated using two genome-specific qRT-PCR primers(Table S1),and they showed similar levels to HB-D2 in the BC2F3plants carrying the WT alleles(Fig.S6).Finally,analysis of RNAseq data available in the wheat expVIP database(https://www.wheat-expression.com/)showed similar expression profiles(Fig.S7).Transcripts of HB-D2,HB-A2,and HB-B2 were detected in different organs and tissues.

3.6.The rHb-A2 and rHb-B2 EMS mutants show phenotypic alterations

A search of the sequenced EMS-mutagenized populations derived from wheat cultivars‘Cadenza’and‘Jimai 22’yielded two mutants for HB-D2.Mutant lines Cadenza1290 and jm_chr1D_217637192 carry the same mutation G217637192A as in M133(Fig.S8).However,seeds of both mutant lines were not available likely due to the reduced fertility of the mutants.

We then searched for mutations in the miR165/166 binding site of HB-A2 and HB-B2.We obtained one‘Cadenza’mutant line,Cadenza1761,for HB-A2,which carries the mutation G580A(G194R;BLOSUM62 score=-2)located in the miRNA165/166 complementary site(Figs.6A,S9).For the B-genome homoeolog HB-B2,we identified the‘Cadenza’mutant line Cadenza0269 that has a single nucleotide change C589T(P197S;BLOSUM62 score=-1)in the miRNA166 complementary site(Figs.6B,S9).Using the genome-specific primer pairs CaA1761 and CaB0269(Table S1),we genotyped 20 plants from each of the selected mutant lines.We identified six plants that were homozygous for the rHb-A2 allele,nine plants homozygous for the WT allele,and five heterozygous plants in the line Cadenza1761.In the mutant line Cadenza0269,we detected four homozygous mutant plants,14 homozygous WT plants,and two heterozygous plants.In a greenhouse experiment,we observed that the plants homozygous for the rHb-A2 mutant allele exhibited paired spikelets and reduced plant height(Fig.6A).But we could not determine whether these homozygous rHb-A2 plants formed upward-curled leaves since the plants homozygous for the WT allele in this mutant line also showed some degree of leaf curling(Fig.6A).

For the plants carrying the rHb-B2 allele,we also observed similar mutant phenotypes(Fig.6B),but the intensity of these phenotypic defects was lower than in the rHb-A2 homozygous mutant plants(Fig.6).As shown in Fig.6B,some of the secondary spikelets developed into needle-like structures in the rHb-B2 mutant.The smaller phenotypic effects of the rHb-B2 mutation relative to rHb-A2 may be related to the smaller disruptive effect of the predicted amino acid changes,since the predicted minimum free energy(mfe)of the interactions between miRNA and the target sites were very similar(Fig.S10).The larger difference in mfe in M133 may have contributed to its stronger phenotype(Fig.S10).

3.7.rHb-D2 transgenic plants exhibit obvious abnormalities

Previous studies have reported that overexpression of the mutated resistant versions of HD-Zip III genes resulted in abnormalities such as rolled leaves and lateral florets from ectopic expression of these genes in corresponding tissues[3,5,13,42].To determine whether rHb-D2 has the same effects,we generated 12 independent T0transgenic plants expressing rHb-D2 driven by the maize UBI promoter.The presence of the rHb-D2 transgene in these plants was confirmed using the primer pair HB-JCX(Table S1).Transcript levels of rHb-D2 were significantly higher(P＜0.001)in all transgenic plants than in the non-transgenic‘Fielder’control plants(Fig.S11).Transgenic plants grown in a growth chamber environment(25 °C day/23 °C night and 18 h light/6 h dark photoperiod)displayed severe leaf rolling and extreme dwarfism(Fig.7A,B).At maturity,some of the transgenic plants(e.g.mT0-09 and mT0-12)showed paired spikelets(Fig.7C).In addition,plants from T1transgenic families mT1-09 and mT1-12 carrying the transgene showed the rolled-leaf phenotype(Fig.7D).

The extreme phenotype of the transgenic plants is likely the result of a combination between the high expression levels generated by the UBI promoter and the reduced mRNA degradation generated by the mutation in the miRNA165/166 binding site in the constitutively expressed the resistant rHb-D2 allele.Taken together,the high-resolution genetic map,the independent mutants, and the transgenic plants confirm that TraesCS1D02G155200(HB-D2)is the causal gene for the abnormal plant architecture associated with APA1 locus in the M133 mutant.

3.8.Natural variation in HB2

Fig.5.HB-D2 mRNA levels in BC2F3 sister plants homozygous for either the WT or mutant alleles.Plants were grown in growth chambers at 25 °C day/23 °C night with an 18 h light/6 h dark photoperiod.Leaves and roots were collected at the two-leaf stage,stems at the heading stage,and young spikes were collected when they were～10 mm in length.Transcript levels were normalized relative to ACTIN gene that was used as endogenous control(n=4).Error bars indicate the SEM.The asterisks indicate the statistical significance of the differences(ns,not significant;*,P＜0.05;**,P＜0.01;and***,P＜0.001).Hb-D2,BC2F3 plants homozygous for the WT allele;rHb-D2,BC2F3 plants homozygous for the mutant allele.

In order to determine the allelic variations in HB-D2,HB-A2,and HB-B2,a total of 15 published reference genomes and 1382 wheat genotypes with exome sequencing data derived from the 500 exomes project[37],the 1000 wheat exomes project(includes 811 wheat genotypes after imputation and filtering)[36],and the 59 exomes project(generated in WheatCAP project)[38],were used for comparative analysis.Our survey found that HB-D2 is highly conserved among different wheat genotypes with only four rare amino acid polymorphisms(Table S4).For HB-A2 and HB-B2,we detected 12 and 8 amino acid polymorphisms,respectively(Table S4).However,no SNPs were observed in the complementary sites for miRNA165/166 in the three HB2 homoeologs among these 1397 wheat accessions,indicating that these sequences are highly conserved.

4.Discussion

4.1.Identification of HB-D2 as the causal gene of the morphological changes in the M133 mutant

Leaf,spike,stem,and root morphologies are complex traits controlled by multiple genes.In rice and other plant species,many genes/QTL associated with these morphological traits have been mapped or cloned[4,7,43,44].Genetic intervals controlling leaf,spike,stem,and root morphological changes have been also reported in wheat,but it was previously difficult to isolate the causal genes due to the large and complex nature of the wheat genome[45–47].Fortunately,the recently released wheat genomes and the availability of new high throughput sequencing tools have made this task feasible.

BSR-seq,a recently-developed method that combines BSA with RNA-seq,has been used to identify major genes in a number of different plant species,including maize[48],wheat[49],cabbage[50],and soybean[51].Using BSR-seq analysis and the publicly available wheat genome sequences,we were able to rapidly map the causal gene region, which included HB-D2(TraesCS1D02G155200).This gene encodes a class III homeodomain-leucine zipper(HD-ZIP III)transcription factor,a gene family that has been associated with morphological changes in other plant species[3,4,11,12].

4.2.Effect of miRNA165/166 resistant alleles in different plant species

The mutations in many gain-of-function mutants of HD-ZIP III genes,such as rev,phb,and phv in Arabidopsis,rld1 in maize,and lf1 in rice,were all found to be located within the miRNA165/166 complementary site[3,4,52].In this study,we detected a single nucleotide substitution(G575A,Fig.4)in HB-D2 that disrupted the base pairing at the 3′end of miRNA165/166.This mutation caused multiple phenotypic changes,the most serious of which included upward-curling leaves,paired spikelets,and dwarfism(Fig.1).The same nucleotide substitution was reported in the miRNA165/166 complementary site of rld1 in the maize Rld1-O mutant[4].

The phylogenetic analysis of the HD-ZIP III proteins from seven plant species showed that HB-D2,RLD1,and LF1 cluster together in the same clade(Fig.S5).In maize,gain-of-function mutations in Rld1 were mainly involved in leaf polarity regulation[4].In rice,a single mutation within the miRNA165/166 binding site of the lf1 gene caused curled leaves and three-floret spikelets[3,5].We noticed that the upward curling of the leaves in the rHb-D2 mutant is similar to the leaf development defects in Rld1-O and lf1 mutant plants,indicating a conserved function for these HD-ZIP III genes in rice,maize,and wheat.However,in addition to the rolled-leaf phenotype,we also observed other morphological changes in rHb-D2 mutant plants that were not observed in the Rld1-O and lf1 mutants,suggesting that these genes may have some degree of functional differentiation in the different plant species.However,phenotypes such as the paired spikelets may be associated with the different architecture of the wheat inflorescence.

Gain-of-function phenotypes due to HD-ZIP III genes caused by mutations in the miRNA binding sites and increased expression are easier to detect than recessive mutations.In Arabidopsis,ectopic expression of the HD-ZIP III genes PHB,REV,and PHV results in the formation of radially symmetrical leaves[40,41,53].Ectopic expression of LF1 in rice causes lateral florets to develop at the axil of the sterile lemma[3].The changes observed in these species are similar to those observed here in wheat,suggesting that the disruption of the miRNA binding site and the resulting increased expression is a common mechanism underlying these phenotypes[3,54].This hypothesis is supported by the similar genotypes generated in this study by the overexpression of the miRNA-resistant allele of HB-D2 in wheat(Fig.7).In rice,overexpression of miR166-resistant versions of HD-ZIP III genes resulted in morphological changes whereas over-expression of the WT versions did not[3,13].

Fig.6.rHb-A2 and rHb-B2 EMS mutant plants showing phenotypic changes relative to plants homozygous for the WT allele.(A)Visual phenotypes of plants homozygous for the WT(Hb-A2)and mutant(rHb-A2)alleles in the line Cadenza1761 and sequencing chromatograms showing the G580A polymorphism.(B)Visual phenotypes of plants homozygous for the WT(Hb-B2)and mutant(rHb-B2)alleles in the line Cadenza0269 and sequencing chromatograms showing the C589T polymorphism.The mutated nucleotides are indicated by blue arrows.Paired spikelets are indicated by red arrows.

4.3.Effect of different mutations and environmental conditions on HB2

We observed plants carrying rHb-D2,rHb-A2,and rHb-B2 mutant plants displayed varying degrees of phenotypic effects(Figs.4,6).The rHb-D2 mutant showed the strongest phenotypic effect followed closely by rHb-A2,whereas rHb-B2 displayed the weakest effect.This gradient of phenotypic effects may be the result of the combination of two different effects:the stronger disruptive effect of rHb-D2 and rHb-A2 on the predicted amino acid changes(smaller BLOSUM62 values)and the larger effect of rHb-D2 on the predicted minimum free energy of the interaction with the miRNA relative to the other two mutations(Fig.S10).The different locations of the mutations may also contribute to the different intensity of the phenotypes.Similarly to the Rld1-O mutant in maize[4],the rHb-D2 and rHb-A2 mutations disrupt pairing at the 3′end of miRNA165/166(Figs.4,S9),whereas rHb-B2 affects a more central part of the miRNA binding site(Fig.S9).

Several HD-ZIP genes are known to be associated with changes in heading date or flowering time,such as sgd2[55]and Roc4[56].Our data demonstrate that the rHb-D2 mutation was also associated with 8.8 days delay in heading time when the plants were grown in a growth chamber at 10 °C with a 12 h photoperiod(Fig.S3).However,the effect on heading time was reduced when the plants were grown under a longer photoperiod and at higher temperature.By contrast,the intensity of some morphological defects(e.g.,rolled leaf and dwarfism)observed in the plants carrying the resistant rHb-D2 allele were milder when plants were grown at 10 °C/12 h photoperiod than when they were grown at higher temperatures and longer photoperiods(Fig.S3).These results suggest that photoperiod or/and temperature may modulate some of the effects of miRNA165/166 or its target genes,a result not reported in previous studies.

Fig.7.Phenotypic effects of rHb-D2 overexpression in transgenic plants of the wheat cultivar‘Fielder’.(A)Phenotypes of the T0 transgenic plants mT0-01,mT0-05,and T0-07 overexpressing rHb-D2 compared to the non-transgenic‘Fielder’control.(B)Representative images of leaves from the mT0-01 transgenic and control plants.(C)Representative images of spikes from the mT0-09 transgenic and control plants.Red arrows indicate the paired spikelets.(D)Phenotypes of the T1 transgenic families mT1-09 and mT1-12 overexpressing rHb-D2 compared to the non-transgenic‘Fielder’control.

Although the changes in heading time and the increases in the number of secondary spikelets(paired spikelets)could be beneficial for wheat improvement,the multiple pleiotropic effects of the miRNA165/166 resistant rHb2 alleles preclude their direct use in agriculture.The study of the targets of this transcription factor may help to dissect these different effects,and separate beneficial and detrimental effects.The validation of HB-D2 as the causal gene affecting the development of leaves,spikes,roots,and stems is a first step in that direction and contributes to a better understanding of the regulation of morphological traits in wheat.

CRediT authorship contribution statement

Dengji Jiang:Data curation,Investigation,Formal analysis,Visualization,Validation,Writing–original draft.Lei Hua:Data curation,Investigation,Formal analysis,Visualization,Validation,Writing–original draft.Chaozhong Zhang:Data curation,Investigation,Formal analysis,Visualization,Validation,Writing–original draft.Hongna Li:Software,Formal analysis,Visualization.Zheng Wang:Resources.Jian Li:Resources.Guiping Wang:Software,Formal analysis,Visualization,Resources.Rui Song:Software,Formal analysis,Visualization.Tao Shen:Software,Formal analysis,Visualization.Hongyu Li:Software,Formal analysis,Visualization.Shengsheng Bai:Software,Formal analysis,Visualization.Yanna Liu:Software,Formal analysis,Visualization.Jian Wang:Software,Formal analysis,Visualization.Hao Li:Resources.Jorge Dubcovsky:Resources,Funding acquisition,Writing–review&editing.Shisheng Chen:Project administration,Methodology,Funding acquisition,Writing–original draft,Writing–review & editing.

Data availability

Raw sequencing data generated in this project has been deposited in NCBI’s Sequence Read Archive(https://www.ncbi.nlm.nih.-gov/sra,Bioproject PRJNA826248).

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

We thank Juan Debernardi(UCD)for his suggestions on the analyses of the miRNA data.Work at SC lab was supported by the Provincial Natural Science Foundation of Shandong(ZR2021MC056 and ZR2021ZD30)and the Open Project Funding of the State Key Laboratory of Crop Stress Adaptation and Improvement.Research at JD lab was funded by Competitive Grant 2022-68013-36439(WheatCAP)from the USDA National Institute of Food and Agriculture.

Appendix A.Supplementary data

Supplementary data for this article can be found online at https://doi.org/10.1016/j.cj.2022.05.002.