Rice (Oryza sativa L.) is an important food crop in the world, especially in Asia. About 1/3 of the world's population relies on rice as the staple food (Wang et al., 2018; Zhao et al., 2018), among which Oryza sativa ssp. japonica and Oryza sativa ssp. indica are the most widely planted. (Zhang et al., 1992; Wang et al., 2018). With the continuous growth of world population, human's demand for rice yield and quality is also increasing (Seck et al., 2012). Rice yield and quality are affected by various biotic and abiotic factors in the external environment, among which Rice blast caused by Magnaporthe oryzae is one of the most important biotic factors. Under suitable environmental conditions, rice is affected by Magnaporthe oryzae from seedling stage to grain filling stage. In resistance to M. oryzae infection, resistant varieties (Sharma et al., 2012; Li et al., 2016) and fungicides (Skamnioti and Gurr, 2009) are usually used. However, these two methods have certain deficiencies and limitations. On the one hand, rice varieties carrying species-specific resistance genes will gradually lose their resistance in the face of M. oryzae with the rapid evolution and mutation. On the other hand, the use of fungicides also presents potential environmental safety risks (Ma et al., 2019). Therefore, the application of broad-spectrum resistance gene is very important in the prevention and control of M. oryzae. At present, many broad-spectrum resistance genes have been identified in rice (Li et al., 2019), such as Pi1 (Hua et al., 2012), Pi9 (Liu et al., 2002; Qu et al., 2006), Pi21 (Fukuoka et al., 2001), Pi54 (Sharma et al., 2010), Pita (Jia et al., 2000), bsr-d1 (Li et al., 2017), etc. Most of the changes in gene expression occurred at the transcriptional stage when plants were infected by pathogenic bacteria (Jain et al., 2017). Therefore, using transcriptomics to study the infection of M. oryzae to rice can effectively identify potential candidate genes with resistance function (Jain et al., 2017). Near Iso-genic lines (NIL) are one of the ideal genetic materials for studying the resistance mechanism of rice to M. oryzae (Mackill et al., 1992). Usually, NIL carrying a broad-spectrum resistance gene is constructed for study, and the corresponding analysis method is usually also simple pairwise comparison. Few studies used the same parents to construct NIL containing different broad-spectrum resistance genes to compare the resistance differences of different broad-spectrum resistance genes (Jain et al., 2019).

With the rapid development of high-throughput sequencing technology, researchers have begun to use multi-sample sequencing to study systems biology. The traditional difference comparison analysis method can only compare two samples and cannot effectively deal with the multi-sample high-dimensional data. The emergence and development of network analysis algorithm has changed this dilemma. Network analysis can effectively process high-dimensional data and extract biologically significant differential expression profiles. Weighted gene co-expression network analysis (WGCNA) is a widely used network analysis method. WGCNA is an unsupervised learning method, which is based on the expression spectrum data of different samples for clustering and can identify co-expression modules, screen candidate biomarkers and identify potential target sites among different samples (Langfelder et al., 2008). WGCNA has been used in many plant studies. Mao et al. (2009) used 1 094 gene chip data to construct Arabidopsis gene co-expression network, and identified 382 core genes. Greenham et al. (2017) applied WGCNA to the transcriptional data of Brassica rapa and constructed the gene co-expression network at different time points under drought stress, revealing the physiological and transcriptional characteristics of Brassica rapa at the early stage under drought stress. Shang and Gao (2020) used WGCNA to analyze transcriptional differences in Arabidopsis thaliana infected with powdery mildew. Qin et al. (2020) identified 15 co-expression modules of genes closely related to potato root drought resistance based on WGCNA, and gene annotation found that most of these core genes were related to drought stress regulation. Chang et al. (2020) identified the co-expression genes related to endogenous abscisic acid in response to the stress of Sclerospora graminicola by WGCNA in foxtail millet. Ma et al. (2019) used WGCNA to identify co-expression modules of plant height and ear height in maize. Thirunavukkarasu et al. (2013) used WGCNA to analyze 528 gene chip data of maize (Zea Mays L.) under waterlogging stress and identified 7 functional modules. Ju et al. (2019) used WGCNA to identify co-expression modules of genes related to internode elongation of cotton fruiting branches. Tai et al. (2018) used WGCNA to analyze transcriptomic data from different tissues of tea (Camellia sinensis) and identified 35 co-expression modules, of which 20 modules were significantly associated with the biosynthesis of catechins, theanine and caffeine. And identified hub genes that regulate the metabolism of the three substances. In addition to gene expression profile data, WGCNA can also be used for other high-throughput data. Pei et al. (2017) slightly improved WGCNA and successfully applied it to proteomics and metabolomics data, obtaining a biologically meaningful explanation. DiLeo et al. (2011) compared the performance of principal component analysis (PCA), batch learning self-organizing maps (BL-SOM), and weighted gene co-expression network analysis (WGCNA) on the data of tomato fruit metabolism. Compared with the previous two methods, WGCNA also provided information about co-expression of metabolites while clustering metabolites based on expression quantity, and successfully identified three co-expression modules related to tomato fruit ripening.

Based on RNA-seq data GSE117030 (Jain et al. 2019) from GEO database, this study used WGCNA to construct gene co-expression network, and conducted correlation analysis between gene co-expression module and processing mode. In order to provide new ideas and clues for further study on the molecular mechanism of rice resistance to M. oryzae infection, NIL carrying different broad-spectrum resistance genes were excavated and co-expression network was constructed.

1 Results and Analysis

1.1 Construction of a weighted gene co-expression network

After the soft threshold β=16 is calculated, the similarity matrix is transformed into adjacency matrix, and then the adjacency matrix is transformed into topological overlap matrix (TOM). In order to eliminate the errors caused by background noise and false correlation, the formula dissTom=1-Tom is used to transform the topological matrix into a heterogeneous matrix. Finally, the clustering generated by function hcluster() is cut by dynamic cutting method. Genes with similar expression patterns are clustered on the same branch, and each branch represents a co-expression module. Different colors represent different modules, among which gray modules represent genes that cannot be assigned to any module, so gray modules are not counted as identified co-expression modules.  Differential genes were clustered by correlation according to their expression levels, and genes with higher correlation were assigned to the same module (Figure 1). Turquoise module gene was the most abundant with 1 779 genes. The number of white module genes was the least, with 104 genes (Table 1).

1.2 Specific module identification

After the different treatments were regarded as quality traits and converted into continuous traits, association analysis was conducted with the modules and correlation heat map was made (Figure 2). It was found that the correlation between the tan and the treatment of PB1+Pi1 inoculated with M. oryzae Mo-nwi-53 24 h was 0.83 (P=6e-07). The correlation between the midnightblue module and treatment of PB1+Pi54 inoculated with M. oryzae Mo-nwi-53 24 h was 0.97 (P=2e-14). Similar expression characteristics can also be seen from the gene expression heat map (Figure 3). Tan module and midnightblue module were selected as target modules for subsequent analysis.

1.3 GO enrichment analysis of the target module

To further explore the function of genes in the target module, agriGO was used to perform GO enrichment analysis on tan module and midnightblue. The results showed (Figure 4) that GO items significantly enriched by genes in tan module were mainly related to cell components such as GO:0005829: cytoplasmic solution; GO:0009579: thylakoid; GO:0005856: cytoskeleton; GO:0044444: cytoplasmic part; GO:0005737: cytoplasm.

1.4 Mining of hub genes and construction of interaction network

According to the network structure results generated by WGCNA, the associated nodes with connectivity (weight value)>0.80 were screened for the construction of weighted gene co-expression network and hub gene mining. After data were imported into Cytoscape, the network was analyzed, and the nodes with the top 10 Outdegree values were screened out as the hub genes of the weighted gene co-expression network (Figure 5). Through functional query, it was found that some hub genes of the weighted gene co-expression network were related to plant disease resistance. For example, LOC_Os01g08710.1 was a WRKY transcription factor (Table 2).

2 Discussion

Life science studies tend to be more and more large-sample studies. Groen et al. (2020) used 15 635 transcripts to study the adaptive mechanisms of rice under natural selection. Traditional analysis methods can not effectively classify genes in large samples and can not dig out the biological significance hidden in the data. The advent of network analysis has changed that. Network analysis can classify complex data, efficiently explore the overall gene expression rules between different samples, and calculate the interaction patterns between genes in samples. WGCNA presents the global expression pattern of genes in the samples by modularizing the classification of genes. In addition, the function of unknown genes can be predicted by using known genes to provide ideas and clues for subsequent biological experiments. Therefore, WGCNA plays an important role in gene chip and transcriptome data of multi-samples, and also has important applications in other omics data processing (DiLeo et al., 2011; Pei et al., 2017).

In this study, based on RNA-seq data of rice near-isogenic lines carrying different broad-spectrum resistance genes, a total of 23 specific modules were identified and 2 modules highly correlated with traits were selected for further study. Jain et al. (2019) used the transcriptome data to pally compare rice near-isogenic lines carrying different broad-spectrum resistance genes and analyze their signal network, and found that the combination of broad-spectrum resistance genes Pi9 and Pi54 can maximize rice blast resistance. R genes such as Pi9, Pita, Pi21 and Pi54 play important roles in rice blast broad-spectrum resistance. The near-isogenic lines constructed based on these R genes could effectively cope with the infection of rapidly mutating M. oryzae. Jain et al. (2019) compared the transcriptional differences between different NILs and susceptible varieties by comparative transcriptome, and found that rice plants with different near-isogenic lines had different SDEL (Significant differentially expressed loci) in the face of M. oryzae infection. These SDEL regulate the metabolism of rice plants, such as starch synthesis and degradation, fatty acid synthesis, phospholipid hydrolysis, phenylpropane synthesis, jasmonic acid metabolism, etc. They did not identify a gene interaction network for M. oryza infection, only pairwise comparisons were made. On this basis, WGCNA was used to analyze the transcriptional differences of different near-isogenic lines after infection with M. oryzae, and the co-expression network of highly related genes with different near-isogenic lines was constructed. Two specific modules were identified by constructing a weighted gene co-expression network, and these two specific modules were highly correlated with NIL carrying different broad-spectrum resistance genes in response to M. oryzae infection, and the corresponding biological significance was revealed. Using connectivity as an indicator, the potential hub genes in specific modules were revealed. The results of this study provide new ideas and theoretical guidance for further understanding the molecular mechanism of rice broad-spectrum resistance genes.

3 Materials and Methods

3.1 Data acquisition

Log in GEO to download transcriptome expression profile data of GSE117030 (Jain et al., 2019), and extract FPKM values using R (Version 3.6.3) software as gene expression profile data.

3.2 Module division

R-package WGCNA (Langfelder et al., 2008) was used for weighted gene co-expression network analysis.  Before WGCNA analysis, genes need to be filtered. goodSamplesGenes(), a function provided by R-package WGCNA, was used for gene and sample filtering to improve the accuracy of the co-expression network. The soft thresholding power was determined according to the principle of scale-free networks. The soft thresholding power provided by software was usually used for subsequent analysis. The clustering tree was constructed according to the correlation of gene expression levels, and the co-expression patterns were identified by dynamic pruning method. The minimum number of genes in modules was set to 50, and the modules with similar expression patterns were merged according to the similarity of module eigenvalues (0.75).

Correlation analysis was conducted between traits and each module, and module-trait heat map was drawn. The closer the absolute value of correlation between a trait and a module is to 1, the more the gene of the module is related to the trait. Co-expression modules with biological significance can be found through the heat map for subsequent analysis.

3.3 Functional enrichment of module genes

agriGO (V 2.0) (Du et al., 2010) was used for GO function analysis of genes in candidate modules. With FDR≤0.05 as the threshold, GO terms that meet this condition are considered as differential enrichment GO terms.

3.4 Construction and visualization of gene co-expression networks

Cytoscape_3.7.2 (Shannon et al., 2003) was used to process the gene co-expression network output by WGCNA to screen out the hub genes and visualize them. Each node in the network represents a gene, and the edge represents the mutual regulatory relationship between genes. Gene co-expression network can help us accurately screen out possible hub genes and make functional prediction of genes with unknown functions by using genes with known functions.

3.5 Data statistics

All statistical analyses in this study are based on R 3.6.3. R software packages ggplot2 and pheatmap were used to draw part of the diagrams, and R software package VennDiagram (Chen et al., 2011) was used to draw Venn diagram.

Authors’ contributions

LX and MLN are the experimental designer and executor of this study. They have completed data analysis and written the first draft of the paper. GLW, LYX, YM, ZSS, HXH, ZYY participated in experimental design and analysis of experimental results; HHC is the architect and principal of the project, directing experimental design, data analysis, paper writing and revision. All authors read and approved the final manuscript.

Acknowledgement

This study was supported by the National Natural Science Foundation of China (31801792, 31960554), the International Cooperative Research Center for Yunnan Green Food (2019ZG00901) and Yunnan Academician Workstation of the Chinese Academy of Engineering of Yunnan Province (2018IC063).

References

Chang G.R., Li R.J., Zhang Q., Zhang Y.M., Hang Y.H., and Zhang B.J., 2020. Identification of co-expression genes related to endogenous abscisic acid in response to the stress of Sclerospora graminicola by WGCNA in foxtail millet, Zhongguo Nongye Kexue (Scientia Agricultura Sinica), 53(16): 3280-3293.

Chen H., and Boutros P.C., 2011, VennDiagram: a package for the generation of highly-customizable Venn and Euler diagrams in R, BMC bioinformatics., 12(1): 35.
https://doi.org/10.1186/1471-2105-12-35
PMid:21269502 PMCid:PMC3041657

DiLeo M.V., Strahan G.D., den Bakker M., and Hoekenga O.A., 2011, Weighted correlation network analysis (WGCNA) applied to the tomato fruit metabolome, PLoS ONE, 6(10): e26683.
https://doi.org/10.1371/journal.pone.0026683
PMid:22039529 PMCid:PMC3198806

Du Z., Zhou X., Ling Y., Zhang Z., and Su Z., 2010, agriGO: a GO analysis toolkit for the agricultural community, Nucleic Acids Res., 38(): W64-W70.
https://doi.org/10.1093/nar/gkq310
PMid:20435677 PMCid:PMC2896167

Fukuoka S., and Okuno K., 2001, QTL analysis and mapping of pi21, a recessive gene for field resistance to rice blast in Japanese upland rice, Theor. Appl. Genet., 103(2): 185-190.
https://doi.org/10.1007/s001220100611

Greenham K., Guadagno C.R., Gehan M.A., Mockler T.C., Weinig C., Ewers B.E., and Mcclung C.R., 2017, Temporal network analysis identifies early physiological and transcriptomic indicators of mild drought in Brassica rapa, eLife, 6: e29655.
https://doi.org/10.7554/eLife.29655
PMid:28826479 PMCid:PMC5628015

Groen S.C., Ćalić I., Joly-Lopez Z., Platts A.E., Choi J.Y., Natividad M., Dorph K., Mauck W.M., Bracken B., Cabral C.L.U., Kumar A., Torres R.O., Satija R., Vergara G., Henry A., Franks S.J., and Purugganan M.D., 2020, The strength and pattern of natural selection on gene expression in rice, Nature, 578(7796): 572-576.
https://doi.org/10.1038/s41586-020-1997-2
PMid:32051590

Hua L.X., Wu J., Chen C., Wu W., He X., Lin F., Wang L., Ashikawa I., Matsumoto T., and Wang L., 2012, The isolation of Pi1, an allele at the Pik locus which confers broad spectrum resistance to rice blast, Theor. Appl. Genet., 125(5): 1047-1055.
https://doi.org/10.1007/s00122-012-1894-7
PMid:22643901

Jain P., Dubey H., Singh P.K., Solanke A.U., Singh A.K., and Sharma T., 2019, Deciphering signalling network in broad spectrum Near Isogenic Lines of rice resistant to Magnaporthe oryzae, Sci. Rep., 9(1): 1-13.
https://doi.org/10.1038/s41598-019-50990-8
PMid:31729398 PMCid:PMC6858299

Jain P., Singh P.K., Kapoor R., Khanna A., Solanke A.U., Krishnan S.G., Singh A.K., Sharma V., and Sharma T.R., 2017, Understanding host-pathogen interactions with expression profiling of NILs carrying rice-blast resistance Pi9 gene, Front. Plant Sci., 8: 93.
https://doi.org/10.3389/fpls.2017.00093

Jia Y., McAdams S.A., Bryan G.T., Hershey H.P., and Valent B., 2000, Direct interaction of resistance gene and avirulence gene products confers rice blast resistance, EMBO J., 19(15): 4004-4014.
https://doi.org/10.1093/emboj/19.15.4004
PMid:10921881 PMCid:PMC306585

Ju F.Y., Zhang S.P., Liu S.D., Ma H.J., Chen J., Ge C.W., Shen Q., Zhang X.M., Liu R.H., Zhao X.H., Zhang Y.J., and Pang C.Y., 2019, Identification of co-expression modules of genes related to internode elongation of cotton fruiting branches by WGCNA, Mianhua Xuebao (Cotton Science), 31(5): 403-413.

Langfelder P., and Horvath S., 2008, WGCNA: an R package for weighted correlation network analysis, BMC Bioinformatics, 9: 559.
https://doi.org/10.1186/1471-2105-9-559
PMid:19114008 PMCid:PMC2631488

Li W.T., Chern M., Yin J.J., Wang J., and Chen X., 2019, Recent advances in broad-spectrum resistance to the rice blast disease, Curr. Opin. Plant Biol., 50: 114-120.
https://doi.org/10.1016/j.pbi.2019.03.015
PMid:31163394

Li W.T., Zhu Z.W., Chern M., Yin J.J., Yang C., Ran L., Cheng M.P., He M., Wang K., Wang J., Zhou X.G., Zhu X.B., Chen Z.X., Wang J.C, Zhao W., Ma B.T., Qin P., Chen W.L., Wang Y.P., Liu J.L., Wang W.M., Wu X.J., Li P., Wang J.R., Zhu L.H., Li S.G., and Chen X.W., 2017, A Natural Allele of a transcription factor in rice confers broad-spectrum blast resistance, Cell, 170(1): 114-126, e115.
https://doi.org/10.1016/j.cell.2017.06.008
PMid:28666113

Li Z.Y., Xia J., Chen Z., Yu Y., Li Q.F., Zhang Y.C., Zhang J.P., Wang C.Y., Zhu X.Y., Zhang W.X., and Chen Y.Q., 2016, Large-scale rewiring of innate immunity circuitry and microRNA regulation during initial rice blast infection, Sci. Rep., 6(1): 1-10.
https://doi.org/10.1038/srep25493
PMid:27150822 PMCid:PMC4858701

Liu G., Lu G., Zeng L., and Wang G.L., 2002, Two broad-spectrum blast resistance genes, Pi9(t) and Pi2(t), are physically linked on rice chromosome 6, Mol. Genet. Genomics., 267(4): 472-480.
https://doi.org/10.1007/s00438-002-0677-2
PMid:12111554

Ma J., Cao Y.Y., Wang L.F., Li J.J., Wang H., Fan Y.P., and Li H.Y., 2020, Identification of gene co-expression modules of maize plant height and ear height by WGCNA, Zuowu Xuebao (Acta Agronomica Sinica), 46(3): 385-394.
https://doi.org/10.3724/SP.J.1006.2020.93021

Ma S.W., Lin S.K., Wang M.L., Zou Y., Tao H., Liu W., Zhang LN.., Liang K.J., Ai Y.F., and He H.Q., 2019, Differential Expression Proteins Contribute to Race-Specific Resistant Ability in Rice, Oryza sativa L., Plants, 8(2): 29.
https://doi.org/10.3390/plants8020029
PMid:30678057 PMCid:PMC6410114

Mackill D., and Bonman J., 1992, Inheritance of blast resistance in near-isogenic lines of rice, Phytopathology., 82(7): 746-749.
https://doi.org/10.1094/Phyto-82-746

Mao L., Van Hemert J.L., Dash S., and Dickerson J.A., 2009, Arabidopsis gene co-expression network and its functional modules, BMC Bioinformatics, 10(1): 346.
https://doi.org/10.1186/1471-2105-10-346
PMid:19845953 PMCid:PMC2772859

Pei G., Chen L., and Zhang W., 2017, WGCNA application to proteomic and metabolomic data analysis, Methods Enzymol., 585: 135-158.
https://doi.org/10.1016/bs.mie.2016.09.016
PMid:28109426

Qin T.Y., Sun C., Bi Z.Z., Liang W.J., Li P.C., Zang J.L., and Bai J.P., 2020, Identification of drought-related co-expression modules and hub genes in potato roots based on WGCNA, zuowu xuebao (Acta Agronomica Sinica), 46(7): 1033-1051
https://doi.org/10.3724/SP.J.1006.2020.94130

Qu S.H., Liu G.F., Zhou B., Bellizzi M., Zeng L., Dai L., Han B., and Wang G.L., 2006, The broad-spectrum blast resistance gene Pi9 encodes a nucleotide-binding site–leucine-rich repeat protein and is a member of a multigene family in rice, Genetics, 172(3): 1901-1914.
https://doi.org/10.1534/genetics.105.044891
PMid:16387888 PMCid:PMC1456263

Shang H.Y., and Gao Z.Y., 2020, Analysis of differential transcriptional genes in Arabidopsis thaliana after powdery mildew infection by WGCNA, Fenzi Zhiwu Yuzhong (Molecular Plant Breeding), 18(20): 6734-6742.

Shannon P., Markiel A., Ozier O., Baliga N.S., Wang J.T., Ramage D., Amin N., Schwikowski B., and Ideker T., 2003, Cytoscape: a software environment for integrated models of biomolecular interaction networks, Geno. Res., 13(11): 2498-2504.
https://doi.org/10.1101/gr.1239303
PMid:14597658 PMCid:PMC403769

Sharma T., Rai A., Gupta S., and Singh N., 2010, Broad-spectrum blast resistance gene Pi-k^h cloned from rice line tetep designated as Pi54, J. Plant Biochem. Biotechnol., 19(1): 87-89.
https://doi.org/10.1007/BF03323441

Sharma T., Rai A., Gupta S., Vijayan J., Devanna B., and Ray S., 2012, Rice blast management through host-plant resistance: retrospect and prospects, Agricultural Research., 1(1): 37-52.
https://doi.org/10.1007/s40003-011-0003-5

Skamnioti P., and Gurr S.J., 2009, Against the grain: safeguarding rice from rice blast disease, Trends Biotechnol., 27(3): 141-150.
https://doi.org/10.1016/j.tibtech.2008.12.002
PMid:19187990

Tai Y.L., Liu C., Yu S.W., Yang H., Sun J.M., Guo C.X., Huang B., Liu Z.Y., Yuan Y., Xia E.H., Wei C.L., anf Wan X.C, 2018,  Gene co-expression network analysis reveals coordinated regulation of three characteristic secondary biosynthetic pathways in tea plant (Camellia sinensis), BMC Genomics, 19(1): 616.
https://doi.org/10.1186/s12864-018-4999-9
PMid:30111282 PMCid:PMC6094456

Thirunavukkarasu N., Hossain F., Mohan S., Shiriga K., Mittal S., Sharma R., Singh R.K., and Gupta H.S., 2013, Genome-wide expression of transcriptomes and their co-expression pattern in subtropical maize (Zea mays L.) under waterlogging stress, PLoS ONE, 8(8): e70433.
https://doi.org/10.1371/journal.pone.0070433
PMid:23936429 PMCid:PMC3735631

Wang W., Hu B., Yuan D.Y., Liu Y.Q., Che R.H., Hu Y.C., Ou S.J., Liu Y.X., Zhang Z.H., Wang H.R., Hua I., Jiang Z.M., Gao X.K., Qiu Y.H., Meng X.B., Bai Y., Liang Y., Wang Y.Q., Zhang L.H., Li L.G., Sodmergen, Jing H.C., Li J.Y., and Chu C.C., 2018, Expression of the nitrate transporter gene OsNRT1. 1A/OsNPF6. 3 confers high yield and early maturation in rice, The Plant Cell, 30(3): 638-651.
https://doi.org/10.1105/tpc.17.00809
PMid:29475937 PMCid:PMC5894839

Wang W., Mauleon R., Hu Z., Chebotarov D., Tai S., Wu Z., Li M., Zheng T., Fuentes R. R., Zhang F., Mansueto L., Copetti D., Sanciangco M., Palis K.C., Xu J., Sun C., Fu B., Zhang H., Gao Y., Zhao X., Shen F., Cui X., Yu H., Li Z., Chen M., Detras J., Zhou Y., Zhang X., Zhao Y., Kudrna D., Wang C., Li R., Jia B., Lu J., He X., Dong Z., Xu J., Li Y., Wang M., Shi J., Li J., Zhang D., Lee S., Hu W., Poliakov A., Dubchak I., Ulat V.J., Borja F.N., Mendoza J.R., Ali J., Li J., Gao Q., Niu Y., Yue Z., Naredo M.E.B., Talag J., Wang X., Li J., Fang X., Yin Y., Glaszmann J.C., Zhang J., Li J., Hamilton R.S., Wing R. A., Ruan J., Zhang G., Wei C., Alexandrov N., McNally K.L., Li Z., and Leung H., 2018, Genomic variation in 3,010 diverse accessions of Asian cultivated rice, Nature, 557(7703): 43-49.
https://doi.org/10.1038/s41586-018-0063-9
PMid:29695866 PMCid:PMC6784863

Zhang Q., Maroof M.A.S., Lu T.Y., and Shen Z.B., 1992, Genetic diversity and differentiation of indica and japonica rice detected by RFLP analysis, Theor. Appl. Genet., 83(4): 495-499.
https://doi.org/10.1007/BF00226539
PMid:24202597

Zhao Q., Feng Q., Lu H.Y., Li Y., Wang A., Tian Q.L., Zhan Q.L., Lu Y.Q., Zhang L., Huang T., Wang Y.C., Fan D.L., Zhao Y., Wang Z.Q., Zhou C.C., Chen J.Y., Zhu C.R., Li W.J., Weng Q.J., Xu Q., Wang Z.X., Wei X.H., Han B., and Huang X.H., 2018, Pan-genome analysis highlights the extent of genomic variation in cultivated and wild rice, Nat. Genet., 50(2): 278-284.
https://doi.org/10.1038/s41588-018-0041-z
PMid:29335547