Abstract
Lodging is one of the major abiotic stresses, affecting the total crop
yield and quality. The improved lodging resistance and its component
traits potentially reduce the yield losses. The section modulus (SM),
bending moment at breaking (M), pushing resistance (PR), and
coefficient of lodging resistance (cLr) are the key elements to
estimate the lodging resistance. Understanding the genetic architecture
of lodging resistance–related traits will help to improve the culm
strength and overall yield potential. In this study, a natural
population of 795 globally diverse genotypes was further divided into
two (indica and japonica) subpopulations and was used to evaluate the
lodging resistance and culm strength–related traits. Significant
diversity was observed among the studied traits. We carried out the
genome-wide association evaluation of four lodging resistance traits
with 3.3 million deep resolution single-nucleotide polymorphic (SNP)
markers. The general linear model (GLM) and compressed mixed linear
model (MLM) were used for the whole population and two subpopulation
genome-wide association studies (GWAS), and a 1000-time permutation
test was performed to remove the false positives. A total of 375
nonredundant QTLs were observed for four culm strength traits on 12
chromosomes of the rice genome. Then, 33 pleiotropic loci governing
more than one trait were mined. A total of 4031 annotated genes were
detected within the candidate genomic region of 33 pleiotropic loci.
The functional annotations and metabolic pathway enrichment analysis
showed cellular localization and transmembrane transport as the top
gene ontological terms. The in silico and in vitro expression analyses
were conducted to validate the three candidate genes in a pleiotropic
QTL on chromosome 7. It validated OsFBA2 as a candidate gene to
contribute to lodging resistance in rice. The haplotype analysis for
the candidate gene revealed a significant functional variation in the
promoter region. Validation and introgression of alleles that are
beneficial to induce culm strength may be used in rice breeding for
lodging resistance.
Keywords: GWAS, association mapping, rice (Oryza sativa L.), lodging
resistance, culm strength, genetic architecture
1 Introduction
Rice (Oryza sativa L.), as one of the main staple food crops, is very
important for food security worldwide ([52]Meng et al., 2021). Among
the various biotic and abiotic stresses, lodging is a major factor,
causing significant reduction in grain yield and quality by decreased
photosynthesis and nutrient transportation ([53]Berry et al., 2004;
[54]Islam et al., 2007). The modern high yielding varieties bred for
long spikes and high grain weight become vulnerable to stem lodging
([55]Hirano et al., 2017). The field management with more fertilizer
consumption and high planting density leads to reduced culm strength
and an increased risk of lodging ([56]Shah et al., 2019). Previously,
reduced plant height with semi-dwarf gene sd1 caused a green revolution
and decreased the risk of lodging in cereals ([57]Peng et al., 1999;
[58]Sasaki et al., 2002). However, the clear genetic architecture of
lodging resistance is still unclear.
There are two main types of lodging in rice, that is, stem lodging and
root lodging. The current study focused on the traits associated with
lodging resistance of rice culm. The significant association of lodging
resistance with culm strength has been reported in rice ([59]Meng et
al., 2021). Culm strength is one of the complex quantitative traits
that are contributed by multiple traits, including culm morphology,
diameter, thickness, plant height, and stem stress bearing capacity,
due to accumulation of macromolecules like lignin and carbohydrates
([60]Kashiwagi and Ishimaru, 2004; [61]Kashiwagi et al., 2006; [62]Yano
et al., 2015; [63]Fan et al., 2018; [64]Sowadan et al., 2018). Culm
lodging can be further classified as culm-breaking and culm-bending
types. Both of the culm strength indicators often occur in
transplanting cultivation. The culm-breaking can be estimated by the
bending moment at breaking (M), which is the measure of culm thickness
calculated as section modulus (SM) and culm stiffness calculated as
bending stress (BS) ([65]Ookawa et al., 2010; [66]Chigira et al.,
2020). Meanwhile, the culm flexibility can also be evaluated as the
coefficient of lodging resistance ([67]Grafius and Brown, 1954).
Some studies and breeding trials have been conducted to improve the
culm strength, but they have not succeeded due to the negative
correlation of culm strength (tallness and thickness) and yield
parameters ([68]Hirano et al., 2017). Exploitation of potential genetic
resources, identification of novel loci governing the culm strength,
and gene pyramiding are the promising solution to improve yield and
strength simultaneously ([69]Meng et al., 2021). Recently, various
genes/QTLs governing lodging resistance and/or its contributing traits
have been reported. The QTLs SCM1 and SCM2 were reported to induce culm
strength. SCM2, similar to the APO1 gene, controls the culm diameter
and culm morphology by encoding an F-box domain containing protein
([70]Ookawa et al., 2010; [71]Rashid et al., 2016). A QTL SCM3 was
reported to control the culm morphology by affecting the strigolactone
signaling pathway ([72]Yano et al., 2015; [73]Meng et al., 2021). The
gnla mutant could improve the lodging resistance by increasing culm
diameter ([74]Tu et al., 2022). The COBRA protein encoding gene BC1 was
reported to contribute the culm strength by cell wall cellulase
microfilaments elongation and cell wall thickening ([75]Li et al.,
2003). The sd1 gene and TUT1 allele es1-1 were reported for short
stature, resulting in lodging resistance. In addition to these genes,
the genetic studies on internode length, grain weight, panicle length,
and tiller angle ([76]Rashid et al., 2022) have also been observed to
contribute to culm strength, leading to lodging resistance in rice.
In recent years, rice researchers widely used the genome-wide
association studies (GWASs) to investigate the complex traits of
agronomic and commercial importance including yield parameters,
flowering time ([77]Huang et al., 2010; [78]Huang et al., 2012;
[79]Rashid et al., 2022), and abiotic stresses like drought
([80]Al-Shugeairy et al., 2015) and salinity ([81]Kumar et al., 2015).
Nonetheless, a few GWAS studies were conducted for lodging resistance
([82]Hu et al., 2013; [83]Meng et al., 2021) in crops. In this
research, a genome-wide association study (GWAS) was conducted for the
SM, M, PR, and cLr. The regression models, such as general linear model
(GLM), mixed linear model (MLM), and compressed mixed linear model
(cMLM), were used to estimate the marker trait association (MTA). The
1000-times permutation test was also used to remove false-positives.
Moreover, the combination of various regression models in GWAS helped
to control false positives ([84]Wu et al., 2016; [85]Misra et al.,
2017; [86]Zhang et al., 2018). This study clearly aimed to explore the
genomic regions associated with lodging resistance in the whole rice
genome. We investigated the genetic variation and heritability among
lodging resistance traits, their statistical and biological
correlation, significantly associated loci, pleiotropic loci, and
candidate genes. The research will not only provide a deep insight into
the genetic basis of lodging resistance but also point toward improving
yield potential in rice breeding programs.
2 Materials and methods
2.1 Plant material and field management
The study was conducted by using a natural population of 795
genetically and geographically diverse rice genotypes selected from a
full population of 3,000 rice genome projects (3KRGP) ([87]Li J. Y. et
al., 2014; [88]Li Z. et al., 2014). The germplasm was selected as a
representative population including the rice mini-core collection from
China Agricultural University (CAU) China and 525 diverse breeding
lines from the Chinese Academy of Agricultural Sciences China as
previously described ([89]Zhao et al., 2018a; [90]Zhao et al., 2018b;
[91]Rashid et al., 2022). The genetic purity and homozygosity of
germplasm was confirmed by growing the plants through a single seed
descent (SSD) approach for 3 years (2013–2016). The experimental plant
material was grown in a randomized complete block design (RCBD) in
three replications for 3 consecutive years. Two experimental blocks at
the bird-net secured experimental farm of China Agricultural
University, Sanya, Hainan (18 ^oN, 109 ^oE), China. Three lines of 1 m
length for each genotype were grown, and seven seeds per line were
maintained. Sampling and phenotyping were performed on three guarded
plants from each line. The average value of three replicated plants
from each block during 3 years ([92]Yu et al., 2007) was used for
analysis. The standard cultural practices as per local area
requirements were kept constant.
2.2 Phenotyping
The lodging resistance of three guarded plants from each of two
experimental blocks was calculated by the standard lodging resistance
parameters, including section modulus (SM), pushing resistance (PR),
bending moment at breaking (M), and coefficient of lodging resistance
(cLr). SM was measured in cubic-millimeters (mm^3) by the formula
([93]Ookawa et al., 2010):
[MATH:
SM=π⋅{a13b1−a2
3b2}<
/mo>32a1, :MATH]
where a[1] is the culm outer diameter of the minor axis in an oval
cross-section, b[1] is the outer diameter of the major axis in an oval
cross-section, a[2] is the inner diameter of the minor axis in an oval
cross-section, and b[2] is the inner diameter of the major axis in an
oval cross-section. The whole plant’s pushing resistance (PR), also
known as bending stress at breaking (BnS), was measured in grams per
centimeter-square (g/cm^2) with a lodging meter at the maturity stage,
according to a method reported previously ([94]Kashiwagi and Ishimaru,
2004). The bending moment (M) of the basal internode at breaking (g/cm)
was calculated as follows as previously described ([95]Ookawa et al.,
2016):
[MATH:
M=Section <
mtext>Modulus (SM)×Bending Stress (BnS). :MATH]
The coefficient of lodging resistance (cLr) value of the plant was
determined in g/cm according to the method described by Grafius and
Brown ([96]Grafius and Brown, 1954). The average value of ten
individual plants was considered the final reading for the succeeding
analysis.
2.3 Germplasm genotyping and population structure evaluation
The second generation of whole-genome high-throughput sequences with
∼13-fold coverage was obtained by Illumina sequencing, and the sequence
reads were aligned against japonica cv. Nipponbare reference genome.
The direct comparison with corresponding regions in the reference
sequence produced 10 million raw SNPs. The SNP screening was performed
with threshold criteria of missing rate more than 20 percent and at
least 5% minor allele frequency. The completed set of selected SNP
markers distributed over the 12 chromosomes of the rice genome were
used to calculate the population structure parameters.
The population structure (Q-matrix), principal components (PCs) and
kinship (K) matrix were calculated by the GAPIT program
([97]http://www.maizegenetics.net/GAPIT). The 795 individuals in the
full population could be divided into two sub-populations as japonica
and indica based on PCs and K-matrix. Later on, the kinship matrix and
the principal components were re-estimated for each sub-population to
be used as covariates in the mixed linear model (MLM) of regression
analysis. To map the known genes and candidate loci, their physical
positions were acquired from the online available rice database
([98]http://rapdb.dna.affrc.go.jp).
2.4 Statistical analysis of phenotypes
The statistical descriptive for phenotypic or morphological data was
calculated by SPSS version 19
([99]http://www-01.ibm.com/software/analytics/spss/). The
F-distribution test was performed to evaluate the significance of
available diversity among genotypes. The analysis of variances (ANOVA)
was performed to estimate the genotypic variance (V[g]) and
environmental variance (V[e]). The interaction variances were estimated
by using the STATISTICA software (StatSoft 1995; Tulsa, OK, United
States). The experimental repeatability (r ^2) and the broad-sense
heritability (h^2) ([100]Sukumaran et al., 2015) were estimated by
using genotypic variance (V[g]), error variance (V[e]), and interaction
variances as previously described ([101]Rashid et al., 2022). All the
variance components were estimated by ANOVA in the F-distribution test.
The frequency distribution–based data normality was confirmed by the
Q-Q plot drawn by the SPSS program and GAPIT scripts in the R-program.
2.5 Genome-wide association analysis
The average value of ten plants for each genotype was used as the final
estimate of each trait. For the marker trait association for lodging
resistance related traits, two regression models, viz, general linear
model (GLM) and compressed mixed linear model (cMLM) were applied to
three sets of populations as (i) whole population, (ii) indica
sub-populations, and (iii) japonica sub-populations. In the GLM model,
the population structure (PCAs) was incorporated as a covariate. The
population structure (Q-matrix) and kinship matrix (K) were included as
fixed and random effects, respectively, to conduct a compressed mixed
linear model (cMLM). The 1,000-permutation test by the mixed linear
model MLM was performed for detection and screening of false SNP
signals.
2.6 Identification of significant QTLs, candidate genes, and pleiotropic loci
The QTL regions were strictly defined by linkage disequilibrium (LD)
blocks of 167 kb and 123 kb left and right genomic regions of
significant SNPs (SNP ± LD) for japonica and indica subpopulations,
respectively, and 170 kb for the whole population, as previously
reported ([102]Huang et al., 2010; [103]Huang et al., 2012; [104]Zhao
et al., 2018a; [105]Zhao et al., 2018b). To estimate a significance
threshold level, we used the formula “−log10 (0.01/effective number of
SNPs with a p-value less than 0.01),” that is, the threshold at a
significance level of 1% after Bonferroni-adjusted multiple test
correction ([106]Zhao et al., 2018b). The QTLs with no annotated loci
were removed. The significant QTLs were further screened by deletion of
non-coding SNPs residing other than CDS or promoter regions of
annotated loci. The significant loci identified for more than one
lodging resistance parameter were considered pleiotropic loci.
2.7 Functional annotation and pathway enrichment analysis
All the candidate loci within the QTL range of pleiotropic loci were
subjected to detect their functional annotation. The candidate genes
list was submitted to ShinyGo 0.76 ([107]Ge et al., 2019)
([108]http://bioinformatics.sdstate.edu/go/) for gene ontological
analysis and the same list was submitted to the Database for
Annotation, Visualization, and Integrated Discovery ([109]DAVID)
([110]Sherman et al., 2022) for functional pathway evaluation based on
the [111]Kyoto Encyclopedia of Genes and Genomes (KEGG) ([112]Kanehisa
et al., 2020) pathway database.
2.8 In silico and in vitro expression estimation
The in silico expression profile of top candidate genes was further
evaluated in the rice gene expression database (RED)
([113]http://expression.ic4r.org/). Total RNA was isolated from the
stem part collected at fourth internode of the strong and weak plants
using the RNA-prep pure plant kit (TIANGEN) and reverse transcribed
using HiScript II QRT SuperMix (Vazyme) ([114]Chun et al., 2020).
Quantitative real-time PCR for selected candidate genes in a
pleiotropic QTL on chromosome 7 was performed using SYBR green mixture
on an ABI 7500 real-time PCR detection system. Rice Ubiquitin gene was
used as an internal control for normalization ([115]Chun et al., 2020).
All primer sequences are listed in [116]Supplementary Table S1.
2.9 Haplotyping of candidate genes
The total gene length of candidate genes were obtained and SNPs in the
promoter and coding gene region were extracted. All the significant and
non-significant SNPs in the promoter and CDS regions of selected loci
were used to evaluate their functional haplotypes ([117]Rashid et al.,
2016). The SNP markers were aligned and arranged in a Microsoft Excel
sheet and saved as a text file in FASTA format. The program DnaSP was
used for haplotype analysis, and the results were visualized by Prism 8
software ([118]Librado and Rozas, 2009).
3 Results
3.1 Trait diversity and heritability among lodging resistance traits
The globally diverse population of 795 genotypes was evaluated for four
lodging resistance traits, viz., SM, M, PR, and cLr. The existence and
extent of natural variation among these traits was estimated. A high
level of divergence among genotypes and a strong correlation among the
traits were observed for all studied traits. An identical and
positively semi-skewed frequency distribution was observed for all
traits in three populations ([119]Figure 1; [120]Table 1), which showed
their suitability for the genome-wide association study. The
statistical descriptive of studied traits were compared among the whole
population, japonica subpopulation, and indica subpopulation. The
ranges and average values of all studied traits for the whole
population were similar to japonica subpopulation, while their ranges
in indica subpopulation were shorter ([121]Table 1). The average value
of SM, M, and PR was higher in japonica subpopulation than that of
indica, which may indicate the greater strength of the available
japonica germplasm. A significant variation among the genotypes of the
natural population was observed by two-way analysis of variance (ANOVA)
and coefficient of variation ([122]Table 1, [123]Table 2). The
genotypes showed a variable phenotypic response in various environments
(years and blocks). There was no variation in SM among various years,
and there was no interaction of genotypes with annual ecological and
metrological conditions ([124]Table 2). The genotypes showed a
significantly variable phenotypic response in different blocks except
for PR and M, while there was a significant interaction of blocks with
genotypes. The significantly high values (0.87–0.97) of the
repeatability parameter (r ^2 ) along with broad sense heritability (h
^2 ) demonstrated the genotypes as the main source of variation for
lodging resistance, which was consistent in all three populations. The
heritability in SM and cLr was relatively low, whereas the genotype
effect ([125]Table 1) as the main cause of variation was confirmed by
the higher repeatability estimate.
FIGURE 1.
[126]FIGURE 1
[127]Open in a new tab
Frequency distribution and Q-Q plots of 795 genotypes of a globally
diverse population for four lodging resistance traits.
TABLE 1.
Statistical descriptive of four culm lodging resistance–related traits
in whole population and japonica and indica subpopulations.
Range Mean ± SD CV (%) h ^2 (%) Vg V(G x E) Ve Skewness Kurtosis r ^2
SM (mm^3) Whole 0.13–23.35 6.23 ± 3.65 58.69 23.10 135.85 11.16 5.88
1.24 2.12 0.88
Jap 0.13–23.35 6.64 ± 4.13 62.17 15.42 46.49 10.49 4.91 1.31 2.11 0.88
Ind 0.44–21.67 6.02 ± 3.37 56.06 15.24 54.31 10.91 6.32 1.09 1.48 0.89
M (g.cm) Whole 273.33–2306.67 824.9 ± 215.11 30.53 84.32 251,015.31
15,572.71 2853.91 0.65 1.36 0.97
Jap 313.33–2306.67 844.35 ± 261.07 30.92 93.7 228,633.5 13,088.87
2694.50 1.16 3.64 0.97
Ind 273.33–1650 815.08 ± 246.77 30.28 87.95 262,527.1 16,770.43 2912.38
0.33 −0.21 0.97
PR (g.cm^2) Whole 13.67–115.33 41.25 ± 12.59 30.53 55.55 912.22 261.21
54.17 0.64 1.37 0.87
Jap 15.67–115.33 42.22 ± 13.05 30.92 73.07 881.47 281.12 47.38 1.16
3.63 0.86
Ind 13.67–82.50 40.75 ± 12.34 30.28 47.59 913.64 248.83 57.69 0.33
−0.21 0.87
cLr (g/cm) Whole 0.15–0.87 0.43 ± 0.14 33.19 24.46 0.12 0.03 0.01 0.51
−0.31 0.89
Jap 0.17–0.82 0.42 ± 0.14 32.97 67.08 0.13 0.03 0.19 0.57 −0.28 0.89
Ind 0.15–0.78 0.44 ± 0.15 33.27 24.72 0.12 0.03 0.50 0.47 −0.32 0.89
[128]Open in a new tab
SD, standard deviation from mean; CV, coefficient of variation; h^2,
broad sense heritability; r ^2, repeatability; V[g], genotypic
variance; V (G x E), genotypes to environment interaction variance;
V[e], environmental variance.
TABLE 2.
Analysis of variance among the 795 genotypes of globally diverse
population for lodging resistance traits.
Trait SOV DF MS Significance
SM Genotypes (G) 794 51.524 ***
Years (Y) 2 15.488 NS
G x Y 794 5.935 NS
Blocks (B) 1 529.751 ***
G x B 794 10.756 ***
PR Genotypes (G) 794 901.352 ***
Years (Y) 2 57.641 ***
G x Y 794 57.461 NS
Blocks (B) 1 330.984 NS
G x B 794 260.496 ***
M Genotypes (G) 794 250,486.96 ***
Years (Y) 2 139,786.52 ***
G x Y 794 3283.60 NS
Blocks (B) 1 19,437.789 NS
G x B 794 15,453.52 ***
cLr Genotypes (G) 794 0.125 ***
Years (Y) 2 0.190 ***
G x Y 794 0.006 NS
Blocks (B) 1 0.337 **
G x B 794 0.030 ***
[129]Open in a new tab
SOV, source of variation; DF, degree of freedom; MS, mean square; ***
significance, p > 0.001; **, p > 0.01; NS, p < 0.05.
The studied traits showed a significant positive correlation among
them. However, the correlation between PR and cLr was nonsignificant
([130]Table 3). The highest correlation was between PR and M. On the
other hand, M and PR showed relatively less correlation with cLr (0.266
and 0.023) and SM (0.293 and 0.608), respectively. A similar
correlation pattern was observed among traits in indica subpopulation,
while in japonica subpopulation, the cLr was observed to be highly
correlated with PR and M (0.799).
TABLE 3.
Pearson correlation coefficient (r) calculated between 10 pairs of
lodging-related traits in structured rice population of 795 accessions.
PR SM M cLr
PR 1 0.339 0.90 0.799
SM 0.293 1 0.339 0.68
M 0.894 0.608 1 0.799
cLr 0.023^ns 0.740 0.266 1
[131]Open in a new tab
Ns, non-significant at p < 0.05; other all values are significant at p
< 0.05; below, diagonal are for whole population, while above diagonal
are for japonica subpopulation.
3.2 Population structure and relative kinship
A genetically, morphologically, and geographically diverse sample of a
large population was used in this study. The whole population was
genotyped and resulted in 3.3 million SNPs after cleaning and filter
screening. These SNPs covered almost 373,245,519 bp of the whole rice
genome with an average SNP density of 12.3 SNPs per kb. The 3.3 million
SNPs from the whole population of 795 genotypes were evaluated to
estimate the kinship matrix and principal components (PCs) and resulted
in a strong population structure. The PC1 revealed the maximum
variation, which divided the whole population into two distinct groups.
The same results were supported by a phylogenetic tree and a kinship
matrix ([132]Supplementary Figure S1). The group-I was composed of 289
japonica accessions, while the group-II framed the 506 indica
genotypes, hence denoted as japonica and indica subpopulations,
respectively. The other structural features were consistent with those
previously reported ([133]Zhang et al., 2009; [134]Zhang et al., 2011).
3.3 Association mapping and QTL analysis
The GWAS was conducted to detect the genomic regions associated with
lodging resistance in rice. Depending upon population structure, GLM
and cMLM regression analysis were adopted. In the first step, the
association of SNP markers to lodging resistance traits was estimated
by GLM and cMLM models for japonica, indica subpopulations, and the
whole population separately ([135]Figure 2, [136]Supplementary Figures
S2-S4). Then, the 1,000-time permutation with randomized phenotypes was
performed by MLM to validate the GWAS results and to remove the false
positives. The threshold–log10 (p) ≥ 4 was set for calling the
significant associations. It is suggested that more association signals
in a common region are caused by linkage disequilibrium (LD) decay,
while the lower LD decay rate can reduce the power of association
analysis for complex traits. The 0.28 and 0.25 LD decay rate with the
167 kb and 123 kb genome-wide R ^2 reduction for japonica and indica
subpopulations has already been reported ([137]Huang et al., 2010;
[138]Huang et al., 2012). Hence, in this study, the genomic region
containing a minimum of three significant SNPs within a 170 kb range
was defined as a QTL. By these criteria, a total of 55, 113, 140, and
67 QTLs were observed for SM, PR, M, and cLr, respectively, which were
defined by 2210, 766, 3548, and 383 significant SNPs ([139]Table 4, and
[140]Supplementary Table S2). Among these QTLs, 15, 80, 70, and 11 QTLs
were commonly observed in various populations and statistical models
([141]Table 4).
FIGURE 2.
[142]FIGURE 2
[143]Open in a new tab
Manhattan plot for bending moment at breaking M, indicating the
regression values for SNP markers in japonica, indica subpopulations,
whole population, and whole population after 1000-times permutation
with the (A) general linear model (GLM) and (B) compressed mixed linear
model (cMLM), the dotted lines showed the threshold level at -log(P) ≥
4.
TABLE 4.
Summary of QTLs identified for lodging resistance traits in three
populations by different regression models.
Traits Whole population Indica sub‐population Japonica sub‐population
Common QTLs Total QTLs
GLM MLM GLM MLM GLM MLM
SM 19 32 8 10 4 6 15 55
PR 13 19 3 9 73 82 80 113
M 18 45 13 18 64 63 70 140
cLr 17 21 8 23 6 7 11 67
Total 67 117 32 60 147 158 176 375
[144]Open in a new tab
3.4 Biological correlation and pleiotropic QTLs
For the studied traits, some QTL regions were commonly detected for
multiple phenotypes, which indicated the biological correlation among
lodging resistance parameters. Among all of the identified QTLs, 33 QTL
regions were found to be overlapped for more than one trait ([145]Table
5). The QTL regions may be considered the potential pleiotropic loci
for lodging resistance–related traits. The PR and M were observed as
the most associated traits revealed by Pearson coefficient of
correlation ([146]Table 3, [147]Table 5) and showed the maximum
biological correlation with 20 overlapping QTLs. The one-third
([148]Chigira et al., 2020) of total 33 QTLs were pleiotropically
governing at least three lodging resistance traits. The maximum number
of pleiotropic QTLs was found on chromosome 1, followed by chromosomes
7 and 11 with five and four QTLs, respectively. There was one genomic
region of 316 kb on chromosome 7 (qSM7-1, qPR7-1, qM7-2, qcLr7-1), that
was significantly associated with all four lodging related traits.
TABLE 5.
List of QTLs with pleiotropic effect, situating at common QTL regions
in the rice genome for lodging resistance traits.
Pleiotropic QTL Chromosome Start position End position Traits QTLs
count Size of QTL (bp) Candidate ORFs
SM PR M cLr
1 1 3,538,854 3,683,583 qSM1-1 qPR1-3 qM1-1 3 144,729 19
2 1 5,566,150 5,744,781 qPR1-4 qM1-3 qcLr1-1 3 178,631 19
3 1 5,849,743 5,971,020 qSM1-3 qM1-4 2 121,277 14
4 1 6,271,684 6,376,084 qSM1-4 qPR1-5 2 104,400 17
5 1 6,581,990 7,070,411 qPR1-6 qM1-5 2 488,421 70
6 1 14,353,502 14,890,707 qPR1-7 qM1-7 2 537,205 73
7 1 24,957,693 27,134,674 qSM1-8 qM1-14, qM1-15 qcLr1-5 4 2,176,981 317
8 1 31,504,595 31,969,329 qPR1-13 qM1-19 qcLr1-6 3 464,734 74
9 1 32,366,242 33,807,002 qPR1-14 qM1-21 qcLr1-7 3 1,440,760 62
10 1 40,397,326 42,294,223 qSM1-12 qPR1-15 qM1-26 3 1,896,897 291
11 2 15,454,980 17,571,479 qPR2-5 qM2-4 2 2,116,499 311
12 2 24,657,034 24,799,865 qPR2-10 qcLr2-5 2 142,831 46
13 3 14,879,518 15,135,921 qSM3-2 qPR3-2 qM3-6 3 256,403 44
14 3 20,883,680 24,068,489 qPR3-10 qM3-13 2 3,184,809 479
15 3 27,631,508 27,677,782 qPR3-11 qcLr3-6 2 46,274 6
16 4 33,314,165 33,457,154 qPR4-12 qM4-12 2 142,989 24
17 5 5,161,253 5,348,092 qSM5-2 qcLr5-2 2 186,839 26
18 6 29,687,352 29,709,960 qSM6-2 qM6-6 qcLr6-3 3 22,608 3
19 7 5,261,202 5,577,829 qSM7-1 qPR7-1 qM7-2 qcLr7-1 4 316,627 52
20 7 7,234,758 7,608,752 qSM7-2 qM7-3 2 373,994 56
21 7 9,435,313 9,925,095 qSM7-3 qM7-4 2 489,782 64
22 7 10,111,738 11,910,229 qSM7-4 qM7-5 2 1,798,491 280
23 7 16,028,149 16,517,487 qPR7-3 qM7-7 qcLr7-3 3 489,338 74
24 8 18,109,039 18,657,995 qPR8-5 qM8-5 2 548,956 74
25 8 18,278,186 20,416,330 qSM8-1 qPR8-5 qM8-6 3 2,138,144 249
26 9 16,303,469 16,907,874 qPR9-4 qM9-8 2 604,405 96
27 10 871,334 1,056,808 qSM10-1 qM10-1 2 185,474 29
28 10 2,641,123 2,829,097 qSM10-2 qM10-2 2 187,974 30
29 11 7,873,690 8,601,084 qSM11-3, qSM11-4 qM11-2 2 727,394 110
30 11 15,327,195 15,530,521 qPR11-8 qM11-3 2 203,326 22
31 11 16,648,356 22,483,866 qSM11-5, qSM11-6 qPR11-9, qPR11-10 qM11-5
qcLr11-3, qcLr11-4, qcLr11-5, qcLr11-6 8 5,835,510 823
32 11 21,216,213 23,048,708 qSM11-8 qPR11-11 qM11-8 3 1,832,495 83
33 12 10,232,335 10,956,543 qSM12-2 qM12-1 qcLr12-1 3 724,208 94
[149]Open in a new tab
3.5 Candidate gene identification
To evaluate the available candidate loci in identified QTLs, the SNPs
in the physical range of annotated loci or their promoters were
observed. A total of 4,031 annotated loci in a total of 30.1 Mb of the
genomic region ([150]Table 5) were revealed for pleiotropic loci in
promoter and coding sequence (CDS) regions ([151]Supplementary Table
S3). The commonly identified QTL on chromosome 7 possesses 59 candidate
genes. As per functional annotations, three out of a total of 59 genes,
viz., OsFBA2, OsFBX222, and LTPL84-protease inhibitor may be considered
the candidate loci for LR-related traits ([152]Supplementary Table S3).
3.6 Functional annotations of candidate genes
All of the candidate loci for pleiotropic QTLs were subjected to
finding the functional annotations. Among the various biological
process-related gene ontological (GO) terms, ‘Transmembrane transporter
activity’ and ‘Localization’ were observed as the top enriched GO
terms. Among the cellular components, ‘membrane’, and among molecular
functions, ‘Biological regulation’ were the top enriched GO terms
([153]Figures 3A,C). Similarly, among the KEGG pathways, ‘Ion binding,’
‘catalytic activity,’ ‘membrane’, ‘localization,’ and ‘transportation’
related pathways were significantly enriched ([154]Figure 3B).
FIGURE 3.
[155]FIGURE 3
[156]Open in a new tab
Top enriched gene ontological (GO) terms (A), top enriched KEGG
pathways (B) annotated for the candidate genes within the genomic
regions of 33 pleiotropic QTLs detected for four lodging resistance
traits, and Venn diagram (C) indicating the overlapping genes among GO
terms including biological processes (BP), cellular components (CC),
and molecular functions (MF).
3.7 In silico and in vitro (qRT-PCR based) expression validation of candidate
loci
We further evaluated the candidate genes on chromosome 7 for their
expression analysis. One strong and thick stem plant, IRAT109, and a
weak and thin stem plant, YueFu, were used to analyze the expression of
these three genes. The results showed a significant increase in
expression levels of OsFBA2 in strong plants, but the expression levels
of LTPL84-proteae inhibitor and OsFBX222 varied non-significantly.
These results were consistent with the database of in silico expression
analysis conducted on tissue-wise expression data extracted from the
Ricexpro database
([157]https://ricexpro.dna.affrc.go.jp/category-select.php). Moreover,
OsFBA2 has the most obvious difference in expression between the two
germplasm ([158]Figure 4). Hence, this gene may not only be considered
a candidate for future studies on lodging resistance but also validated
the authenticity of results from the current study.
FIGURE 4.
FIGURE 4
[159]Open in a new tab
Gene expression analysis of three candidate genes in a pleiotropic QTL
on chromosome 7 of the rice genome. NS; no significant difference, ** p
< 0.01 based on two-tailed Student’s t-tests.
3.8 Haplotype identification
A total of 463 genotypes and 10 SNPs (two SNPs in CDS and eight SNPs in
promoter regions) were used to find the possible haplotypes of
candidate genes. Consequently, the total gene pool could be classified
with the availability of four haplotypes. Among these four haplotypes,
Hap1 and Hap2 represented the indica sub-population, while Hap3 and
Hap4 were mainly distributed in the japonica sub-population
([160]Figure 5A). Among the four haplotypes, Hap3 was prominent, with
significantly high values of M and SM ([161]Figures 5B–D). Hence, the
SNP marker at 5248026 bp physical position on chromosome 7 was
considered a functional SNP to induce lodging resistance ([162]Figure
5).
FIGURE 5.
[163]FIGURE 5
[164]Open in a new tab
Haplotype identification of OsFBA2. (A) Candidate genomic region and
haplotypes of OsFBA2 from 463 rice cultivars, where Tej; temperate
japonica population, Trj; tropical japonica population, Ind; indica
population. (B–E) Phenotypic variation comparison among various
haplotypes for bending moment at breaking (B,D), section modules (C,E).
NS; no significant difference, * p < 0.05, ** p < 0.01 based on
two-tailed Student’s t-tests.
4 Discussion
Lodging is one of the major abiotic stress factors affecting yield loss
of rice, especially in developing Asian countries ([165]Khush, 1997).
Historically, the direct and fundamental method to improve the breeding
material for complex quantitative traits was morphological
marker-assisted selection, where the frequency of favorable alleles
could be increased within a specific population over multiple cycles of
selection. Since the pre-green revolution, many breeding efforts have
been made to improve lodging resistance with improved yield in rice
([166]Meng et al., 2021). Nonetheless, the lodging resistance in rice
crop showed negative correlation to some of the yield parameters as
large spike for high yield made the genotypes prone to lodging
([167]Hirano et al., 2017), dwarfism limits the total yield potential
([168]Islam et al., 2007; [169]Rana et al., 2021), and high fertilizer
application with dense planting made the genotype vulnerable to lodging
in windy and rainy seasons ([170]Shah et al., 2019). Hence, the
identification of potential QTLs and/or genes associated with lodging
resistance traits and their manipulation in breeding material is the
best way to achieve sustainable improvement. It has been reported that
the culm strength for lodging resistance could be contributed by stem
diameter (SD), stalk bending strength (SBS), and the pushing resistance
in field conditions ([171]Ishimaru et al., 2008; [172]Hu et al., 2013),
as a significant negative correlation of culm lodging rate with these
three traits has been observed. Understanding the genetic basis of
these traits will help to induce lodging resistance in the breeding
gene pool.
In this study, a geographically and genetically diverse population of
795 rice genotypes was used to uncover the genetic architecture of culm
lodging resistance in rice. The genotype heterogeneity within a
genotype was removed, and purity was confirmed for 3 years by the SSD
approach, and then phenotypes were investigated for 3 years. The
wide-range, optimum normal frequency distribution and significant
phenotypic diversity in the population for studied traits showed the
polygenic control of traits. A significant genetic potential of
genotypes for lodging resistance was observed by high heritability
values. The high-resolution genotyping with 13X sequencing and 12.3
SNPs per kb markers density enhanced the power of genome-wide
association analysis by improving mapping resolution and imputation
efficiency ([173]Huang et al., 2010). Hence the deep sequence data from
a significant number of landraces acquired the maximum genetic
diversity, which was a key target to reveal in GWAS ([174]Rashid et
al., 2022).
A genome-wide association study is an effective tool to dissect the
genetic basis of complex traits like lodging resistance and its
components. However, it is always pretentious due to the occurrence of
false positives and true negatives ([175]Atwell et al., 2010). In
different studies on human and plant genomes, two factors, as
estimation of population structure and genetic control, were observed
to be suitable to avoid the detection of false-positive signals
([176]Devlin and Roeder, 1999; [177]Marchini et al., 2004; [178]Yu et
al., 2006). In the current study, the whole population was divided into
indica and japonica subpopulations on the basis of population
structure, Q-matrix, kinship-matrix, and PCs. Then, the marker trait
association was separately and recurrently evaluated for the whole
population and indica and japonica subpopulations with multiple
regression modals including GLM, MLM, and cMLM. This approach has
previously been reported to reduce false positives ([179]Yang et al.,
2010; [180]Yang et al., 2011). The detection of a clear population
structure as in indica and japonica subpopulations, the use of multiple
regression models for association estimation, and the 1000-time
permutation test significantly reduced the detection of false
positives. Furthermore, GWAS resolution can be judged by LD estimation
([181]Flint-Garcia et al., 2003). The 167 kb and 123 kb LD decay rate
in japonica and indica subpopulations has been reported in previous
studies ([182]Huang et al., 2010; [183]Huang et al., 2012). Thus, in
this study, the threshold criteria of at least three SNPs in a 170 kb
LD block were followed for defining the QTLs to enhance the power of
GWAS resolution. Eventually, we mapped the 375 genomic regions
significantly associated with lodging resistance traits in rice. The
recurrently identified genomic regions for different traits showed
their significant association with lodging resistance. These regions
may be candidates for further gene cloning and marker-assisted breeding
studies.
Previous genetic research studies on lodging traits used the GWAS
approach and successfully identified the major QTLs explaining 65.7% of
the variance for maximum load to breaking moment and critical mass
([184]Hu et al., 2013). The identification of SCM1, SCM2, SCM3, and
SCM4 governing culm strength endorsed the contribution of SM, M, PR,
and cLr to lodging resistance ([185]Ookawa et al., 2010). Hence, we
looked for the pleiotropic loci for these traits and found 33
overlapping QTLs for more than one trait. A strong correlation among
traits was also observed from the phenotypic data. To confirm the
analysis, a pleiotropic QTL on chromosome 7 governing all four studied
traits was evaluated for candidate genes. Among the 59 open-reading
frames in this QTL range, three genes, including F-box domain protein
genes OsFBA2, OsFBX222, and a LTPL84-proteae inhibitor genes, were
identified. We further evaluate the qRT-PCR based in vitro expression
of candidate genes in thick-culm and thin-culm plants. The results
showed the significantly high expression levels of LTPL84-proteae
inhibitor and OsFBA2 in strong plants, which confirmed the role of
these genes in culm strength and lodging resistance. The F-box domain
protein genes like SCM2 have already been known to contribute to the
culm strength in rice. We further evaluate the qRT-PCR based in vitro
expression of candidate genes in strong and weak genotypes. A
significant reduction in the expression level in weak genotypes
confirmed the role of these genes in culm strength and lodging
resistance. Hence, it could be proposed as a candidate gene for further
investigation. The identification of candidate genes and their
expression-based validation proved the accuracy and reliability of
analytical methodology in current research. The isolation and
manipulation of candidate genes in pleiotropic loci may be beneficial
to enhance the culm strength, leading to lodging resistance in rice.
Furthermore, using the same approach for other pleiotropic and
non-pleiotropic loci, multiple candidate genes could be extracted
([186]Supplementary Table S3).
The GO analysis of candidate genes from all pleiotropic QTLs supported
the results for lodging resistance as the gene ontological terms for
biological processes ‘transmembrane transport activity’ and
‘localization’ in cellular component ‘membrane’ were the most enriched,
while ‘biological regulation’ was the most enriched molecular function.
Similar research studies as strigolactone signaling ([187]Sang et al.,
2014; [188]Yano et al., 2015) and the lignin accumulation ([189]Yoon et
al., 2015) to enhance the culm strength have been reported. These
results were reinforced by the identification of the ‘Biosynthesis of
Cofactors,’ ‘Sphingolipid Metabolism,’ and ‘Pyrimidine Metabolism’
pathways, which had their role in cellular localization and required
for cell proliferation ([190]Siddiqui and Ceppi, 2020).
The pyramiding of the candidate genes with other genes of agronomic
importance may enhance the overall yield and culm strength in rice.
Nonetheless, due to the epistatic nature of some genes, it may be
necessary to find the most favorable alleles of candidate genes. Hence,
we performed the haplotype analysis to reveal the functional allele for
the candidate gene OsFBA2 and observed that Hap3 was the most promising
haplotype in japonica sub population. These genes have not been
reported to contribute to culm strength. Hence, we speculate that
OsFBA2 is the most likely candidate gene for pleiotropic QTL-19 on
chromosome 7 ([191]Table 5). Further studies such as gene cloning and
functional analysis for these two and the other identified candidate
genes in this study will be helpful in understanding their genetic
mechanism to enhance culm strength.
Data availability statement
The original contributions presented in the study are publicly
available. This data can be found here: PRJEB6180.
Author Contributions
MR, YoZ, YaZ, HZ, ZL, and YL conceived and designed the experiments,
performed the experiments, analyzed the data, prepared figures and/or
tables, authored or reviewed drafts of the manuscript, and approved the
final draft. FA, RA, YP, HM-DA, DL, and ZZ validated the research,
review-edited the manuscript. HZ and ZL supervised the research and
provided the resources.
Funding
This work was supported by grants from the Ministry of Science and
Technology of China (No. 2021YFD1200502), the National Natural Science
Foundation of China (No. 32072036), the Open Project of Guangxi Key
Laboratory of Rice Genetics and Breeding
(2022-36-Z01-KF032018-15-Z06-KF10). The funders had no role in the
study design, data collection, and analysis, the decision to publish,
or the preparation of the manuscript.
Conflict of interest
The authors declare that the research was conducted in the absence of
any commercial or financial relationships that could be construed as a
potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors
and do not necessarily represent those of their affiliated
organizations, or those of the publisher, the editors, and the
reviewers. Any product that may be evaluated in this article, or claim
that may be made by its manufacturer, is not guaranteed or endorsed by
the publisher.
Supplementary material
The Supplementary Material for this article can be found online at:
[192]https://www.frontiersin.org/articles/10.3389/fgene.2022.960007/ful
l#supplementary-material
[193]Click here for additional data file.^ (10.1KB, xlsx)
[194]Click here for additional data file.^ (376.8KB, xlsx)
[195]Click here for additional data file.^ (208.8KB, jpg)
[196]Click here for additional data file.^ (249KB, jpg)
[197]Click here for additional data file.^ (9.8KB, xlsx)
[198]Click here for additional data file.^ (246.7KB, jpg)
[199]Click here for additional data file.^ (97.6KB, pdf)
References