extfieldsofstudy: Medicine Biology

text:
Mapping quantitative trait loci for yield-related traits and predicting candidate genes for grain weight in maize

Quantitative trait loci (QTLs) mapped in different genetic populations are of great significance for marker-assisted breeding. In this study, an F2:3 population were developed from the crossing of two maize inbred lines SG-5 and SG-7 and applied to QTL mapping for seven yield-related traits. The seven traits included 100-kernel weight, ear length, ear diameter, cob diameter, kernel row number, ear weight, and grain weight per plant. Based on an ultra-high density linkage map, a total of thirty-three QTLs were detected for the seven studied traits with composite interval mapping (CIM) method, and fifty-four QTLs were indentified with genome-wide composite interval mapping (GCIM) methods. For these QTLs, Fourteen were both detected by CIM and GCIM methods. Besides, eight of the thirty QTLs detected by CIM were identical to those previously mapped using a F2 population (generating from the same cross as the mapping population in this study), and fifteen were identical to the reported QTLs in other recent studies. For the fifty-four QTLs detected by GCIM, five of them were consistent with the QTLs mapped in the F2 population of SG-5 × SG-7, and twenty one had been reported in other recent studies. The stable QTLs associated with grain weight were located on maize chromosomes 2, 5, 7, and 9. In addition, differentially expressed genes (DEGs) between SG-5 and SG-7 were obtained from the transcriptomic profiling of grain at different developmental stages and overlaid onto the stable QTLs intervals to predict candidate genes for grain weight in maize. In the physical intervals of confirmed QTLs qKW-7, qEW-9, qEW-10, qGWP-6, qGWP-8, qGWP-10, qGWP-11 and qGWP-12, there were 213 DEGs in total. Finally, eight genes were predicted as candidate genes for grain size/weight. In summary, the stable QTLs would be reliable and the candidate genes predicted would be benefit for maker assisted breeding or cloning.

different populations. For example, the trichome density trait was confirmed with four recombinant inbred lines (RIL) populations of A. thaliana 13 , the stay green trait was confirmed with two RILs populations of sorghum 14 , and the kernel length trait of barley was confirmed by different RILs population 15 .
Grain weight is one of the most important yield-related trait in crops. It is of great significance to clone the genes controlling grain weight and then to clarify their molecular genetic mechanism. Great achievements have been made in genes/QTLs identification and dissection for grain size and grain weight in many crops, such as rice [16][17][18][19][20] , soybean 21,22 , wheat 23,24 . Especially many genes related to grain weight or grain size in rice, including GS3 25 , GS5 16 , qGL3 26 , GW2 18 , GW8 27 , GS2 28 , qGW7/GL7 29 , have been successfully cloned. The evidence showed that the grain size trait was regulated by multiple signaling pathways, and the main pathway included ubiquitin-proteasome degradation pathway 30 , phytohormone signaling pathway and G protein independent pathway. Earlier studies showed that crop's yield was greatly influenced by grain size and grain weight. Although great achievements on the genes controlling grain size and grain weight have been made in maize in recent years 5,[31][32][33] , it was relatively low compared to rice and Arabidopsis thaliana. It is necessary for maize yield-related traits to confirm stable QTLs, to decrease functional gene number in stable QTL intervals and to predict candidate genes. All these works would provide basis for cloning functional genes and marker-assisted breeding.
The purposes of this study were: (1) to identify QTLs for yield-related traits with an F 2:3 population from SG5 × SG7; (2) to compared these QTLs detected in this study with the QTLs identified in other populations, including an F 2 population from the same parents as the F 2:3 used in this study and recent studies; (3) RNA-seq technology was applied to identify transcriptional variations between maize inbred lines SG5 and SG7 subjected to grain weight; (4) to identify the DEGs related to grain weight between SG5 and SG7 and to predict candidate genes.

Results
phenotype evaluation of the mapping population. The phenotypic data of seven yield traits, i.e., ear length (EAL), ear diameter (EAD), cob diameter (CD), kernel row number (KRN), ear weight (EW), grain weight per ear (GWP) and 100-kernel weight (KW), were collected from the F 2:3 mapping population in 2016. The mean values of the seven traits were shown in Table 1. The phenotypic values of the two parents SG5 and SG7 were different in all the seven traits. The seven yield traits all displayed in bell-shaped normal distribution (Supplemental Fig. 1). Table 2 listed the pearson's correlation coefficients and the significance tests between every two of the seven observed traits. The closes relation occurs between GWP and EW (0.975).
QtL analysis using a high-density linkage map. CIM and GCIM procedures were applied to identify QTLs associated with the seven yield traits based on a high-density linkage map. The map constructed in previous study 34 , had 3305 bin-markers. Using CIM, a total of thirty-three QTLs were detected for the seven yield traits. These QTLs were distributed on all of the 10 maize chromosomes. Among them, four QTLs controlled KW trait and were located on maize chromosomes 3, 7, 8, and 9; five related to EAL trait were located on chromosomes 1, 2, 3, and 5; four associated with EAD trait were located on chromosomes 1, 4, and 10; four controlling CD trait were located on chromosomes 1, 3, and 6; five controlled KRN trait and were located on chromosomes 2, 3, 4, and  9; four were associated with EW trait and located on 2, 5, 7, and 9; and seven controlling GWP trait were located on chromosomes 2, 3, 5, 7, and 9. QTL mapping result from CIM procedure is shown in Fig. 1 and .6%, respectively. The logarithm of odds (LOD) scores ranged from 3.0 (qKW-12) to 7.4 (qEAD1-1). For GCIM, a total of fifty-four QTLs were detected. Among these QTLs, fourteen were also detected by CIM (Table 4), three KW QTLs were located on chromosomes 7, 8 and 9, two EAD QTLs on chromosomes 1, and 2, nine CD QTLs on chromosomes 1, 2, 4, 5, 6, and 7, thirty-four KRN QTLs on over all chromosomes except chromosomes 6 and 9, two EW QTLs on chromosomes 2 and 7, and four GWP QTLs on chromosomes 5, 7, and 9. Their LOD scores ranged from 3.1 (qKW-7) to . The GCIM approach seems more powerful in detecting small effects QTLs, especially for early generation population. The related information is summarized in Table 4 .

Confirmation of QTLs in different generation materials.
It is essential to confirm QTLs before used in marker assisted breeding 35 . We compared the QTLs mapped with the F 2:3 population with those detected using other populations, including an F 2 generating the same parents as the mapping population in this study 34 . For the 33 QTLs mapped by CIM in this study ( Table 3), eight of them located on chromosomes 1, 2, 3, and 7 were consistent with those detected in the F 2 population (Fig. 2). There were fifteen major QTLs overlapped with the those identified in previous studies [36][37][38] (Fig. 2), and they were located on chromosomes 1, 2, 3, 5, 7 and 9. For the fifty-four QTLs detected by GCIM, five of them(qKW-7, qEAD1-1, qCD-3, qKRN2-1, qKRN-1) were also mapped in the F 2 population from the same parents, and twenty one were identical to the reported QTLs in other studies ( Table 4). The QTL detected in more than one mapping populations is considered to be stable QTL.
To sum up, for all the QTLs detected by CIM and GCIM methods, thirty of them were stable QTLs, and fourteen and twenty-two were confirmed by CIM and GCIM, respectively. The detailed information about confirmed QTLs were listed in Tables 3 and 4. The confirmed QTLs for grain weight include one major QTL for KW on chromosome 7 (qKW-7, mapped by CIM and GCIM), two major QTLs for EW on chromosome 7 and 9 respectively (qEW-9, mapped by CIM and GCIM; and qEW-10, mapped by CIM only), two major QTLs for GWP on chromosome 2 and 7 respectively (qGWP-6, mapped by CIM only; qGWP-10, mapped by CIM and GCIM), two major QTLs for GWP on chromosome 5 (qGWP-8 and qGWP-13, detected by CIM and GCIM respectively), and three major QTLs for GWP on chromosome 9 (qGWP-11 and qGWP-12, mapped by CIM; and qGWP-14, mapped by GCIM).
comparison of QtL regions. In this study, a total of thirty-three QTLs and fifty-four QTLs were respectively mapped by using CIM and GCIM methods for the seven yield-related traits. However, more than half of these QTLs were not stable, implying that these traits were controlled by multiple minor genes. The closely www.nature.com/scientificreports www.nature.com/scientificreports/ connected QTLs might be one locus, such asqEAD1 and qEAD1-1 on chromosome 1, qCD-3 and qCD3-1 on chromosome 1, qKRN-2 and qKRN-2-1 on chromosome 2. In addition, some of the detected QTLs in this study tended to have pleiotropism phenomenon. For example, qEAD-1 and qCD-3 were located on the same position on chromosome 1, but they were related to EAD and CD respectively. qEW-8 and qGWP-9 were located on the same position on chromosome 5, but they were related to EW and GWP respectively. qEW-10 and qGWP-11 located on the same position on chromosome 9 were related to EW and GWP respectively. And qKW-7, qEW-9 and qGWP-10 were located on the same position on chromosome 7, but they were related to KW, EW and GWP respectively.
Identification of candidate DEG for grain weight. RNA-seq procedure was carried out on an Illumina NovaSeq instrument. Totally, 18 RNA samples from SG5 and SG7 were collected form grains at three developmental stages (three biological replicates per stage) and used for RNA-seq analysis. The Pearson correlation analysis for RNA-seq data between samples was showed in Supplemental Table 2 and Fig. 2. A total of 932,738,890 with average 51, 818,827 clean reads on each chromosome were obtained after we filtering and performing quality control against the raw reads. (Supplemental Table 1). The clean reads were mapped onto the maize B73 genome by TopHat v. 2.0.12 39 software. For each replicate, 40,277,620-62,763,408 reads were obtained. Of these reads, 81.97-86.23% were mapped to the reference genome as unique and multiple matches (Supplemental Table 1). Based on these mapped reads, DEGs, novel transcripts, alternative splicing and events were detected. The fragments per kilobase of exon per million fragments mapped (FPKM) model was used for calculating all genes expression values. The threshold of corrected P-value 0.05 and log 2 (Fold change) of 0.5 (absolute value) were set as the thresholds for pair-wise comparison to detect DEGs between SG5 and SG7. To analyze the functions of DEGs, Gene Ontology function enrichment was conducted by Blast swiss prot database (website). The top significant biological process (BP), cellular component (CC), and molecular function (MF) terms were shown in Supplemental Fig. 3. In the BP group, DEGs are enriched in carbohydrate metabolic and lipid metabolic processes. In the MF group, DEGs are enriched in catalytic activity, hydrolase activity, oxidoreductase activity and so on. In the CC group, DEGs enrichment do not reach up to a significant level.
To decrease the candidate DEG number, we focused on those overlapped within the physical intervals of those confirmed QTLs related to grain weight. These QTLs were located on chromosomes 1, 2, 5, 7, and 9, including qKW-7, qEW-9, qEW-10, qGWP-6, qGWP-8, qGWP-10, qGWP-11, qGWP-12, qGWP-13 and qGWP-14. There were more than 1000 protein-coding genes in the mapped physical intervals in total, but after removing those genes that were not DEGs between parents, only 213 genes left and applied for further comparative genomics analysis and candidate gene prediction (Supplemental Table 3). candidate gene prediction. In ubiquitin-proteasome degradation pathway, grain size/weight was regulated by the interactions among three kinds of protein enzymes, i.e., ubiquitin protein ligase E3, ubiquitin-binding enzyme E2 and ubiquitin-activating enzyme E1 30 .The functional mechanism of DA2 40 in Arabidopsis and GW2 18 in rice in regulating grain size were both involved in the ubiquitin-proteasome degradation pathways. Previous www.nature.com/scientificreports www.nature.com/scientificreports/ studies have proved that abscisic acid (ABA) plays functions during seed development, and ABA-deficient Arabidopsis mutants produce seeds with increased size and mass 41 . MADS-box transcription factor genes have also been proven to play critical role in regulating grain size/weight indirectly 42 . FEM111/AGL80 is such a kind of gene in Arabidopsis 43 . Among the candidate genes in the intervals of qGWP-6, there were two genes GRMZM2G097089 and GRMZM2G158191 encoding E3 ubiquitin protein ligase. At qGWP-8 interval, one gene GRMZM2G169994 encoding E3 ubiquitin protein ligase was located. At the physical interval of the three identical QTLs qKW-7, qEW-9 and qGWP-10, there were two genes, one was GRMZM2G113039 that encodes E3 ubiquitin protein ligase and the other was GRMZM2G134480 that encodes ubiquitin-activating enzyme E1. At the physical interval of the two consistent QTLs qEW-10 and qGWP-11, two genes GRMZM2G057959 and GRMZM2G128953 were located and encoded ABA receptor MADS-box transcription factor respectively. And at qGWP-12 interval, one gene GRMZM2G036697 that encodes probable E3 ubiquitin-protein ligase was located.
These eight genes were chosen as candidate genes for grain weight in maize for the future prospects (Table 5).

Discussion
High yield is a permanent objective for maize breeders. Studies on QTL mapping and gene dissection for yield-related traits have become a research focus in maize in recent years 37,[44][45][46] . It is of great significance for maize breeder to map QTLs, to predict candidate genes, to clone functionlal genes and to clarify gene's genetic mechanism for yield-related traits, especially for grain weight. Compared to rice and Arabidopsis thaliana, the studies on maize grain weight candidate genes were relatively slow. Stable QTLs are useful for marker-assisted breeding. QTLs mapped in one population might well be not detected in another population. Thus, it is critical important for QTLs to be confirmed to rule out statistical anomalies while used in marker assisted breeding 7 . To validate QTLs, the first scheme is to confirm them in other mapping populations, the second scheme is to confirm them in different generations from the same crossing, and the third is to confirm QTLs using the same population evaluated in multiple locations and in multiple years. Mapping QTLs in early generations is that the QTLs with beneficial effect for breeding lines in early generations can be transferred to the high generations directly. In addition, the QTLs detected in early generation are more than those in late generation, especially for those minor QTLs 47 . In this study, an F 2:3 population developed from the cross between SG-5 and SG-7 was used for mapping QTLs for seven yield-related traits. In terms of phenotypic segregation issue in F 2:3 , Zhang et al. (2004) proposed to combine the F 2 phenotypes with the F 2:3 average phenotypes to further increase the power of QTL mapping 48 .
In this study, the F 2:3 phenotypic value were the mean value of eight F 2:3 individuals of F 3 progeny derived from F 2 plant selfing. The confirmed QTLs for yield-related traits are considered as stable QTLs and could be applied in marker assisted breeding, gene cloning and function analyzing, etc. Compared to low-density linkage map, Hori et al. 49 indicated that the higher-density linkage maps were more beneficial for QTL mapping, that is, the markers tightly linked to target QTL are more effective in marker assisted breeding 49 . Furthermore, it is possible to separate two closely linked QTLs by a high-density map with high resolution 49 . In this study, a high -density linkage map developed from an earlier study result 34 was applied in QTL mapping. However, the confidence intervals of identified QTLs always were not narrow enough to predict candidate genes directly. Gene density is within a wide range of 0.5-10.7 genes per 100 kb in maize genome 50 . It is of great significance to study how to decrease the number of candidate genes located on confirmed QTL intervals. For the purpose, a common approach is to develop advanced generation population such as NIL(near isogenic line) population, thus to narrow down the QTL confidence interval significantly, even to clone the QTL or gene based on primary mapping results 51,52 . However, the construction of NILs is time -consuming and tedious, which limit the usage of large amount of objective QTLs in marker assisted breeding to a certain extent. In this study, the DEGs between parents were obtained from transcriptomic profiling of grains at different developmental stages and overlaid onto the confirmed QTL intervals to predict candidate genes for grain weight. Based on these DEGs, the candidate genes in the physical intervals of the ten confirmed QTLs (qKW-7, qEW-9, qEW-10,qGWP  were decreased from over 1000 to 213 (Supplemental Table 3). Comparative genomics analysis was carried out to further predict candidate genes. Finally, a total of eight genes that might be involved in ubiquitin-proteasome degradation 30 , phytohormone signaling and transcription factor pathways were chosen as candidate genes controlling grain size and grain weight in maize.
In this study, a total of ten QTLs associated with grain weight were confirmed as stable QTLs, and eight candidate genes were predicted. All these results would be basis for cloning related functional genes and marker-assisted breeding.

Methods
Development of the F 2:3 mapping population and field trails' investigation. The F 1 hybrid seeds were obtained from an intraspecific cross between two inbred lines SG5 and SG7 in Liupanshui, Guizhou in 2013 summer. There are significant differences in yield-related traits for the two inbred lines ( Table 1) www.nature.com/scientificreports www.nature.com/scientificreports/ High density linkage map, QtL analysis and validation. Methods of extracting genomic DNA, sequencing genotype, grouping sequence data, identifying single nucleotide polymorphisms (SNPs) and constructing high-density bin map were exhibited in our previous study 34 . QTL were detected by two methods: one is CIM method included in QTL Cartographer v. 2.5, and the other is genome-wide composite interval mapping (GCIM) method 53 . GCIM method was applied by using random model with maximum likelihood (ML) function. For CIM and GCIM mapping procedures, the walking speed was 1 cM and the logarithm of odds (LOD) threshold value was set up as 3.0. The position of a significant QTL was determined by its LOD peaks. The positive or negative additive effect of a QTL indicated that the increase or decrease in phenotypic value for a trait is provided by the alleles from SG5 or SG7. The graphic of QTLs positions on 10 maize chromosomes were drawn by MapChart 2.32 software 54 . The physical intervals between the QTLs identified in F 2:3 population were compared with the intervals of QTLs mapped using a F 2 population previously 34 , and Those QTLs both detected from F 2 and from F 2:3 populations show overlapped physical intervals will be regarded as confirmed QTLs.
RnA sample preparation. The two maize inbred lines SG5 and SG7 were grown at the Panxian Maize Breeding Station, in Hainan, China in November 2016. Ear shoot were covered before silking. Some parental plants were hand pollinated when the length of corn silk was about 5 cm. Kernels were sampled from SG5 and SG7 on the 5-, 10, and 15-the days after pollination. Each sample consisted of three biological replicates in parallel. All samples were collected and frozen immediately with liquid nitrogen and stored in refrigerator at −80 °C. The total RNA of kernels at different growing stages was extracted by using TRIzol reagent (Invitrogen). One percent of agarose gels were used for RNA degradation and contamination. The NanoPhotometer ® spectrophotometer (IMPLEN, CA, USA) was used for checking RNA purity. Qubit ® RNA Assay Kit in Qubit ® 2.0 Flurometer (Life Technologies, CA, USA) was used for measuring RNA concentration. The RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA) was used for assessing RNA integrity. The Illumina NovaSeq platform was then applied to RNA-seq.
illumina sequencing and data analysis. A total of 18 samples with three repeats were collected and sequenced at the Illumina NovaSeq platform. Raw reads with fastq format were firstly handled by in-house perl scripts. Clean reads were then obtained after deleting reads containing adapter and ploy-N and removing reads of low quality in raw data. The GC content, Q20 and Q30 of the clean reads were calculated. High quality clean data were then carried out for further downstream analyzing. Maize reference genome and correlated files of gene annotation were downloaded directly from website (https://www.maizegdb.org/). Bowtie v. 2.2.3 and TopHat v2.0.12 39 were used for building reference genome index and aligning paired-end clean reads to the reference genome, respectively. The reads mapped to every gene were counted by HTSeq v. 0.6.1.
For each gene, the expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs (FPKM) was calculated by analyzing the gene length and mapped reads. FPKM is a widely accepted methodfor evaluating gene expression levels based on sequencing depth effect and gene length of the read count simultaneously 55 . The DEGSeq R package (v. 1.20.0) was used for analyzing differential expression genes between two conditions. The P value was corrected by using the Benjamini & Hochberg method. The threshold of corrected P-value 0.05 and log 2 (Fold change) of 0.5 (absolute value) were considered as significantly differential expression. The GOseq R package was used for analyzing Gene Ontology (GO) enrichment of DEGs. The GO term with corrected P-value less than 0.05 was considered as DEG.
Screening of candidate DEGs associated with QTLs for grain weight. In this study, the obtained RNA-seq data were used to explore DEGs between parental lines SG5 and SG7. Pair-wise comparison of transcriptomes between SG5 and SG7 was conducted for detecting DEGs. The DEGs were overlaid onto the physical intervals of confirmed QTLs to predict candidate genes for grain weight in maize.  Table 5. Predicted candidate genes for grain weight. h Start physical position of the gene; m End physical position of the gene; n Gene lenth; p LogFC: Log 2 ratio, number of folds the gene is differentially expressed in RNAseq; q RCP1/RCP2: different of readcounts between P 1 and P 2 ; Positive sign of logFC indicates gene transcript expressed high in SG5 while negative sign indicates gene transcript expressed high in SG7.


extfieldsofstudy: Biology Medicine

provenance: AgriculturalAndFoodSciences,2019-2020-Agricultural And Food Sciences-2019.jsonl.gz.gz

text:
Validation of reference genes for gene expression studies in tartary buckwheat (Fagopyrum tataricum Gaertn.) using quantitative real-time PCR

Quantitative real-time reverse transcriptase polymerase chain reaction is a sensitive technique for quantifying gene expression levels. By implementing three distinct algorithms (geNorm, normFinder and BestKeeper), we have validated the stability of the expression of seven candidate reference genes in tartary buckwheat, including FtSAND, FtCACS, FtExpressed1, FtGAPDH, FtActin, FtEF-1a and FtH3. In this study, the results indicated that FtCACS and FtSAND were the best reference genes for ‘abiotic cotyledons’, FtExpressed1 and FtEF-1α were the best reference genes for aluminium treatment, FtCACS and FtExpressed1 performed the best for the immature seed stage, FtCACS was best for the abiotic treatment, and FtH3 appeared to be the most suitable reference gene for the abiotic treatment in hypocotyls and all samples in this study. In contrast, FtActin and FtGAPDH are unsuitable genes. Our findings offer additional stable reference genes for gene expression research on tartary buckwheat at the immature seed stage and under abiotic treatment.

INTRODUCTION
Quantitative real-time reverse transcription polymerase chain reaction (qRT-PCR) has become the most prevalent quantification method used in assays of gene expression on account of its specificity, accuracy, efficiency and high sensitivity (Jiang et al., 2014). The data from qRT-PCR can be analysed using absolute or relative quantification. Absolute quantification yields an exact gene copy number by transforming quantification cycles (Hruz et al., 2011) into a standard curve. Relative quantification relies on internal control genes as references to present the qRT-PCR data of the target genes (Ginzinger, 2002). Relative quantification has become a more widely used method in gene expression assays because most researchers mainly focus on the analysis of differences in gene expression . To normalise qRT-PCR data, a reference gene is needed, and it should be universally valid, with a similar expression level across all feasible cells, tissue samples and experimental treatments (Huang et al., 2014). However, ideal reference genes do not actually exist (Remans et al., 2014). The selection of the appropriate reference gene is a critical step in controlling the variability of samples when using a sensitive qRT-PCR technique (Tong et al., 2009). If the reference gene is appropriate, the discrepancies that may exist in terms of the initial sample amount, RNA integrity, RNA recovery and the efficiency of cDNA synthesis will be eliminated. Statistical algorithms, such as geNorm (Vandesompele et al., 2002), normFinder (Andersen, Jensen & Orntoft, 2004) and BestKeeper (Pfaffl et al., 2004), have been developed to help to select appropriate reference genes.
Most studies about gene expression in bacteria , yeast (Cankorur-Cetinkaya et al., 2012) and mammals (Terzi et al., 2010) now include reference gene validation (Chapman & Waldenstrom, 2015). However, the validation of reference genes in plants has received very little attention, and housekeeping genes tend to be used as references without any appropriate validation (Gutierrez et al., 2008). , Czechowski et al. (2005 were the first to present a list of stably expressed Arabidopsis genes under a large range of experimental conditions, and the evidence clearly showed that several genes are expressed more stably than traditional reference genes. Since then, suitable reference genes for gene expression studies have been reported in higher plants such as flax (Huis, Hawkins & Neutelings, 2010), soybean (Hu et al., 2009, tomato (Lovdal & Lillo, 2009), Pisum sativum (Die et al., 2010), Brachypodium distachyon (Hong et al., 2008 and carrot (Tian et al., 2015). However, reference gene stability is not consistent across experimental conditions and plant species. Consequently, it is necessary to find additional reference genes for different conditions and species.
Tartary buckwheat (Fagopyrum tataricum Gaertn.) belongs to the Polygonaceae family (Kim et al., 2013), and as an important functional food material, it has a relatively high flavonoid content. Tartary buckwheat is capable of thriving in regions with poor soil or harsh climates (Kim et al., 2009). The results of genome sequencing also show that tartary buckwheat has a remarkable ability to cope with highly variable environmental stress, including drought, salinity, UV-B and cold (Zhang et al., 2017). The stress resistance of tartary buckwheat is mainly due to its abundant flavonoids (Suzuki, Honda & Mukasa, 2005). In particular, tartary buckwheat is a naturally aluminium (Al) tolerant species . Abiotic stresses and flavonoid metabolism regulate the expression of these genes in plants at both transcriptional and post-transcriptional levels. There are many reports on tartary buckwheat under various experimental conditions in which gene expression is normalised to a reference gene (H3) for semi-quantitative RT-PCR or qRT-PCR (Bai et al., 2014). However, to date, there is no systematic strategy to analyse tartary buckwheat reference genes at the immature seed stage or under abiotic stress.
In this work, we aimed to evaluate the potential use of different reference genes for internal normalisation to more accurately measure the expression level of genes of interest in tartary buckwheat. Seven candidate reference genes were selected, and the stability of their expression was assessed in tartary buckwheat at the immature seed stage and under different abiotic stress treatments. Evaluating the stability of the expression of candidate reference genes depends on statistical analysis. Three different statistical software programs (geNorm, normFinder and BestKeeper) were used to calculate the variability of the expression of the candidate genes and determine which were the most suitable.

Plant materials and treatments
Tartary buckwheat ('Xiqiao No. 2') seeds were grown in the field on a farm (XTBG; 29 59′ N, 102 59′E; 800 m elevation) at Sichuan Agricultural University, Ya'an, Sichuan, China (Li et al., 2012). Tissues, including roots, stems, leaves, flowers, immature seed 1 (seed formation started) and immature seed 2 (seeds in the milk) were collected at the immature seed stage (Gupta et al., 2011). The 7-day-old seedlings were stressed with saline or drought by adding 100 mM NaCl or 20% PEG 2000, respectively, to the medium. The 7-day-old seedlings were exposed in a chamber at 4 C with a 16 h photoperiod for cold treatment (Gao et al., 2016). UV-B treatment was conducted under UV-B (302 nm, 0.1 mW/cm 2 ) in a chamber. After 0, 2, 4, 6, 12 and 24 h of treatment, all stressed seedlings were collected and separated into cotyledons and hypocotyls. For the Al treatment, the samples were processed according to a previous report (Zhu et al., 2015). Root tips (zero to two cm) and basal roots (two to four cm) were sampled under both -Al and +Al conditions. All samples were collected in two biological replicates, and RNA was extracted immediately.

Total RNA isolation and cDNA synthesis
Total RNA was isolated from various samples with an RNAout 2.0 kit (Tiandz, Beijing, China) according to the manufacturer's instructions. To remove trace DNA from samples, total RNA extractions were treated with RNase free DNase I. The RNA integrity was detected using 2% agarose gels. A Bio-RAD smart spec TM plus spectrophotometer was used to determine the RNA purity. cDNA was synthesised with a PrimeScript TM RT reagent kit and gDNA Eraser (Perfect Real Time; TaKaRa, Dalian, China).

Selection of candidate reference genes and design of qRT-PCR primers
Potential homologues of the seven candidate genes were identified from the transcriptome sequencing data of tartary buckwheat ('Xiqiao No. 2') (Yao et al., 2017). Primers were designed using Primer Premier 5.0. All primers used in this research are listed in Table 1. The specificity of the amplification was assessed based on the presence of a single band of the expected size in a 1.5% agarose gel following electrophoresis and a single peak in the qRT-PCR melting curve.

Quantitative real-time PCR
The qRT-PCR procedure was designed in accordance with MIQE guidelines (Bustin et al., 2009). qRT-PCR was executed in a CFX96 Real Time PCR system (Bio-Rad, Hercules, CA, USA) with a SYBR Premix EX Taq kit ( TaKaRa, Dalian, China) in a total reaction volume of 15 mL that included 10 mL of SYBR Green mix, primers at 0.5 mM each and one mL of cDNA. The amplification conditions were as follows: 95 C for 30 s and 40 cycles of 95 C for 5 s and 60 C for 20 s. A melting curve from 60 to 95 C was used to verify the specificity of the PCR amplification. All genes were amplified from cDNA. 10-fold serial dilutions of cDNA samples from young leaves were used to establish the standard curves to calculate the amplification efficiency of each primer pair.

Statistical analysis
Three software programs (geNorm v3.5, the Excel add in of normFinder v0.953 and BestKeeper v1) were used to analysis the stability of reference gene expression across all experimental sets (Andersen, Jensen & Orntoft, 2004;Pfaffl et al., 2004;Vandesompele et al., 2002). SPSS v17.0 was used to calculate the span of Cq values for each gene by drawing a box-whisker plot, and the expression levels of FtSTAR and FtDFR were showed by mean ± standard deviation (SD).

RNA solution and quality
A series of 58 samples from tartary buckwheat were divided into six different groups. 'Abiotic cotyledons' was composed of cotyledons from all stress-treated samples. 'Abiotic hypocotyls' included hypocotyls from all stress-treated samples. 'Abiotic total' was composed of 'abiotic cotyledons' and 'abiotic hypocotyls'. 'Al treatment' comprised all samples that were treated with Al. 'Immature seed stage' consisted of six tissues (roots, stems, leaves, flowers, immature seed 1 and immature seed 2) at the immature seed stage based on the state of the seed. Finally, 'total' included all the samples in this study. Protein and organic pollutants were isolated and removed from all samples via RNA extraction. The total RNAs, with A 260 /A 280 ratios of 1.8-2.0, were reverse transcribed into cDNA as templates for qRT-PCR detection.

Expression profiles of candidate reference genes
We selected the best reference genes (among the seven candidate genes) from the seven candidate genes in six different groups (sequencing data was shown in Fig. S1). The amplification efficiencies (E) of all reactions ranged from 95.5% to 107.5% and were calculated from standard curves with good linear relationships (R 2 > 0.99). Single-peak melting curves were obtained for all qRT-PCR amplifications (Fig. S2) (Fig. 1).

Software: geNorm, NormFinder and BestKeeper
After a simple comparison of the raw Cq values, the three software programs were used to further analyse the stability of the expression of the seven candidate genes. The genes were ranked in descending order in terms of the stability of their expression in each of six groups (Table 2), and the details of the cotyledons and hypocotyls under abiotic stress are provided in Table S1. GeNorm performed a stepwise exclusion of the most unstable gene and then recalculated M until only two genes remained, and these two genes had the most stable expression (Vandesompele et al., 2002). The gene with the lower M-value was considered to have the most stable expression. An M-value limit of <1.5 was  suggested by geNorm. All the results from geNorm were lower than 1.5 in 'abiotic cotyledons' and 'total'; FtCACS and FtH3 (M = 0.33, 0.60) were chosen as the most stable genes, while FtH3 and FtSAND (M = 0.65, 0.55) were identified as the most stable genes in the 'abiotic hypocotyls' and 'abiotic total' groups, respectively, In 'immature seed stage', FtCACS and FtExpressed1 (M = 0.16) were recognised as the most stable genes. Finally, FtEF-1a and FtExpressed1 (M = 0.16) were the most stable genes in 'Al treatment'.
To determine the optimal number of reference genes required for accurate normalisation, geNorm was used to calculate the pairwise variation (Vn/Vn+1) between the sequential normalisation factors (NFs) (NFn and NFn+1). As suggested, a threshold value of 0.15 was adopted. As depicted in Fig. 2, pairwise variation analysis indicated that the ideal number of reference genes may be different for the different groups. For instance, only two genes are necessary for normalisation for 'immature seed stage', 'abiotic cotyledons' and 'Al treatment', but the pairwise variation of the other three groups was above the threshold of 0.15. The pairwise variations for cotyledons and hypocotyls under different abiotic treatments are provided in Fig. S3.
Unlike geNorm, normFinder (Andersen, Jensen & Orntoft, 2004) depends on a variance estimation approach, which allows the comparison of inter/intra-group variation. Genes with the lowest average expression stability values are the most stable. In the normFinder analysis, FtH3 and FtCACS were the most stable genes in 'abiotic cotyledons', and FtSAND and FtExpressed1 were the most stable genes under 'Al treatment'. In addition, the most stable reference genes in 'immature seed stage' were FtSAND and FtGAPDH, while the other three groups had the same two most stable reference genes (FtSAND and FtCACS). The least stable gene in all groups was FtActin.

Reference gene validation
To validate the availability of a reference gene, the expression levels of FtSTAR under Al treatment and of the FtDFR gene under UV treatment were determined using the seven candidate reference genes for normalisation. For Al treatment, STAR has a conserved response in plants, and this has been shown in previous reports, such as reports on rice STAR2 (Huang et al., 2009), Arabidopsis ALS (Larsen et al., 2005) and tartary buckwheat FtSTAR2 (Zhu et al., 2015), whose expression increased after exposure to Al. The expression of FtSTAR2 was greatly increased (Fig. 3), which was reinforced by the Al-induced expression of STAR2. The use of the most favourable reference gene (FtExpressed1) resulted in the greatest variation, resulting in increases of approximately 1.67-fold in root tips and 6.24-fold in basal roots, and the use of the other most stable gene (FtEF-1a) resulted in increases of approximately 1.42-fold and 6.19-fold, respectively. Finally, the use of the least stable gene (FtActin) resulted in increases of approximately 1.90-fold in root tips and 2.48-fold in basal roots. These genes also showed a pattern of increased expression, regardless of whether stable or unstable reference genes were used for normalisation. However, the relative expression using the most stable genes showed the higher fold increases compared with the unstable genes. Therefore, FtExpressed1 and FtEF-1a are the most suitable genes under Al treatment. Dihydroflavonol 4-reductase (DFR) catalyses the reduction of dihydroflavonols to leuco-anthocyanins and is a key enzyme in the biosynthesis of anthocyanins (Yuan et al., 2007). The contents of anthocyanin increased under UV treatment in tartary buckwheat hypocotyls (Eguchi & Sato, 2009). In hypocotyls under UV treatment, the relative quantification of FtDFR using the most stable genes (FtCACS, FtH3 and FtActin) for normalisation exhibited similar expression patterns but different expression levels. When other less stable genes were used for normalisation, the results obviously differed from those with the most stable genes (Fig. 4). The expression levels of the FtDFR gene under UV treatment showed similar trends after normalisation to stable genes in hypocotyls. Thus, we selected the most stable reference genes for gene expression normalisation in tartary buckwheat.

DISCUSSION
SYBR Green I is a fluorescent reporter dye used in qRT-PCR; when it binds double-stranded DNA, its fluorescence increases nearly 1,000-fold (Morrison, Weis & Wittwer, 1998). The problems in qRT-PCR caused by the variability of RNA templates and inappropriate data normalisation are obvious and widely known but disregarded (Nolan, Hands & Bustin, 2006). To address these questions, proper reference genes should be used. Thus, validating candidate reference genes and selecting stable reference genes are becoming general interests for researchers. Candidate reference genes can be selected from reference gene validation papers that use the same experimental conditions or closely related species (orthologies) or by mining transcriptomic data for stably expressed genes or traditional housekeeping genes (Hruz et al., 2011). The seven candidate reference genes tested in tartary buckwheat include three genes that are stably expressed in common buckwheat (SAND, CACS and Expressed1), three commonly used reference genes (GAPDH, Actin and EF-1a) and one unique gene (H3). Neither high (>30) nor low (<15) Cq values are recommended based on general guidelines (Wan et al., 2010). In this study, the Cq values of the reference genes ranged from 14.98 (FtGAPDH) to 29.65 (FtActin), suggesting variable expression. Although the raw Cq value comparison can provide a rough estimate of the stability of gene expression, it is not sufficient to accurately evaluate the expression patterns of reference genes.
Scientists have developed several approaches to select stable reference genes for normalisation. However, to date, there is no consistent algorithm that can be used to test gene stability. These software programs have their own merits and drawbacks: NormFinder can avoid the misinterpretations caused by the artificial selection of co-regulated genes (Andersen, Jensen & Orntoft, 2004); geNorm is based on the assumption that none of the analysed genes are co-regulated (Vandesompele et al., 2002); and geNorm can estimate the fewest number of reference genes needed for accurate normalisation (Stamova et al., 2009). It seems possible that more reliable controls can be obtained using a combination of different algorithms (De Almeida et al., 2010). We used three software programs that were based on different statistical approaches to assess seven candidate genes in tartary buckwheat. The distinct statistical algorithms are likely to generate inconsistent rankings of stability. After integrating the outcomes of the three programs above in this study, we recommend the reference genes in tartary buckwheat as follows: FtCACS + FtSAND for 'abiotic cotyledons'; FtExpressed1 + FtEF-1a for 'Al treatment'; FtCACS + FtExpressed1 for 'immature seed stage'; FtCACS for 'abiotic total'; and FtH3 for 'abiotic hypocotyls' and 'total'. The H3 gene is the most commonly used internal gene for normalisation in tartary buckwheat (Luo et al., 2016). Our results indicate that FtH3 is a stable gene under abiotic treatment. However, it is not stable across different organs or after Al treatment. In common buckwheat (Demidenko, Logacheva & Penin, 2011), CACS and Expressed1 were validated as the most stable genes in the development and fruit stages. FtCACS and FtExpressed1 being the most stable reference genes in the immature seed stage of tartary buckwheat supports the statement that the orthologues of identified reference genes could serve the same purpose in other species. In tomato, similar results were obtained in a study on reference gene selection (Exposito-Rodriguez et al., 2008).
GAPDH and Actin are the most generally used reference genes for the analysis of gene expression in various plant species (Kumar et al., 2011). However, our analysis indicates that GAPDH and Actin are not reliable genes for comparative expression analysis. The cause of these fluctuations at the gene expression level is probably that GAPDH and Actin have several biological functions, such as participating in the glycolytic pathway and other processes (Stürzenbaum & Kille, 2001). Actin supports the cell and determines its shape, and it also takes part in other cellular functions (Kravets, Yemets & Blume, 2017). Although Actin is a stable reference gene in developmental studies, it is not stable under various conditions.

CONCLUSIONS
As far as we know, our studies take advantage of three software (geNorm, NormFinder and BestKeeper) to analysis the stability of seven candidate reference genes for the first time in tartary buckwheat using qRT-PCR. The three software identified slightly different genes as most suited for normalisation, prompting us to merge the data. The result showed that the expression of GAPDH or Actin is unstable across all samples. We also provide a list with the stable reference genes in six group and under certain conditions.

ADDITIONAL INFORMATION AND DECLARATIONS Funding
This work was supported by the National Natural Science Foundation of China (31500289). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Grant Disclosure
The following grant information was disclosed by the authors: National Natural Science Foundation of China: 31500289.


text:
The Effects of Different Lead Pollution Levels on Soil Microbial Quantities and Metabolic Function with/without Salix integra Thunb. Planting

: Background and Objectives: Salix integra Thunb., a fast-growing woody species, has been used in phytoremediation in recent years. It has the potential to accumulate high amounts of lead (Pb) in its growth, however, its effects on soil microbial community structure and function during its phytoextraction processes are not well understood, especially at different pollution levels. Materials and Methods: In our study, we set unplanted and planted Salix integra in areas with four levels of Pb treatments (0, 500, 1000, and 1500 mg/kg). After six months of planting, the rhizospheric soil, bulk soil, and unplanted soil were collected. Soil properties and microbes participating in nitrogen and phosphorus cycling were measured, following standard methods. Microbial metabolic functions were assessed using a Biolog-ECO microplate. Results: The bacteria (nitrogen-ﬁxing bacteria, ammonifying bacteria, inorganic phosphorus-solubilizing bacteria, and nitrosobacteria) all increased in the 500 mg/kg treatment and decreased in the 1500 mg/kg treatment compared with the 0 mg/kg treatment, especially in rhizospheric soil. The microbial metabolisms decreased along with the increase of Pb levels, with the exception of the rhizospheric soil with a 500 mg/kg treatment. The metabolic patterns were relative to the pollution levels. The utilization of carbohydrates was decreased, and of amino acids or fatty acids was increased, in the 500 mg/kg treatment, while the opposite occurred in the 1500 mg/kg treatment. The values of soil properties, microbial quantities, and metabolic activities were higher in rhizospheric than bulk soil, while the differences between bulk and unplanted soil were different among the different Pb treatments. The soil properties had little effect on the microbial quantities and metabolic activities. Conclusions: S. integra planting and Pb levels had an interactive effect on the microbial community. In general, S. integra planting promoted microbial quantities and metabolic activity in rhizospheric soil. Lower Pb pollution increased microbial quantities and promoted the utilization of amino acids or fatty acids, while higher Pb concentrations decreased microbial quantities and metabolic activities, and promoted the utilization of carbohydrates.

Introduction
Lead (Pb) is considered seriously harmful for human health because of its toxicity, persistence, bioaccumulation and biomagnification throughout the food chain [1], and high degree of accumulation in soil [2]. In the past five decades, more than 800,000 tons of Pb have been released into the environment globally, most of which has accumulated in soil and thus caused serious heavy metal pollution [3]. In addition, the problem of the soil environment in industrial and abandoned mining land is significant and serious. Therefore, it is a matter of great importance to remove Pb from soil in order to maintain a safe food chain and healthy environment. In general, physicochemical remediation technologies include soil excavation, soil cleaning, electrokinetic remediation, and extraction. However, these methods all have the problem of secondary pollution.
In the past three decades, phytoremediation technology has been proposed as a cost-effective and ecofriendly method for cleaning up polluted sites without secondary pollution [4]. Among the different techniques of phytoremediation, phytoextraction and phytostabilization are particularly relevant to alleviating metal pollution in soil [5]. Hyperaccumulators have also been strongly suggested as advantageous for phytoremediation. Hyperaccumulating plants exhibit foliar metal concentrations that are a minimum of 2-3 orders of magnitude higher than foliar concentrations found in other plants growing in normal soils, and at least one order of magnitude greater than the range in foliar concentrations of plants growing on metalliferous soils (established thresholds for hyperaccumulators growing in their natural habitats include 100 µg/g for Cd, Se, and Tl; 300 µg/g for Co, Cu, and Cr; 1000 µg/g for Ni, Pb, and As; 3000 µg/g for Zn; and 10,000 µg/g for Mn) [6]. However, the application of hyperaccumulators for phytoremediation is often limited by their slow growth rate and low biomass yield [7][8][9]. Trees, especially fast-growing woody species, such as Populus [10] and Robinia pseudoacacia L. [11], with their large biomass and deeper, more integrated root systems, have provided a unique means for the deep phytoremediation of soil [12,13].
Phytoremediation technology is based on the interactions of plants and microbes [14]. Root exudates stimulate the growth of specific bacterial and fungal populations, and change the microbial metabolic activity in rhizospheric soil [15][16][17]. In return, microbes modify the chemical composition of root exudates, the bioavailability of heavy metals, and soil physicochemical properties by their metabolic activities [18][19][20]. For example, plant colonization in mine tailings dramatically remodels the structure and function of the soil microbial community [19,21,22]. Zn pollution leads to decreased levels of basal respiration and ammonia-oxidizing bacteria, while Thlaspi caerulescens growth increases the values of substrate-induced respiration and total bacteria [23]. Similarly, planting of R. pseudoacacia also has a significant influence on bacterial communities in Pb/Zn contaminated soils [11]. Poplar phytomanagement regimes have led to decreased bioavailability of Zn and Cd, in concert with changes in the microbial communities [24]. Meanwhile, many fungi could alleviate heavy metal toxicity to plants and influence heavy metal accumulation and transportation [25][26][27][28][29]. Although advanced genetic editing of plants and bacteria may soon enhance phytoremediation even further [30], the current lack of reports about the impacts of metal contamination on microbial metabolic functions during phytoremediation makes their study of utmost importance.
In recent years, willow species, which can be grown intensively for use in energy production, have been suggested for use in the remediation of metal contaminated soils worldwide [9,31]. Much work has been conducted to evaluate the ability of Salix integra Thunb. to absorb heavy metals [32][33][34][35][36], especially Pb, and accumulation of up to 456 mg kg −1 has been achieved in the shoots [37]. However, the effects of S. integra planting on the soil microbial community are not well understood, especially at different pollution levels, which limits our understanding of the phytoremediation processes of S. integra. Therefore, to better understand the phytoremediation mechanism of S. integra, we set unplanted and planted S. integra in areas with four levels of Pb treatments (0, 500, 1000, and 1500 mg/kg) in our study. We hypothesized that (1) S. integra planting would promote microbial metabolic activities and quantities of functional microbes in rhizosphere soil; (2) the effects of different Pb pollution levels on microbial metabolic activities and quantities of functional microbes would be different; and (3) microbes and microbial metabolic functions would be affected by the soil properties.

Study Site and Plant Material
The study site was located in the experimental station of Hebei Agricultural University, with a temperate climate (Baoding city, Hebei province, China. 38 • 45 21" N, 115 • 24 37" E). The mean annual temperature is approximately 13.0 • C and the annual precipitation is about 532 mm. The soil is a typical meadow cinnamon soil (a transition type from cinnamon soil to meadow soil). Salix integra is a shrub of the family Salicaceae. One year cutting seedlings of S. integra with the same growth vigor were selected for use in the study.

Experimental Set Up
The experiment was carried out by excavating four square plots below the ground level. The size of each plot was 4 m (length) × 4 m (width) × 0.6 m (depth). Perfluorinated ethylene-propylene copolymer (PEP) plastic sheets were used to line the bottom of the four plots, which played an anti-seepage role. Finally, the soils were backfilled correspondingly. This method avoids the disadvantages of pot experiments, in which soil compaction can occur easily and thereby limit plant growth. A rain shed with a transparent plastic sheet was built above the study area to control the experiment.
In our study, four treatment regimens of Pb were conducted in each of the plots. They were 0 mg/kg (CK), light pollution 500 mg/kg (LT), medium pollution 1000 mg/kg (MT), and high pollution 1500 mg/kg (HT), respectively. Pb was injected into the soil as Pb(NO 3 ) 2 solution in June 2017. The solution volume was calculated by the water holding capacity. The soils were subjected to aging and equilibration under natural conditions for three months. During the aging period, the soil moisture was maintained at approximately 75% of the maximum water holding capacity. Meanwhile, the soils were ploughed frequently. In each treatment, we set the unplanted area and planted S. integra area at the same time. The different treatment areas were separated by a PVC sheet. The layout of the study site is shown in Figure 1. The initial soil properties were as follows: organic carbon (2780 mg/kg), total nitrogen (240 mg/kg), available nitrogen (22.28 mg/kg), total phosphorus (510 mg/kg), available phosphorus (10.54 mg/kg), total potassium (8080 mg/kg), available potassium (278.4 mg/kg), cation exchange capacity (8.46 cmol/kg), and pH (8.14).
One year cutting seedlings were planted in May 2017. Eight trees were planted in each treatment. The control and treatment groups were all under the same management conditions (watering, loosening, weeding, etc.) during the study period.

Soil Collection
After six months growth, the rhizospheric soil, bulk soil, and unplanted soil were collected in November, in three repetitions under each treatment (sampled in plots 1, 2, and 4). The following method was utilized for rhizospheric soil: the roots were recovered from the soil cores (the diameter 10cm, and the depth 20-30cm), then gently shaken to remove any loosely adhering soil, and after that placed in a sterile plastic bag and shaken dramatically. The soils from five plants were pooled together, thus giving a mixed sample of rhizospheric soil as one replication. The bulk soils and unplanted soils were all obtained by soil-drilling. The bulk soils were collected around the tree, and the unplanted soils were randomly collected from the unplanted area. Soils from five random spots were mixed as one replication. Each replication sample was divided into five parts for the measurement of soil nutrients, pH, cation exchange capacity, quantities of soil culturable microbes relative to nitrogen and phosphorus cycling, and community level physiological profiling (CLPP). After removal of crop residues, the soil samples used for measuring culturable microbes and CLPP were transported to the laboratory immediately in sterilized plastic and stored at 4 • C. The tools were cleaned and disinfected between different treatments when sampling the soil.

Soil Properties Analysis
The soils were air dried and sieved to a diameter of <2 mm. The organic carbon (SOC) was measured by the potassium dichromate heating method after the soils were digested with H 2 SO 4 . Available nitrogen (AN) was analyzed by the alkali diffusion method. Available phosphorus (AP) was analyzed by Mo-Sb colorimetry (using a Vis spectrophotometer, 722S, Jinhua China) after the soil was extracted with sodium bicarbonate. Available potassium (AK) was analyzed with atomic absorption spectrophotometry (AA-680, SHIMADZU) after the soil was extracted with CH 3 COOH. The pH value was measured in deionized water (1:2.5 m/V) using a pH meter. The cation exchange capacity (CEC) was measured by the method of barium chloride-sulfuric acid exchange [38].

Community Level Physiological Profiling
CLPP is a useful tool for assessing microbial community metabolic function and functional diversity in both terrestrial and aquatic environments [41][42][43]. CLPP, also known as carbon source utilization patterns of soil microbial communities, is assessed using a Biolog-ECO micro-plate (Category NO. 1506 U.S.A.). One plate, containing 31 carbon sources and one blank (in triplicate), was used in each sample. The soil was sifted to remove impurities. A subsample of 10 g (dry weight equivalent) of each soil sample was suspended in 90 mL of 0.1% sterile sodium chloride (NaCl) solution, and then shaken for 30 min (200 r/min) at room temperature. The supernatant was diluted 10 3 times, and the diluted soil suspension (150 µL) was added to each well of a Biolog-ECO micro-plate. The micro-plate was incubated at 25 • C in darkness. Optical density (OD) was measured with SpectraMax384 (Molecular Devices, U.S.A.) at 590 nm (color and turbidity) and 750 nm (turbidity) wavelengths every 24 h until no more growth in the OD value could be observed. The procedure lasted up to eight days to ensure the carbon utilization saturation phase was reached in all samples [44].
For all samples, absorbance readings at 144 h were identified as the time point for further CLPP data analysis. Data were analyzed to determine metabolic diversity indices, according to Sofo's study [45]. The indices examined were as follows: average well color development (AWCD) that were utilized to reflect microbial metabolic activity [46], richness index (S), Shannon's diversity index (H'), Shannon's evenness index (E), Simpson's diversity index (D), McIntosh's diversity index (U), and McIntosh's evenness index (Z).

Statistical Analyses
Each data point represented the results of triplicate sample analysis, and each value was expressed as the mean ± standard error. A two-way ANOVA was applied to compare the impact of the two factors (Pb concentration and different sources of soil) on significant differences of indices among different treatments. The normality of the analyzed features was checked with the Shapiro-Wilk test. Depending on the distribution of the estimated parameters, either the least significant difference (LSD) in the parameter test or the Kruskal-Wallis test in the non-parameter test was used to check for significant differences among the treatments. In all analyses, differences were considered significant at a p value < 0.05. The above statistics were performed with the package SPSS Statistics 19.0 (IBM, Armonk, NY, America).
Principle component analysis (PCA) was performed using the data matrix (11 samples and 31 variables) of AWCD normalized carbon source utilization patterns data to further assess the differences among different treatments. Redundancy analysis (RDA) was conducted to test the relationship of soil properties and microbial communities. The analyses of PCA and RDA were all performed by the software CANOCO (Windows version 4.5) (Biometris-Plant research international, Wage-ningen, the Netherlands).

Soil Physicochemical Properties
The interactive effect (PbC × DS) of Pb concentration effect (PbC) and different sources of soil (DS) was significant on AN, AP, AK, and CEC, but not on pH ( Figure 2). The values of CEC were increased in LT and MT and decreased in HT, especially in rhizosphere (Figure 2d). The contents of AP and AK were strongly inhibited in Pb treatments no matter whether in rhizosphere, bulk, or unplanted soil. The inhibitory effects were higher with the increase of Pb concentration (Figure 2b,c), while the AN was promoted when compared with CK ( Figure 2a) and was highest in MT.
The effects of S. integra planting varied based on different Pb treatments. In general, the values of AK, AN, AP, and CEC were higher in rhizosphere than bulk soil, regardless of CK, LT, MT, or HT. However, the trends were complex when compared with unplanted soil. For example, the AK of unplanted soil was lowest in MT and HT, while it was highest in CK and LT. The effects of both Pb pollution and S. integra planting on pH were not significant.

The Quantities of Soil Culturable Microbes
The interactive effect of PbC × DS was significant on the four kinds of bacteria (Figure 3a-d). Compared with the CK, the quantities of IPB were increased in LT and decreased in MT and HT, whether in rhizosphere, bulk, or unplanted soil (Figure 3b). The AMB and NFB showed a similar trend only in rhizosphere soil, and they increased in LT and decreased in MT and HT when compared with the CK. However, the AMB decreased in bulk soil of Pb pollution treatments and increased in unplanted soil of HT, and the NFB decreased only in bulk soil of HT and increased in unplanted soil of LT compared with the CK (Figure 3a and d). The quantities of NB were increased in the MT and LT of rhizosphere soil, and decreased only in the unplanted soil of HT (Figure 3c). The four kinds of bacteria were all highest in rhizosphere, while the difference between bulk and unplanted soils varied among the different treatments.

Metabolic Activity of Microbial Communities
The dynamic patterns of AWCD values are shown in Figure 4a-e. The AWCD values were relatively small with 24 h of incubation, which indicates that the carbon sources were not utilized during that period. After that, the AWCD value increased gradually as time went on. The AWCD values changed with time in a similar trend among the different treatments. S. integra planting promoted the metabolic activities of microbial communities, which were all higher in rhizospheric soil, and the differences increased with the increase of Pb concentrations from LT to HT. The AWCD values of unplanted soil were higher than those of bulk soil in the LT. However, the differences between bulk and unplanted soil were small in MT and HT.
The carbon source utilization rates were basically stable at 144 h, so the AWCD values for 144 h were used for further analysis (Figure 4e). There was no significant difference between LT and CK in rhizospheric soil, and LT and MT in bulk soil. The others were reduced with increased Pb concentrations.

Utilization of Specific Carbon Sources
Thirty one carbon sources in Biolog-ECO plates were divided into four categories, according to the microbial metabolic pathway of three major nutrients. The categories were carbohydrates and their derivatives (CG, 12 kinds), amino acids and their derivatives (AG, 6 kinds), fatty acids and lipids (FG, 5 kinds), and metabolic mediates and secondary metabolites (MG, 8 kinds) [47].
Pb concentrations, carbon sources, and their interactive effect significantly affected the AWCD (Figure 5a-c). The utilization of specific carbon source types changed with different Pb levels, as compared with the CK. The CG, AG, and MG showed no significant differences between CK and LT, but four categories were all decreased in HT in rhizospheric soil (Figure 5a). The four categories were all decreased as compared with the CK in bulk soil, and the adverse effects were higher with increased Pb concentrations (Figure 5b). The highest utilization of specific carbon sources were FG in CK and CG in HT, while the FG, CG, and AG showed no significant differences in LT and MT in rhizospheric soil. FG was highest in LT and MT, and FG, CG and AG showed no significant differences in CK and HT in bulk soils. The FG and CG were highest in LT, CG was highest in HT, and FG, CG and AG showed no significant differences in MT in unplanted soil (Figure 5c). However, MG was lowest in all the treatments. Furthermore, we compared the relative utilization percentage of the specific carbon source in each treatment (Figure 5d). We found that the CG utilization was decreased both in rhizospheric soil and bulk soil, while AG was increased in rhizospheric soil, and FG was increased in bulk soil when compared LT with CK. In contrast, the CG utilization was increased both in rhizospheric soil and bulk soil, when MT and HT were compared with CK. Decreases of AG and MG were found, respectively, in rhizospheric and bulk soil under both MT and HT. The CG utilization was increased with the increase of the Pb concentrations in unplanted soil.

Microbial Community Metabolic Profiles
The carbon source utilization patterns (CSUPs) of microbial communities were visualized using principal component analysis (PCA). According to the results (Figure 6), the variance contribution rate of the first principal component (PC1) was 42.5%, and the second principal component (PC2) was 19.0%. All the samples were distributed along the positive direction on the PC1 axis.
The microbial CSUPs formed four groups. From right to left, the first group contained the rhizosphere of CK, LT, and MT and the bulk soil of CK, while the rhizosphere of HT was close to the second group, which included the bulk and unplanted soils of LT. The third and fourth groups contained the bulk and unplanted soils of MT and HT, respectively. This result indicates that the soil microbial community had a unique model of carbon utilization in rhizosphere soil, while the bulk and unplanted soils were distinctive in relation to the contamination levels.

Diversity Indices of Soil Microbial Communities
The interactive effect of PbC × DS was significant on S, U, and Z (Figure 7a,c,d). LT increased the values of S, and MT and HT decreased the values of U in rhizospheric soil as compared with CK. The values of S and U were all decreased in bulk soil, when compared with CK. The adverse effects were higher with increasing Pb concentration in unplanted soil, no matter the values of S, U, or Z. The effects of S. integra planting promoted the values of S, U, and E to be, in general, higher in rhizospheric soil than bulk soil. However, the trends were complex when compared with unplanted soil.
When the interactive effect of PbC × DS was not significant, we just compared the significant main effect, for example E and H' (Figure 7b,e). The values of E and H in rhizospheric soil were higher than in bulk and unplanted soils. Neither PbC nor DS had a significant effect on the values of D.

Relationships of Microbial Quantities, Metabolic Profiles and Soil Chemical Properties
The relationships of soil properties relative to microbial quantities and metabolic profiles were examined by RDA analysis (Figure 8a,b). Eigenvalues of RDA showed that axes 1 and 2 explained 44.5% and 14.5% of the variance of the microbial quantities data, respectively (Figure 8a). CEC values were significant (F = 4.169, p = 0.012) in explaining the variations of microbial quantities, and could explain up to 49.6% of the variation. The contributions of other factors to the observed variation were as follows: AN (22.5%), AK (10.6%), pH (6.4%), and AP (10.8%). NB had a significant negative correlation with AP (p = 0.045), and had a significant positive correlation with AN (p = 0.015). IPB and NFB had significant positive correlations with CEC (p = 0.015) and AK (p = 0.033), respectively.
As to the metabolic activity, the eigenvalues of RDA showed that RDA1 and RDA2 only explained 15.6% and 2% of the varieties of soil microbial carbon utilization rate profiles (Figure 8b). The soil properties made no significant contribution to the variation of metabolic activity. CG was negatively correlated with all factors, and had significant correlation with AK (p = 0.007) and CEC (p = 0.023).

Effects of S. integra Planting and Pb Pollution on Soil Properties and Microbial Quantities
The four kinds of bacteria were all increased in LT and decreased in HT compared with CK, especially in rhizospheric soil. A higher Pb concentration suppressed microbial growth, such that certain sensitive taxa were extinguished by heavy metal poisoning. For example, the relative abundance of Firmicutes decreased from 88.5% (CK) to 12.1% (Pb pollution treatment) during phytoremediation by Trifolium repens L. [48]. Zhu et al. [48] also suggested heavy metal contamination could play a more important role in shaping bacterial diversity. The bacterial growth rate and biomass in soils was found to consistently diminish along the Cd gradient, regardless of planting Sedum alfredii [49]. Similarly, Golebiewski et al. [50] found that species richness, diversity, and bacterial communities in soils with low levels of Pb contamination were higher than in highly polluted soils. This study found that microbial quantities were increased in lower Pb concentrations, which conflicts with former studies. We inferred that the root exudates were stimulated by a lower Pb concentration, which attracted more bacteria in return. Our study concerned the culturable bacteria only, and because of the limitation of the cultivation approach for enumeration [51], it is not guaranteed that these results can be generalized to non-cultivable microbes.
In general, AP and AK were decreased with the increase of Pb concentrations, while AN was increased. There are two reasons for this, with the first being that the chemical effect of adding Pb had significant effects on ion exchange in the soil. For example, Pb 2+ reacted with PO 3 4− to form Pb 3 (PO 4 ) 2 (insoluble matter), especially in alkaline soil, which led to the decrease of AP [52]. Another reason may be that heavy metal pollution induced changes in the microbial community structure [53][54][55] and thereby affected nutrient cycling, especially available phosphorus and potassium [50,[56][57][58][59]. However, the trends of microbial quantities and nutrients in response to Pb concentration were different in our study. The reason for this may be that the nutrient cycles required microbes with different functions, for example, the microbes involved in the nitrogen cycle, including nitrogen-fixing bacteria, ammonifying bacteria and nitrifying bacteria. Consistent with the RDA analysis, the explanatory variables of AN, AP, and AK had a relatively small impact on microbial quantities. The reason for this may be that we just measured four types of microbes, which were only a small proportion of the microbes involved in N and P cycling. At the same time, adding Pb 2+ could change other soil properties including redox potential and cation/anion exchange capacity [60]. It partly removed cations (Na + , K + , Ca 2+ , Ca2+, etc.) from the soil solution that led to the increase of CEC [60,61]. However, the CEC was decreased in the HT, which may be related to the balance of soil properties to the higher Pb pollution. In general, S. integra planting increased the values of soil properties and microbial quantities to a greater degree in rhizosphere than bulk soil, across all Pb treatments. However, the trends were complex when compared with unplanted soil. These were consistent with previous studies' indications that microbial biomass clearly increased in rhizospheric soils compared to bulk soil [16,[62][63][64]. Similarly, T. caerulescens growth increased microbial respiration rates in the Zn-and Cd-contaminated soil [23]. Planting of R. pseudoacacia and Populus also had a significant influence on the rhizospheric bacterial community structure in heavy metal contaminated soils [11]. These findings could be partly explained by the reduction of metal concentration in the rhizospheric soil due to plant uptake. Yang et al. [62] reported that phytoextraction by Sedum alfredii significantly reduced the available and total Cd of 0-2 and 2-4 mm soils. In addition, there were more root exudates in the rhizosphere environment, which may lead to higher microbial populations and stimulate the growth of specific bacterial and fungal populations [15][16][17]. Furthermore, soil organic carbon sources excreted by roots were reported to be vital in determining rhizospheric soil microbial compositions [65]. The values of soil properties that were increased also contributed to the increase of root exudates and microbial quantities in rhizospheric soil. In contrast, Epelde et al. [23] reported lower values of total N and extractable K + observed in planted versus unplanted soils because of plant uptake. The root exudates, such as carboxyl, phenolic, hydroxyl, and carbonyl groups, led to negative charges on the soil surface, and cations could be partly removed from the soil solution, thereby increasing the CEC in the rhizospheric soil [60,66].
The difference between bulk soil and unplanted soil seen with S. integra planting not only had direct effects on rhizospheric soil, but also had indirect effects on bulk soil. As in the earlier studies, the effects of the rhizosphere should not be limited to specific volumes [62,67]. For instance, the rhizosphere effects of Sedum alfredii on the microbial community structure and metabolic function were not limited to the root surface, notwithstanding that the effect decreased as the distances increased [62]. Similarly, the interactive effects of pioneer plants and harsh soil environmental conditions remodel the specific bacterial communities, both in rhizosphere and bulk soil, in mine tailings [19]. The reasons for this may be that the root exudates could be transported from rhizospheric soil to bulk soil by active transport (animal and fungal mycelium) or passive transport (concentration difference and water action) [68].

Effects of S. integra Planting and Pb on Microbial Metabolic Profiles
S. integra planting promoted microbial metabolic activities and diversity in rhizospheric soils, while Pb pollution suppressed microbial metabolic activity and diversity. The inhibitory effect was increased with increasing Pb concentration, with the exception of LT in rhizospheric soil. This finding was consistent with the results regarding microbial quantities in our study. Many researches have indicated that excessive metals in soil are toxic to microbes, and that the ecological functions performed by microbes could be seriously inhibited [69][70][71].
The soil microbial community showed a unique pattern of carbon utilization in rhizospheric soil samples. Bulk and unplanted soil samples showed carbon utilization relative to the contamination levels. The microbial community increased the utilization of AG and MG in rhizospheric soil that underwent low Pb treatment, and increased the utilization of CG under higher Pb pollution. Therefore, we inferred that the dominant species of the microbial community changed. When subjected to heavy metal stress, some microbial populations become extinct, while certain resistant community members survive and form the basis of a new community [72]. Plant-microbes-environment balance and interact with each other during phytoremediation. This process is very complex.
Some researchers have shown that the composition and quantities of root exudates change a lot when plants respond to heavy metal stress [60,73,74]. Plants can be selected for beneficial microbial communities in their rhizosphere by releasing certain metabolites targeting specific microbes [75]. AG and CG were the major components of root secretion. Amino acids and organic acids, which are released from plant roots, have been shown to enhance phytoextraction of heavy metals by increasing their bioavailability [75][76][77][78]. Alanine and proline activated Pb and increased the Pb uptake ability of a Pb-accumulating ecotype of Sedum alfredii [79]. The literature about the effect of carbohydrates is relatively sparse, however, it is suggested that they could immobilize Pb 2+ and Cu 2+ and prevent them entering into the plant [80]. Dissolved organic matter, which is released from the rhizosphere of S. alfredii, could markedly decrease Zn and Cd adsorption and enhance their mobility by forming soluble organic matter-metal complexes [81]. Of course, detailed information about the components of root exudates from S. integra require further study to reveal the mechanism for regulating microbial metabolism activity.
Environment factors only explained 14.8% of the variance of soil microbial carbon utilization rate profiles. Only AK and CEC were significantly negatively correlated with CG, which further shows that the environmental factors determined by us are not the main reason for the metabolic differences of the microbial community. Therefore, root exudates and other factors, including heavy metal content, should be taken into consideration.

Conclusions
S. integra planting increased the microbial quantities and metabolism activity in rhizospheric soil. The effect of S. integra growth was observed not only in the rhizosphere, but also in bulk soil, although the range of influence requires further study. Microbial quantities and metabolic patterns were relative to the Pb pollution levels: a low concentration increased the microbial quantities and a high concentration decreased it, to different degrees. The microbial metabolisms were decreased with the increase of Pb levels, with the exception of the rhizospheric soil of LT. The utilization of CG was decreased, and of AG and FG was increased in the low Pb treatment, while the trends were opposite to this in the high concentration condition. Therefore, we believe that S. integra develops its own strategy to survive in different levels of polluted soil by changing root exudates, and in doing so further regulates the microorganism community. In return, microbes activate different defense mechanisms to eliminate Pb stress. Therefore, advancing the current understanding of the interaction of the root-microbial community will be essential in the future, for example by conducting metabonomic research to analyze root exudates and microbial metabolism by Gas Chromatography-Mass Spectrometer (GC-MS) or nuclear magnetic resonance (NMR).


text:
Contribution of Rainwater to the Irrigation Requirement for Paddy Cultivation at Tanore Upazila in Rajshahi, Bangladesh

Groundwater extraction for irrigation in Bangladesh has caused groundwater depletion, especially in the Northern region. As such, shifting reliance from groundwater to surface water is one of the solutions to mitigate this problem. This study aims at investigating the contribution of effective rainfall to the total consumptive use requirement of rice cultivated in Tanore, Rajshahi, in Bangladesh. The prospect of rainwater harvesting using ponds and its contribution to the consumptive use requirement of rice was also determined. Effective rainfall, temperature, monthly percentage of bright sunshine hours, and consumptive use factor for rice data were collected from Bangladesh Meteorological Department (BMD) and Barind Multipurpose Development Authority (BMDA). Blaney-Criddle method was used to calculate total consumptive use requirement of rice. Analysis results showed that from June to September, there is no requirement of irrigation due to high rainfall in this period. It was revealed that total consumptive use requirement of rice in Tanore averaged at 436 million m3 of water between the period of 2005 and 2012. Effective rainfall contributed to about 38% of the total consumptive use in this period. It was also found that another 5% of the total consumptive use requirement could be supplied by rainwater harvesting using the ponds of Tanore.

Introduction
Bangladesh is one of the most densely populated countries in the world and it had a population of about 120 million in its 147 570 km 2 (14.57 million ha) of land in 2011 1 and now the population has increased to 160 million. The country is suffering from food shortage. As a result, agricultural land is cultivated more than two to three times in a year requiring more irrigation water than the past. With the increasing use of groundwater for irrigation, the annual extraction of groundwater is far higher than the net average recharge from natural sources. It can lead to lowering of the groundwater table 2 increasing the cost in extraction of groundwater. Future crop production may be hampered by this depleted groundwater table. 3,4 Aeschbach-Hertig and Gleeson 2 collected and showed the evidence that groundwater across the world is on decline. Excessive extraction for agricultural use has been identified to be the main reason behind this phenomenon. 5 Konikow 6 confirmed that United States is also experiencing groundwater depletion. It has been confirmed by plenty of studies that groundwater across the globe is depleting at an alarming rate. 5,[7][8][9][10] It has been reported that the rainfall amount in the northeastern region of Bangladesh is progressively declining. 11 Alam et al 12 warned that western part of Rajshahi district in Bangladesh might experience chronic to severe drought in the future. It has been found that groundwater in Rajshahi district was receding at a rate of 0.23 m/year during the period between 2000 and 2014. 13 Moreover, a study on Tanore Upazila in Rajshahi revealed that groundwater in Tanore is depleting at an average rate of 0.22 m/year. 14 The impact of climate change on water requirement of rice was investigated by Shahid. 15 It was concluded that paddy crop itself will not require extra amount of water, but the warmer climate will require more water to be pumped, which will ultimately aggravate the already stressed groundwater resources of Bangladesh. Mainuddin et al 16 estimated that irrigation water requirement can increase by 1% to 5% in 2050 based on type of crop and time of sowing. To reduce this problem, the excess rainwater yielded as surface runoff can be harvested in ponds and other reservoirs and may contribute to the consumptive use of crops in dry season.
A study conducted in the South-Eastern part of Bangladesh revealed that about 13.94% to 100% of the total consumptive use requirement could be fulfilled by effective rainfall for different parts of the year for rice cultivation. 17 The study used Potential Evapotranspiration/Precipitation Ratio method, United States Department of Agriculture (USDA), Soil Conservation Service (SCS) method, and Renfro equation to estimate the effective rainfall percentage of the total consumptive use requirement. A recent study reconfirmed the depletion of groundwater table decline in Rajshahi district and suggested potential measures to release the stress on groundwater. One of the suggestions was to Boro rice transplantation by November so that late monsoon rain can contribute to the consumptive use requirement. 18 Another study investigated the opportunity of rainwater harvesting (RWH) in Godagari, Rajshahi. 19 The study showed that in Boro season, rice cultivation was fully dependent on irrigation water. It also depicted that there were 15 ponds located in Godagari, which could be used to store runoff and supply irrigation water during dry period. Water requirements for various crops in the Barind region of Bangladesh were determined using the available data from 1985 to 2013. Blaney-Criddle formula was employed to determine the consumptive use requirement of the crops. During this period, the maximum value of consumptive use was found to be 27.88 cm for rice in April.
In Faisalabad, Pakistan, study showed that the rainwater contribution for wheat and maize was around 7.5% and 15.5% per year. 20 The study used CROPWAT computer model to determine consumptive use requirement in that region. Wisser et al 21 performed a global scale study on quantifying the significance of water captured by small reservoirs for crop production. It was reported that 623 to 1122 km 3 of water can be supplied for irrigation worldwide through the water captures in local small reservoirs. An increase of 35% in global cereal production is expected due to this supply of irrigation water. However, it was also mentioned that 20% of the water stored in the reservoirs can get evaporated before use. Another study suggests that climate change will reduce the global crop production by 9%, which can be balanced by harvesting 25% of the rainwater to supply for irrigation. 22 The study suggested to adopt a combination of water management strategies to enhance crop production. A model was developed to assess hydro-economic benefit of RWH system for irrigation in India. 23 Impact of on-farm reservoir (OFR) irrigation system was simulated and compared with the rainfed agriculture. It was revealed that the water use efficiency of the OFR system was 82% more efficient than the rainfed system. Total gain value due to the OFR irrigation system was estimated to be 74% and 14% greater than the rainfed system in the dry and wet season, respectively. Future water requirement and availability in terms of blue and green water was estimated by Rockström et al. 24 Green water is the soil moisture obtained through rainfall used by the plants and blue water is the water available in the surface and ground. It was reported that in 2050, 36% of the world population will have to tackle green and blue water shortage and 59% of world population will encounter blue water scarcity. It was also mentioned that currently, there are many countries who are blue water short, but they have sufficient green water, which combined with blue water is sufficient to produce enough food. It was suggested that proper management of local green water resources can increase resilience to climate change. A countrywide study on Pakistan revealed that effective rainfall contribution to consumptive use requirements of crops grown in Pakistan varied from 13.03% in the Rabi season (October to April) and 21.31% in the Kharif season (May to September) to almost 100% in different parts of Pakistan. 25 The study used data from 58 meteorological stations across Pakistan covering irrigated plains. Among the methods employed by the study involved Renfro equation, US Bureau of Reclamation method, USDA method, and SCS method. Effective rainfall contribution is expected to increase in the future, especially in the Northern part of Pakistan. 26 On the other hand, in case of drought event, severity is going to be intensified in the arid and semi-arid region of Pakistan. 27 Under this circumstance, storing of water for the dry spells may prove to be a viable option for smooth water supply. Djaman et al 28 studied the water requirement of rice in Senegal. It was found that based on location and season, water requirement of rice in Senegal varied between 632 and 929 mm. 28 On the other hand, rainfall amount in Senegal is less than 190 mm, which accounts for 20% to 30% of the total water requirement of rice.
As mentioned earlier, northern part of Bangladesh is suffering from declining rainfall and depleting groundwater table. Hence, this region needs immediate attention and research should be executed to solve the water problem in this region. With an aim of addressing the issue, this study has chosen Rajshahi district. The approach of the study is to find out the water requirement of the main crop, ie, paddy produced in this region, and find alternative solution to meet the water demand of the crop. There are studies on global scale as well as countrywide studies. However, studies on micro scale or localized studies are rare. Localized studies are more important because rainfall amount and crop water requirement for the same crop vary from region to region. This study aimed at investigating the contribution of effective rainfall to the total consumptive use requirement of rice in Tanore Upazila, Rajshahi, Bangladesh. Storing runoff water by RWH and supplying in the dry season is a prospective solution to the groundwater depletion problem in the region. Hence, the contribution of rainfall runoff through RWH is also estimated in this study.

Study area
The study area is located in the Barind tract situated in the Northern part of Bangladesh (Figure 1). Tanore is a part of greater Rajshahi district. The co-ordinates of Tanore are 24° 29′ N to 24° 43′ N in latitudes and 88° 24′ E to 88° 38′ E in longitude. Population of Tanore is about 173 495. Total area of Tanore is 295.8 km 2 (29580 ha). Total cultivable land is about 226.56 km 2 (22 656 ha). Nearly 43% of the people are involved in agricultural activities. With the expansion of irrigation in Tanore Upazila in the mid-1980s, a revolutionary change occurred in the agricultural sector of Tanore. Introduction of high-yielding varieties (HYV) paddy combined with groundwater supply for irrigation made it possible to increase the yields of rice. Single cropped agricultural lands are now producing three crops in one agricultural year after the start of groundwater extraction. The cropping intensity in Tanore is about 262% which is high compared with the national average of 180%. 1 Rice in Bangladesh is grown in three distinct seasons, namely Boro ( January to June), Aus (April to August), and Aman (August to December). 29,30 Estimation of consumptive use and field irrigation requirement Evaporation (E T ) or consumptive use (U) is the total amount of water lost from a cropped or irrigated land due to evaporation from the soil and transpiration by the plants or used by the plants in building up of plant tissue. It is usually expressed as a depth (cm, mm) over the area. 3 Blaney-Criddle formula was used to determine the consumptive use requirement of rice. 31,32 Considering the available data, this method was used. Other methods, eg, Penman equation could not be used because of lack of required data. The Blaney-Criddle formula is as follows: where U is the seasonal consumptive use (cm), t is the mean monthly temperature (°C), p is monthly percentage of hours of bright sunshine (of the year), k is the monthly consumptive coefficient, and f is the monthly consumptive use factor The field irrigation requirement is the amount of water to be supplied by the crop by irrigation; it is determined by subtracting the effective rainfall from the crop consumptive use requirement.
Effective rainfall is the amount of rainwater generated as runoff after deducting the infiltration loss. To determine effective rainfall, a runoff coefficient (C), which is below 1, is multiplied to the actual rainfall depth. The value of runoff coefficient depends on soil type. Sandy soil will have lower value of C and clayey soil will yield higher C value. Hence, P e can be determined as follows: where C is the runoff coefficient and P is the actual depth of rainfall (cm).
The value of C was considered to be 0.8 for Tanore as suggested by Barind Multipurpose Development Authority (BMDA). Field irrigation requirement (FIR) is the difference between consumptive use and effective rainfall: where FIR is the field irrigation requirement and P e is the effective rainfall.

RWH calculation
In estimating the volume of water from RWH from the areas where land is not cultivated, rational formula was used as follows: Q CAP = where Q is the volume of water from RWH yield (m 3 ), C is the coefficient of runoff of 0.8 for Tanore (according to BMDA), A is the land area not cultivated (m 2 ), and P is the rainfall depth (m). The amount of cultivated land in Tanore is 22 656 × 10 4 m 2 (22 656 ha). The remaining land which is not cultivated is 6924 × 10 4 m 2 (6924 ha) and hence, A = 22656 × 10 4 m 2 .

Data collection
Rainfall data, monthly average bright sunshine hours, and temperature data of the study area were collected from the Bangladesh Meteorological Department (BMD). Number of ponds, pond dimensions, total cultivable area of the study, and consumptive use factor for rice in different months were collected from the BMDA.
Available dataset for the meteorological parameters mentioned above ranges from 2005 to 2012. The data were monitored, sorted, and bias corrected, ie, the quality of the data was maintained by the BMD. It records three hourly data for rainfall and hourly data for temperature. Then the data are bias corrected and averaged for monthly period. The monthly averaged data were provided by the BMD. There is no missing data  Figure 2 shows temperature variation throughout the year from 2005 to 2012. December and January are the winter season and the temperature is lowest in these months. On the other hand, summer persists in April and May and the temperature rises near 40 C in these months. The temperature never went below 20 C even in the winter season, which is not the typical case in Bangladesh. In most part of the year, temperature in Tanore remained above 30°C.

Results and Discussions
Rainfall variation in different months of the year is depicted in Figure 3. June and July are considered monsoon in Bangladesh and monthly total rainfall depth up to 49.2 cm occurred in this period. Tanore Upazila experiences the lowest amount of rainfall during winter season (December and January) and the trend starts in November and continues until February. This has great implication on rainwater contribution to the irrigation requirement, ie, the winter season can be considered as the driest period in a year and irrigation water requirement will be high. Monthly total rainfall depth varied from 0 cm in winter to 49.2 cm in monsoon. Almost 92% of the total rainfall occurred between May and October in the period of 2005 to 2012.
Monthly percentage of bright sunshine hours (p) measured by BMDA are illustrated in Figure 4 in the Whisker-Box plot. The value of p is lowest in July -median value is 4.4% due to heavy rainfall in this month. The median value of p in August is 5.2%, which also experiences heavy rainfall. The months of December and January, which are in the winter season also have low p values -medians are 4.8% and 5.6%, respectively. As p value is directly proportional to the consumptive use requirement, a lower p value will reduce the consumptive use requirement.    December -0.95. A higher k value will yield a higher consumptive use requirement. Figure 6 demonstrates variation of monthly consumptive use. Due to low temperature, consumptive use was lowest in December. Lower percentage of bright sunshine hours in July helped in lowering the consumptive use requirement. On the other hand, consumptive use was highest in April and May due to high temperature in this period. The consumptive use requirement varied between 7.7 and 26.7 cm.
Field irrigation requirements after deducting the effective rainfalls from the consumptive uses are given in Figure 7. It can be clearly observed that between June and September, rainfall amount is excessive than the consumptive use requirements and hence, FIR is negative. On the other hand, from February to April, FIR is high -average is near 18 cm. Although consumptive use is not high in the winter, low amount of rainfall keeps the consumptive use and FIR almost same -around 10 cm. From the analysis, it can be stated that lands are needed to be irrigated in 8 months of a year. Consumptive use requirement of the remaining 4 months can be fulfilled by rainfall alone.

Existing water reservoir capacity
In the study area, there are total of 34 existing ponds having a total volume of 432 753 m 3 . These ponds can be used for storing rainwater and runoff and supply water for irrigation in other times of the year. Pond dimensions and water reservoir capacity of these ponds are shown in Table 1. Figure 8 shows the percentage of consumptive use requirement that rainwater contributed from 2005 to 2012 and the rainfall runoff water that could be used to supply for irrigation purpose in non-raining season. The total consumptive use varied between 411 and 450 mm 3 of water per year and field irrigation water requirement was between 164 and 304 mm 3 per year. Consumptive use water requirement and field irrigation water requirement averaged at 436 and 213 mm 3 of water per year, respectively. Effective rainfall contributed 33% to 46% of water to the consumptive use requirement of rice between 2005 and 2012. The average contribution of effective rainfall was about 38%, whereas in Faisalabad, Pakistan, rainfall accounted   20 On the other hand, for other parts of Pakistan, rainfall contribution was found to be between the range of 13.03% and 100%, depending on the season and location. 25 In Senegal, it was found that rainfall contributes less than 30% of the evapotranspiration requirement of rice. 28 The stored runoff water could supply 4% to 9% of consumptive use requirement. The average water that could be stored from runoff and supplied for irrigation was 22 mm 3 , which accounts for roughly above 5% of total consumptive use. On an average, the effective rainfall and stored runoff in the reservoir ponds could supply 43% of the total consumptive use requirement between 2005 and 2012 in Tanore Upazila.

Conclusions
This study has shown the contribution of effective rainfall and potential contribution of RWH in total consumptive use requirement of rice in Tanore, Rajshahi. Monthly effective rainfall depth in Tanore varied between 0 and 49.2 cm. It has been revealed that from June to September, there is no requirement of irrigation water as amount of rainfall exceeds the consumptive use requirement in these months. Effective rainfall contributed from 33% to 46% of total consumptive requirement of rice between the period of 2005 and 2012. Another 4% to 9% of consumptive use requirement can be met by storing rainfall runoff water in ponds located in Tanore Upazila. It is expected that irrigation water supply from stored runoff in the ponds if implemented can release the stress on groundwater at least to some extent. Now, the managers can take the decision of implementing the RWH in local ponds and increasing the capacity of local reservoirs by constructing more ponds. This will not only help solving the groundwater depletion problem, but also enhance resilience of irrigation water supply system during dry spells and benefit in the long run. Researchers can take this study as an example to quantify effective rainfall contribution to the crop water requirement as well as potential of RWH in other regions of Bangladesh, especially in groundwater depletion vulnerable areas. This study used only one method to estimate water requirements of rice in the region. We suggest using other methods too in estimating water requirements of rice and other crops. It will help to pinpoint the exact water requirement of crops and reduce errors. Future research can focus on variability of climate in the region and crop water requirements under those conditions.

Author Contributions
MRH checked for the accuracy of the data, analyzed the reservoir capacity and rainwater harvesting potential and collected detailed information on the study area. MN calculated the consumptive use requirement, rainwater contribution and developed the overall methodology of the work. AAM analyzed the rainfall data, assisted in developing the methodology and drawing the conclusions.