Xinghi Yng,Boxun Nong,Cn Chen,Junrui Wng,Xiuzhong Xi,Zongqiong Zhng,Yu Wei,Yu Zeng,Rui Feng,Ynyn Wu,Hui Guo,Hifeng Yn,Yunto Ling,Shuhui Ling,Yong Yn,Dnting Li,*,Guofu Deng,*

a Guangxi Key Laboratory of Rice Genetics and Breeding,Rice Research Institute,Guangxi Academy of Agricultural Sciences,Nanning 530007,Guangxi,China

b Guangxi Key Laboratory for Polysaccharide Materials and Modifications,School of Marine Sciences and Biotechnology,Guangxi University for Nationalities,Nanning 530006,Guangxi,China

c Biotechnology Research Institute,Guangxi Academy of Agricultural Sciences,Nanning 530007,Guangxi,China

d Sugarcane Research Institute,Guangxi Academy of Agricultural Sciences,Nanning 530007,Guangxi,China

e Microbiology Research Institute,Guangxi Academy of Agricultural Sciences,Nanning 530007,Guangxi,China

Keywords:Rice OsNPF3.1 Functional analysis Evolutionary analysis Natural variation

ABSTRACT The overuse of nitrogen(N)fertilizer in fields has increased production costs and raised environmental concerns.Increasing the N use efficiency(NUE)of rice varieties is crucial for sustainable agriculture.Here we report the cloning and characterization of OsNPF3.1,a gene that controls rice NUE.An amino acid mutation in the OsNPF3.1 coding region caused different NUEs in wild and cultivated rice.OsNPF3.1,which is expressed mainly in the aerial parts of rice,also affects rice plant height,heading date,and thousand-grain weight.The OsNPF3.1 protein is located in the plasma membrane.When OsNPF3.1 was subjected to artificial selection,two naturally varying loci were associated with NUE,of which OsNPF3.1Chr6_8741040 differed between indica and japonica rice.OsNPF3.1 can be used as a new target gene for breeding rice varieties with high NUE.

1.Introduction

N is the mineral element most essential to plants and a limiting factor for plant growth and development.In rice production,massive use of N fertilizer has been used to achieve high yield[1].However,overuse of N not only increases production costs,but pollutes the agricultural environment[2,3].Thus,there is an urgent need to improve rice N-use efficiency(NUE).

Nitrate(NO3-)and ammonium(NH4+)are the main N forms that plants obtain from the soil[4].Nitrate absorption and transport proteins are mainly divided into two categories:NRT1s and NRT2s.NRT1s are composed of low-affinity nitrate transporters and are NITRATE TRANSPORTER 1/PEPTIDE TRANSPORTER family(NPF)members.OsNRT1 was cloned in rice using the homologous sequence of Arabidopsis AtNRT1.1[5].OsNRT1 has two different splicing forms:OsNRT1.1a and OsNRT1.1b.OsNRT1.1a functions only under high-N conditions,whereas OsNRT1.1b can also promote N uptake under low N conditions[6].Three paralogous genes were named OsNRT1.1A,OsNRT1.1B,and OsNRT1.1C based on their amino acid sequence similarity to AtNRT1.1[7].Overexpressing OsNRT1.1A increased NUE[8].A single-base variation in NRT1.1B accounted for the difference in NUE between indica and japonica subspecies[9].The NADH/NADPH-dependent nitrate reductase gene OsNR2,which interacts with OsNRT1.1B,promotes nitrate uptake in indica rice[10].OsNPF6.1 is a dual-affinity nitrate transporter,and the transcription factor OsNAC42 activated OsNPF6.1,thereby increasing rice NUE[11].

Most NRT2s are high-affinity nitrate transporters.Among the members of the NRT2s,OsNRT2.1,OsNRT2.2,and OsNRT2.3a need the assistance of OsNAR2.1 to complete the transport of nitrate[12].OsNAR2.1 interaction with OsNIT1 and OsNIT2 affected the absorption of nitrate and ammonium in rice roots[13].OsNRT2.3 produced two transcripts:OsNRT2.3a and OsNRT2.3b[12],of which OsNRT2.3a functions in long-distance transportation of nitrate from roots to aerial parts[14],and OsNRT2.3b increased the pH buffering capacity of rice[15].In contrast to that of all other NRT2 genes,the expression of OsNRT2.4 responds to changes in auxin level[16].

Ammonium is the main form of N source for rice.Rice has 12 ammonium transporters,which are divided into five subfamilies:OsAMT1–5[17].OsAMT1;1,OsAMT1;2,and OsAMT1;3 all have ammonium transporter activities.OsAMT1;1 is constitutively expressed in stems and roots,OsAMT1;2 is expressed specifically in roots,and OsAMT1;3 is expressed specifically in stems[18].OsAMT2.1 is constitutively expressed in roots and overground parts,whereas the expression of OsAMT3.1 is relatively weak[19].

Regulatory factors also affect NUE in rice.DEP1 interacts with RGA1 and RGB1 to regulate the N response of rice[20].Mutation of ARE1 delayed rice plant senescence and increased its tolerance for N starvation[21].The growth regulator GRF4 is a positive regulator of N metabolism that promoted N absorption,transportation,and assimilation in rice[22].NGR5 is a positive regulator of rice growth and development in response to N,and the highlevel accumulation of NGR5 and GRF4 proteins promoted the absorption and utilization of N fertilizer[23].The GRF4-MYB61 regulatory module integrates pathways controlling cellulose synthesis and NUE[24].Most of the nitrate and ammonium transporters identified to date are responsible for N absorption and transport in the roots,and little is known about the transport and utilization of N in the aerial parts[9,14,25–27],where N exerts its biological function and affects economic yield.

The aim of the present study was to fine map the nitrogen use efficiency gene OsNPF3.1 using a near-isogenic line,and to validate the biological function of OsNPF3.1.The superior allele of OsNPF3.1 was mined by natural population,and will provide genetic resources for cultivating resource-efficient rice varieties that are low-cost,high-yield,and environmentally friendly.

2.Materials and methods

2.1.Plant materials

The NUE of NIL19(GH998)is higher than that of NIL15(NIL-13B4)(Fig.S1)[28],both of which were planted in a paddy field.NIL15 and NIL19 were used for transgenic experiments.An F2progeny was generated by crossing the two lines.A major QTL qNUE6 was mapped in 266.5-kb region on chromosome 6 and qNUE6-interval heterozygous plants were screened with insertion/deletion(InDel)markers ID10 and ID22.2300 F3plants were generated.Seeds were planted in the experimental field of the Rice Research Institute of Guangxi Academy of Agricultural Sciences in March 2018(Nanning,Guangxi,China,22.85°N,108.26°E).Urea fertilizer was applied at 391.3 kg ha-1,potassium chloride at 319.4 kg ha-1,and phosphorus pentoxide at 750 kg ha-1.

2.2.Trait measurement

The heading dates of progenies,gene-edited lines,and RNAi plants were recorded at heading stage.After the beginning of the ear,the group is visually inspected.When the group have 50% of the rice ear heading,it was recorded as the heading stage.Plant height was measured at maturity stage.Thousand-grain weight was measured with an SC-G automatic seed test analysis instrument(Wseen,Hangzhou,Zhejiang,China)after the grain was dried to 12.5% water content.The phenotypes of 154 rice varieties from all around the world were measured as reported by Tang et al.[11].Rice NUE was determined as previously described[28,29].

2.3.Fine mapping

In the previous study,we mapped qNUE6 to a 266.5-kb interval.Using resequencing data,we developed 29,519 InDel markers on chromosome 6 in rice.These markers were used to genotype the F3population,and recombinant plants were identified phenotype and verified genotype in the progeny population.We identified recombination sites based on phenotypic and genotypic analysis results,and the recombination site was used to locate qNUE6.

2.4.Gene annotation and Sanger sequencing

A rice reference genome Oryza sativa v7.0(https://rice.uga.edu/)was used to annotate predicted genes in the candidate genomic region.Using NIL15 and NIL19 DNA as templates,LOC_Os06g15370 was PCR-amplified.The PCR products were purified and sequenced and mutation sites were characterized.

2.5.Synteny analysis

The synteny analysis homologous genes was performed around the rice OsNPF3.1 and Arabidopsis AtNPF3.1 loci.An image of the homologous gene collinearity was acquired with CoGE(https://genomevolution.org/coge/GEvo.pl).

2.6.RT-qPCR

Plant RNA preparation(TIANGEN,Beijing,China)was used to extract RNA from plant tissues.The reaction system and reagents were as in the previous study[30].Actin3 was used as an internal reference gene[20].The qPCR primers are listed in Table S1.qPCR was performed with a CFX96 fluorescent quantitative PCR instrument(BIO-RAD,Hercules,CA,USA),and the 2-ΔΔCTmethod[31]was used to calculate the gene expression level.

2.7.Gene editing

CRISPR-GE software[32]was used to design targets in the LOC_Os06g15370 coding region.According to the parameters,we chose two targets:sequence 1,GGCTGCTGCTGAACGCGCTGGG(+),and sequence 2,CGTCGGGGAGCCCCTTCACGCGG(+).The sgRNA(single guide RNA)sequences were inserted into the pYLCRISPR/Cas9 vector.The validated plasmid vector was transformed into Agrobacterium EHA105 competent cells,which were used to infiltrate and transform the rice line NIL19.Transgenic positive plants were identified by hygromycin screening and verified by sequencing the target region.The edited lines with amino acid mutations were transplanted into the greenhouse,and homozygous knockout mutants were identified by sequencing.

2.8.RNA interference

Two pairs of primers were designed for the sense and antisense fragments of rice LOC_Os06g15370 with Premier 5.0 software(https://www.premierbiosoft.com/primerdesign/index.html).The primer sequences are shown in Table S1,as are the primer sequences of the loop fragment between the sense and antisense fragments.RNAprep pure plant total RNA extraction kit(TIANGEN,Beijing,China)was used to extract total RNA from NIL19 and reverse transcription to obtain cDNA.The sense fragment(219 bp)of LOC_Os06g15370 was cloned by using NIL19 cDNA as the template and P1 as the primer,the antisense fragment of LOC_Os06g15370(219 bp)was cloned with primer P2,and the loop fragment(200 bp)was amplified using the original vector pBWA(V)HS-GH998(modified from pCAMBIA1300)as the template and P3 as the primer.The three amplified fragments were combined in the same system for purification and recovery.The recovered fragments were digested and ligated with the pBWA(V)HSGH998 original vector.The reaction system(20 μL)contained pBWA(V)HS-GH998 vector 4 μL,purified fragment 4 μL,BsaI/Eco31I restriction enzyme 1 μL,T4 DNA ligase 1 μL,buffer 2 μL,and ddH2O 8 μL.The reaction was incubated at 37 °C for 20 min,followed by 37 °C for 10 min and 20 °C for 10 min for 5 cycles(after each cycle,the reaction was held at 37°C for 20 min),and finally at 80°C for 5 min.The constructed plasmid was verified by sequencing and used for genetic transformation of NIL19.The detailed procedures are described in[33].

2.9.GUS(β-glucuronidase)staining

The target fragment with a full length of 2282 bp was obtained by PCR amplification using heat-resistant DNA polymerase and rice NIL19 DNA as the template. The PCR product was digested with AarI restriction enzyme. The enzyme digestion system containing PCR product 4 μL, AarI 1 μL, buffer 2 μL, and ddH2O 13 μL was incubated at 37 °C for 1 h. The vector pBWA(V)HG-GH998 (modified from pCAMBIA1300) was digested with BsaI/Eco31I. The vector digestion system containing vector pBWA(V)HG-GH998 4 μL, Bsa I/Eco31I 1 μL, buffer 2 μL, and ddH2O 1 3 μL was incubated at 37 °C for1 h .The digested PCR product and the digested vector were purified and recovered together, and then ligated with T4 DNA ligase. The ligation system containing PCR product and digested vector mixture 2.5 μL, T4 DNA ligase 1 μL, buffer 1 μL, and ddH2O 5.5 μLwas incubated at 20 °C for 1 h. The constructed plasmid was verified by sequencing and transformed NIL19 with Agrobacterium EHA105 [33].

GUS staining of rice tissues followed the manufacturer’s instructions for GUS staining solution(Leagene,Beijing,China).

2.10.Subcellular co-localization

The target gene was amplified with primer OsNPF3.1SL(Table S1),and the PCR product was loaded on a 1% agarose gel for separation by electrophoresis.The gel band containing the target gene was excised and transferred to a 2.0 mL centrifuge tube.The target fragment was recovered with a gel extraction kit(Omega Bio-tek,Norcross,GA,USA)and stored at-20 °C for later use.The plasmid pBI211-GFP was digested with SacI and XbaI.After digestion,the product was detected by agarose gel electrophoresis.The gel band corresponding to the large vector fragment was cut out and recovered.The target DNA fragment was ligated to the vector and transformed into DH5a competent cells.Colony PCR was used to identify positive clones,and Plasmid Mini Kit I(Omega Bio-tek)was used to extract the DNA from pBI221-GFP-OsNPF3.1 and pBI221-GFP-osnpf3.1 clones with correct sequences.Preparation of protoplasts from Nipponbare seedlings,and the prepared plasmid was added to the rice protoplast,and the fluorescent plant material was examined under a Olympus FV1000 confocal laser scanning microscope(Olympus,Tokyo,Japan),excited at 488 nm with an argon-ion laser,and separated the different wavelengths of light with a 545-nm spectroscope.GFP was detected at 505–530 m,and chlorophyll autofluorescence signal was detected above 560 m.The empty vector pBI221-GFP was used to transform rice protoplasts as a control.

2.11.SNP calling and phylogenetic and genetic diversity analysis

Genome sequences of 250 rice accessions from previous studies were retrieved from 3K(https://snpseek.irri.org/),National Center for Gene Research,Chinese Academy of Sciences,Genetic Strain Stock Center,Plant Genetics Laboratory,National Institute of Genetics.A rice reference genome(Oryza sativa v7.0)was retrieved from the Rice Genome Annotation Project(https://rice.uga.edu/).The genome sequence data of 168 O.sativa indica and 59 O.sativa japonica accessions were retrieved from the National Center for Biotechnology Information(NCBI).The sequences of 23 O.rufipogon accessions were retrieved from the NCBI BioProject PRJEB19404(https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJEB19404;ERR2240123,ERR2240125,ERR2240126,ERR2245548-ERR2245557), PRJDA39855 (https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJDA39855;DRR000347,DRR000348)and PRJDB2009 (https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJDB2009;DRR001183–DRR001190).Raw reads were trimmed with Trimmomatic v0.39(https://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-Src-0.39.zip)with the following parameters:ILLUMINA-CLIP:2:30:7,SLIDING WINDOW:2:20 MINLEN:60;and then the reads were aligned to Oryza sativa v7.07 using BWA-MEM v0.7.17(https://sourceforge.net/projects/bio-bwa/).SNP calling and filtering were performed with GATK(The Genome Analysis Toolkit)v3.8[34](QD(Qual by Depth)＜2 or a quality score＜30 or FS(Fisher Strand)＞60 or MQ(Mapping Quality)＜40).The SNP(Single Nucleotide Polymorphisms)information from the 20 kb upstream and downstream of LOC_Os06g15370 was used to generate a maximum likelihood phylogeny using MEGA v10.13[35]with 1000 bootstrap replications.The online tool iTOL v5.6.3(https://itol.embl.de)was used to display the phylogenetic tree and its annotation.Nucleotide diversity(π),population differentiation,and Tajima’s D statistics of OsNPF3.1 gene were calculated with vcftools v0.1.15(https://sourceforge.net/projec ts/vcftools/files/)for indica,japonica and wild rice groups,with a 1000 nt sliding window and 500 nt steps.

2.12.Haplotype network

Variable sites of OsNPF3.1 in the 250 samples were used for allelic identification.Haplotype-frequency data were processed with DnaSP v6.12.03[36],and median-joining networks were visualized with PopART[37].

2.13.PARMS genotyping

The rice SNP-Seek Database(https://snpseek.irri.org/)was used to identify SNPs in LOC_Os06g15370,and information for 31 SNPs in 3024 varieties were obtained.The variation information of OsNPF3.1 was recorded.Penta-primer amplification refractory mutation(PARMS)[38]primer sets were designed for the two loci with the most frequent mutations.The primer sequences are presented in Table S1.The reaction system and three fluorescencedetection methods are described by Gao et al.[38].

2.14.Association analysis

The basic scenario of a compressed mixed linear model implemented in the R package Genomic Association and Prediction Integrated Tool(GAPIT)[39],was adopted for the association analysis of QTL(quantitative trait locus)-flanking markers with four traits for the 154 sequenced accessions.The genomic association analysis employed the mixed linear model(MLM)and the general linear model(GLM).A medium-throughput SNP genotyping by PARMS[38]was used to acquire SNPs.To minimize the possibility of type II errors in QTL detection,a relatively low threshold was adopted for QTL regions with supporting evidence from multiple traits.

2.15.Statistical analysis

CorelDRAW X8(https://www.coreldrawchina.com/)was used to plot fine mapping.Bar graphs and scatter diagrams were plotted using Graphpad Prism 8.(https://www.graphpad.com/guides/prism/8/user-guide/)Sanger sequencing data were analyzed using software Chromas 2.6.5(https://technelysium.com.au/wp/chromas/).ClustalX 2.1(https://www.clustal.org/clustal2/)was used to perform gene sequence comparison.

Fig.1.Map-based cloning of OsNPF3.1.(A)Schematic fine mapping of qNUE6.The vertical lines represent the sites of key InDel markers.Structure of the OsNPF3.1 gene.Blue boxes,gray boxes and black line represent exon,untranslated regions and intron,respectively.Chr.,chromosome.(B)Sequence comparison of OsNPF1 in NIL15 and NIL19.Yellow marker represents the mutation distinguishing NIL15 and NIL19.(C)Synteny analysis around rice OsNPF3.1 and Arabidopsis thaliana AtNPF3.1 with the image generated by CoGE(https://genomevolution.org/coge/GEvo.pl).

3.Results

3.1.Map-based cloning of rice NUE gene OsNPF3.1

A total of 2300 F3plants were developed using DNA markers to screen for plants with heterozygous qNUE6 segment.There are 53 recombinants between markers ID10 and ID22,and their progeny were planted.Twenty seeds were planted for each recombinant line,and the progenies were tested for recombination.We identified 17 plants with fragment intersection points marked as ID28 and ID41.Finally,we mapped qNUE6 in the 8,719,820–8,767,839 bp(48.02 kb)interval(Fig.1A).A rice reference genome Oryza sativa v7.0(https://rice.uga.edu/)was used to annotate candidate genome regions,and eight predicted genes were found(Table S2).

There were 20 variations(Fig.S2),among which the third-exon G＞A mutation caused a change of threonine to alanine(Fig.1A,B),and the mutation was located within the conserved domain(pfam00854)of LOC_Os06g15370.

The sequence similarity between LOC_Os06g15370 and Arabidopsis AT1G68570(AtNPF3.1)was 78.3%.There was synteny between LOC_Os06g15370 and AT1G68570(Fig.1C).David et al.[40]found that AtNPF3.1 expression was upregulated under limited N nutrition.

We identified LOC_Os06g15370 as a candidate gene of qNUE6.Given that Léran et al.[41]named the NPF family members,we named LOC_Os06g15370 as OsNPF3.1.

3.2.Function analysis of OsNPF3.1

Using NIL19 sequence as the template,we designed a singleguide RNA(sgRNA)targeting the PTR2 domain in the third exon of OsNPF3.1:target 1,GCGCTGCTGCTGAACGCGCTGGG;target 2,CGTCGGGGAGCCCCTTCACGCGG (Fig.2A).The sgRNAs were cloned into the pYLCRISPR/Cas9 vector.The constructed plasmid was amplified with primer Pbw2.After confirmation of the sgRNA sequence,the plasmid was transformed into competent Agrobacterium EHA105.NIL19 callus was transformed with Agrobacterium EHA105[33],and 11 positive plants with kanamycin resistance were obtained.

To identify the mutation of OsNPF3.1 in the edited lines,the target region in T0plants was amplified with primer PTR2 and then sequenced.Twelve edited lines carried mutations in the exon region,and there were 10 mutant genotypes(Table S3),each of which carried an amino acid frameshift.In the transgenic T1plants,the OsNPF3.1 target sites of these plants were sequenced and homozygous mutants were obtained.The offspring of these homozygous lines were used for subsequent experiments.

In the edited line GE1,a T base insertion at the third exon 414 resulted in a frameshift mutation,creating a stop codon(TGA)and premature translation termination(Figs.2B,S3).

GE1,NIL15,and NIL19 were then planted in a greenhouse.The N absorption and use efficiency(NAUE)differed between GE1(6.1%)and NIL19(33.4%)(Fig.2C).The plant height,heading date,and thousand-grain weight of GE1 and NIL19 also differed.There was no significant difference in NAUE,plant height,or heading date between GE1 and NIL15(Fig.2C,D,E,H).However,the thousandgrain weights of GE1 and NIL15 were different(P=9.86×10-5;Fig.2D,E).

After the construction of the OsNPF3.1 RNA interference(RNAi)vector,we transformed the NIL19 callus with pCambia1300::RNAi-OsNPF3.1 to obtain transgenic plants.For the transgenic lines,we extracted RNA from the leaf sheath at the 4-leaf stage and measured the expression level of OsNPF3.1 using RT-qPCR.The expression of OsNPF3.1 in transgenic line RI1 was significantly decreased(Fig.2I).NAUE of RI1(6.4%)and NIL19(33.4%)were significantly different(Fig.2C).The plant height,heading date,and thousand-grain weight of RI1 were consistent with those of NIL15(Fig.2F–H).The thousand-grain weight of RI1 and NIL15 differed(P=0.0208;Fig.2F,G).

Fig.2.Functional identification of OsNPF3.1.(A)Schematic representation of OsNPF3.1 locus,red,blue and purple boxes represent 5′-untranslated region(5′UTR),exons and 3′UTR respectively;black lines represent introns;the single guide RNA(sgRNA)with protospacer adjacent motif(PAM)sequences are at the bottom.(B)CRISPR/Cas9-induced mutations in T2 plant.Yellow box represents the nucleotide mutation locus in GE1.(C)Determination of N absorption and use efficiency NIL15,NIL19,GE1,and RI1.(D)Comparison of plant heights,heading dates and thousand-grain weights among NIL15,GE1,and NIL19.(E)Statistical analysis of plant height,heading date,and thousandgrain weight in NIL15,NIL19,and GE1.(F)Comparison of plant heights,heading dates and thousand-grain weights among NIL15,RI1 and NIL19.(G)Statistical analysis of plant height,heading date,thousand-grain weight in NIL15,NIL19,and GE1.(H)Flowering stage NIL15,GE1,NIL19,and RI1 plants grown under 180 kg ha-1 N in greenhouse.(I)Expression of OsNPF3.1 in NIL15,NIL19 and RI1.(J)Fresh weight of GE1 and NIL19 at 30,40,50,60,70,80,90,100,and 110 d.(K)Dry weight of GE1 and NIL19 at 30,40,50,60,70,80,90,100,and 110 d.Values are means±SD(n=10).P-values from Student’s t-test are indicated.Scale bars,10 cm.*,P＜0.05;**,P＜0.01;***,P＜0.001.

We also investigated the dynamic changes in fresh and dry weight of GE1 and NIL19 throughout the growth and development stages,and the difference in biomass production was significant.At day 100,the fresh weights of GE1 and NIL19 were 80.0 g and 169.1 g and the dry weights was 20.0 g and 41.1 g,with a very significant difference between them(P＜0.001;Fig.2J,K).

3.3.The expression pattern of OsNPF3.1

Sampling started from the 20th day after germination of the seeds,and subsequently at 40,60,80,100,and 110 d.The RNAs of NIL15 and NIL19 were extracted as follows:roots,leaf sheaths and leaves at seedling and tillering stage;roots,culms and leaves at booting stage;and roots,culms,leaves and panicles at flowering stage.At seedling and tillering stages,OsNPF3.1 was expressed mainly in culms and leaves(Fig.3A,B);at booting stage,OsNPF3.1 was expressed mainly in culms(Fig.3C,D);and at flowering stage,OsNPF3.1 was expressed mainly in culms and panicles(Fig.3E,F).

We also transformed NIL19 with vector pOsNPF3.1::GUS,and performed GUS staining of the transgenic plants using NIL19 as a control.At seeding stage,nodes were stained with GUS staining solution(Fig.3G).At tillering and booting stage,leaves,leaf sheaths,and culms were stained(Fig.3H–J).However,OsNPF3.1 expression was not observed in roots(Fig.3G,K;At flowering stage,young panicles were stained(Fig.3L).

Fig.3.Expression profiles of OsNPF3.1.(A–F)RT-qPCR results showing the expression of OsNPF3.1 in NIL15 and NIL19 during the growth and development.(G–L)Representative images of GUS expression driven by the native OsNPF3.1 promoter,and using NIL19 as control.Seeding(G);leaf(H);culm(I);leaf sheath(J);root at booting stage(K);panicle(L).Values are means±SD(n=3).P-values from Student’s t-test are indicated.WT,NIL19;TP,transgenic plant.Scale bars,1.5 cm in(G),0.5 cm in(H–J),2 cm in(K),and 6 cm in(L).Arrows indicate the differential parts of GUS staining.

3.4.Subcellular localization of OsNPF3.1-GFP fusion protein

Total RNA was isolated from plants of both OsNPF3.1(NIL19)and osnpf3.1(NIL15).Reverse transcription was performed with oligo(dT)primers to obtain cDNAs.The 35S:OsNPF3.1-GFP and 35S:osnpf3.1-GFP vectors were transiently transformed into rice protoplasts to express OsNPF3.1-GFP and osnpf3.1-GFP fusion proteins.The subcellular localization of the fusion protein was observed under a confocal laser scanning microscope.The localization difference between OsNPF3.1-GFP and osnpf3.1-GFP was on the cell membrane(Fig.4A–C).

3.5.Evolutionary analysis of the OsNPF3.1 locus

A 40,000-bp genomic region covering OsNPF3.1 was acquired from 250 rice accessions in 3K and NCBI database,including 168 from indica,59 from japonica,and 23 from O.rufipogon(Table S4).The 40,000-bp region contained a 20,000-bp upstream region,a 4590 bp coding region,and a 20,000-bp downstream region.A 1000-bp sliding window with a step size of 500 bp was used to calculate the overall nucleotide diversity(π)and Fst within or between the three rice groups(indica,japonica and O.rufipogon).As shown in Fig.5A,compared to O.rufipogon,both indica and japonica showed greatly reduced genetic diversity,suggesting that OsNPF3.1 resides in a selective-sweep region and acts as a breeding target.The level of genetic difference between indica and O.rufipogon(Fst=0.228)was greater than that between japonica and O.rufipogon(Fst=0.196),and there was an obvious genetic distinction between indica and japonica,which had a relatively high level of population differentiation(Fst=0.335)(Fig.5B).Correspondingly,selection on the OsNPF3.1 gene and departure from neutrality were also tested at segregating nucleotide sites using Tajima’s D in each group.Tajima’s D value for the indica group was 2.204 and was significantly different from both of the other two groups,japonica(-0.276)and O.rufipogon(-0.172).These results indicated that the OsNPF3.1 gene in the indica group was strongly selected during rice domestication,whereas OsNPF3.1 has been subjected to purifying selection in japonica and O.rufipogon.We then performed phylogenetic analysis.Nearly all japonica and most indica accessions with several wild rice clustered together(Fig.5C),while some other indica accessions separated into a clade with their wild ancestors,indicating that OsNPF3.1 was selected independently,and they fixed different haplotypes respectively.To assess the genetic relationships among taxa we calculated the numbers of haplotypes and constructed a phylogenetic tree.Mutation sites of OsNPF3.1 in the 250 samples were identified,and only three haplotypes were found(Fig.5D).Consistent with the phylogenetic tree,indica accessions had two distinct haplotypes(Hap 1 and Hap 2),one of which contained almost all japonica samples.Each haplotype had its own wild population.The samples in haplotype Hap 3 were collected from india and Myanmar and had heterozygous sites.We speculate that OsNPF3.1 allele with the functional mutation originated in and diverged from O.rufipogon during indica rice domestication,leading to higher nitrate-use efficiency than japonica,but more in-depth research should be performed.

Fig.4.Subcellular localization of OsNPF3.1-GFP fusion protein in rice protoplasts.Confocal scanning(GFP),plasma membrane(PM)-mCherry signal,bright field(BF)and merged micrographs of rice protoplast cells transformed with a 35S:OsNPF3.1-GFP construct(A),35S:osnpf3.1-gfp construct(B),and 35S:GFP vector control(C).Scale bars,10 μm.

3.6.OsNPF3.1 natural variations affects rice NUE

The OsNPF3.1 allele variations were retrieved from the Rice SNPSeek Database in the 3K database,and a total of 31 SNPs were obtained(Table S5).The PARMS method was used to genotype the 154 rice varieties,and 18 SNPs could be available for subsequent analysis.We performed association analysis by genotypes(Table S6)and phenotypes(Table S7)of the 131 indica and 23 japonica rice varieties(Fig.S4).There were two SNPs that might be associated with NUE under the GLM and MLM models(Fig.6A,B):OsNPF3.1Chr6_8741040and OsNPF3.1Chr6_8742153.The genotype of the Nipponbare genome at OsNPF3.1Chr6_8741040is GG.The GG genotype was represented by 19 japonica varieties and the AA genotype by 127 indica varieties(Fig.6C).The GG is a natural variant that is insensitive to N.The genotype of Nipponbare genome in OsNPF3.1Chr6_8742153was TT.There were eight rice varieties with the TT genotype and 146 with the GG genotype(Fig.6D).The TT genotype is a natural variant that is insensitive to N.We used the 3K data to analyze the variation of OsNPF3.1Chr6_8741040.The A type accounted for 57.9%,which mainly distributed in indica rice(ind1A,ind1B,ind2,ind3,indx);G type accounted for 40.6%,distributed mainly in japonica rice(japx,subtrop,temp,trop)and aus.At OsNPF3.1Chr6_8742153,the G type accounted for 92.5% and was widely distributed in 13 subpopulations.The T type accounted for 6.9%,distributed mainly in the temp type(https://snpseek.irri.org/).

4.Discussion

4.1.OsNPF3.1 is an NRT1/PTR family gene and affects rice nitrogen use efficiency

In rice,many genes/QTL related to NUE have been mapped on chromosome 6(Table S8),but only a few genes were cloned.Overexpression of OsPTR9 increases ammonium uptake in rice,promotes lateral root formation,and increases yield[42].The nitrogen-responsive negative regulator OsLBD37 and OsLBD39 can directly bind to the OsTCP19 promoter to regulate the expression of OsTCP19,and as a transcription factor,OsTCP19 can inhibit the expression of DLT,a key component in the brassinolide signaling pathway,thereby regulating rice tillering[43].In this study,OsNPF3.1 was confirmed to increase nitrogen use efficiency in rice.

OsNPF3.1 belongs to the NPF members[41].In rice,OsNPF2.2[26],OsNPF2.4[25],OsNPF4.1[27],OsNPF4.5[44],OsNPF6.1[11],OsNPF6.3[8],OsNPF6.5[9],OsNPF7.1[45],OsNPF7.2[46],OsNPF7.3[47],OsNPF7.4[45],OsNPF7.7[48],OsNPF8.9[5],and OsNPF8.20[58]are NUE-related genes.The functional study of these NPF genes will help to analyze the genetic regulation network of rice NUE,which will greatly promote the genetic improvement of rice NUE.

Fig.5.Evolutionary analysis of OsNPF3.1.(A)DNA sequence diversity of the genomic region surrounding OsNPF3.1 in three groups.The red,blue and green lines indicate site nucleotide diversity(π)for Oryza sativa indica,O.sativa japonica and O.rufipogon accessions,respectively.Y-axis,average π value;x-axis,Chr.6:chromosome 6.The position of OsNPF3.1 is as indicated.(B)The diversity(π)and genetic distance(Fst)based on OsNPF3.1 and flanking regions(～40 kb)across the groups,where color indicates phylogenetic group;radius of π indicates genetic diversity;and dashed line length indicates Fst value between two groups.(C)Phylogenetic tree of OsNPF3.1 and flanking regions(～40 kb)from 250 accessions.(D)Haplotype network of OsNPF3.1.Each haplotype group is represented by a circle,and circle size represents the number of accessions within the haplotype.Blue,O.sativa indica;red,O.sativa japonica;green,O.rufipogon.

4.2.OsNPF3.1 affects heading date of rice

The functional analysis of OsNPF3.1 showed that this gene affected rice heading date.Previous study[28]have shown that the near isogenic lines NIL15 and NIL19 differ only on chromosomes 6,8,9,and 10.The difference on chromosome 6 is between 4.3 and 22.1 Mb,and the critical gene for rice heading stage Hd1 is located 9336359–9338643 bp on chromosome.It is thus possible that the difference in growth period between NIL15 and NIL19 was caused by the Hd1 allele mutation.Yano et al.showed that the second exon of the Hd1 allele in Kasalath had a 2-bp deletion,resulting in an early termination[49].The Hd1 protein of Kasalath was shorter than that in Nipponbare because of the lack of a Cterminal region.Sanger sequencing was used to sequence the Hd1 genes of NIL15,NIL19,Kasalath,and Nipponbare,and there was no difference between NIL15 and NIL19 in the functional sites on the Hd1 second exon(Fig.S5).Therefore,the difference in growth period between near-isogenic lines is caused not by Hd1 but by OsNPF3.1.

4.3.OsNPF3.1 is expressed in the aerial parts of rice

The above-ground part of the plant is the place where the amount of nitrogen exerting its biological function is most important[17].Many genes involved in N are expressed in root[17,43],but only a few genes have been found to be involved in the transport of N in the aerial part of rice[17].OsNPF2.4 is a pH-dependent,low-affinity nitrate transporter gene,and is expressed mainly in rice epidermis,xylem parenchyma,and phloem sieve tube companion cells[25].OsNPF4.1 is highly expressed in the phloem of young panicles and affects panicle length[27].The expression of OsNRT1.1B is induced by nitrate and is expressed mainly in root hairs,epidermis,and vascular tissue[9].In the present study,GUS staining showed that OsNPF3.1 was expressed mainly in rice culms and young panicles,but not in roots.This result suggests that OsNPF3.1 may affect in N transport in the aboveground part of rice.

4.4.OsNPF3.1 underwent artificial selection

The analysis of 250 materials showed that the nucleotide polymorphism of OsNPF3.1 in japonica rice was significantly reduced.Cultivated rice in Asia is divided into two main subspecies,japonica and indica rice,which have different characteristics in morphology,development,and physiology.There are differences between indica and japonica subspecies in nitrate absorption and transport.Indica rice has a higher nitrate absorption and transport efficiency than japonica rice.Hu et al.[9]found that OsNRT1.1B had nitrate transport activity under both high and low N conditions,and OsNRT1.1B activity was higher in indica than in japonica rice.Gao et al.[10]found that OsNR2 encoded nitrate reductase,and its allelic mutation accounted for the difference in nitrate assimilation and NUE between the indica and japonica subspecies.OsNR2 in indica showed stronger nitrate reductase activity.

Fig.6.Superior-genotype analysis of OsNPF3.1.Measurement of the four traits of OsNPF3.1Chr6_8741040(A)and OsNPF3.1Chr6_8742153(B).X-axis,genotype of SNP;Y-axis,trait value;Genotyping of 154 rice varieties using PARMS markers OsNPF3.1Chr6_8741040(C)and OsNPF3.1Chr6_8741040(D).X-axis,FAM fluorescence signal value;Y-axis,HEX fluorescence signal value.

Although some genes related to nitrogen uptake,transport and assimilation have been identified in rice Several N-use genes have been identified[50,51],their evolutionary history remains to be clarified.The NUE for indica varieties is 30%–40% higher than that for japonica varieties.Only three genes,NRT1.1B,ARE1,and NR2 have so far been shown to have functional divergence in NUE between the two rice subspecies[9,10,21].Hu et al.[9]suggested that NRT1.1B was subjected to artificial selection during indica domestication,leading to higher NUE.Assessment of NRT1.1B orthologs in the Oryza genus showed that NRT1.1B-indica is a later-derived allele.Hap 1 in indica retained only one genotype(T)from its direct ancestor Oryza rufipogon-I,which has two genotypes(C and T),whereas SNP1 in japonica retained the only genotype(C)from its direct ancestor O.rufipogon-III,suggesting that NRT1.1B-indica has undergone directional selection.Moreover,NRT1.1B-indica had been demonstrated that a single nucleotide difference in NRT1.1B strongly influenced both NUE and yield.For the OsNR2 gene,haplotype analysis based on three SNPs indicated that OsNR2 of 9311(O.sativa L.spp.indica)and Nipponbare(O.sativa L.spp.japonica)mainly represent allelic variation between indica and japonica.OsNR2-indica was grouped with the O.rufipogon alleles in the phylogenetic analysis,indicating greater divergence between OsNR2-japonica and the presumed ancestral alleles.As shown in the ka/ks analysis,OsNR2-japonica was under positive selection during indica–japonica differentiation,even as a weak allele in NR activity[10].A study[10]indicated that directional selection has driven the evolutionary divergence of the indica and japonica OsNR2 alleles.In our study,nucleotide diversity and neutrality analysis suggested that OsNPF3.1 might be undergoing directional selection.In particular,the neutrality test detected selection signatures for indica alleles(Tajima’s D=2.204 for indica),which is evidence of selective sweeps.The haplotype network and phylogenetic tree revealed the grouping of indica OsNPF3.1 s together with those from O.rufipogon accessions,suggesting greater divergence of japonica OsNPF3.1 from the presumed ancestral wild rice OsNPF3.1.Taken together,our results suggest that OsNPF3.1 was subjected to artificial selection during indica domestication,subsequently leading to higher NUE.By using the 3K database,we found that OsNPF3.1Chr6_8741040and OsNPF3.1Chr6_8742153were significantly associated with rice NUE. The OsNPF3.1Chr6_8742153locus was significantly different between indica and japonica rice,suggesting that the difference in this locus may account for the difference in NUE between indica and japonica rice.

5.Conclusions

OsNPF3.1 affects rice NUE,plant height,heading date,and thousand-grain weight.OsNPF3.1 is expressed in the aerial parts of rice,with highest expression in culms and young panicles.OsNPF3.1 is located on the cell membrane.OsNPF3.1 was strongly selected during rice domestication.Its natural variation affects the difference in NUE between indica and japonica rice.OsNPF3.1 is a critical gene for high NUE breeding in rice.

CRediT authorship contribution statement

Xinghai Yang:Conceptualization,Funding acquisition,Validation,Writing–original draft.Baoxuan Nong:Investigation,Data curation.Can Chen:Validation,Formal analysis.Junrui Wang:Software.Xiuzhong Xia:Formal analysis.Zongqiong Zhang:Project administration.Yu Wei:Resources.Yu Zeng:Investigation.Rui Feng:Methodology.Yanyan Wu:Formal analysis.Hui Guo:Validation.Haifeng Yan:Validation.Yuntao Liang:Methodology,Resources.Shuhui Liang:Investigation.Yong Yan:Conceptualization,Supervision.Danting Li:Funding acquisition,Writing–review & editing.Guofu Deng:Conceptualization,Supervision.

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

This work was supported by the National Natural Science Foundation of China(32060476 and 31860371),Guangxi Department of Science and Technology(AA22068087-4),Guangxi Natural Science Foundation of China(2015GXNSFAA139054,2018GXNSFAA138124,and 2020GXNSFAA259041),Guangxi Ministry of Science and Technology(AB21238009),Special Fund of Local Science and Technology Development for the Central Guidance(ZY21195034),and Guangxi Academy of Agricultural Sciences(2021JM04,2021JM49,2021YT030,QN-25,and QN-35).

Data availability

The RNA-seq data that support the findings of this study have been deposited into the National Center for Biotechnology Information(NCBI)Sequence Read Archive(SRA)with accession code PRJNA772895(SRA no.from SRR16493776 to SRR16493793).

Appendix A.Supplementary data

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