To elucidate the genetic architecture underlying phenotypic divergence is essential to the understanding of ecological adaptation and speciation. Two wild rice species (O. rufipogon and O. nivara) are a progenitor-daughter species pair with ecological divergence and provide a unique system for studying ecological adaptation/speciation. Here, we constructed a high-resolved linkage map and conducted a quantitative trait locus (QTL) analysis of 19 phenotypic traits using an F₂ population generated from a cross between the perennial O. rufipogon and annual O. nivara. We identified 113 QTLs associated with interspecific divergence of 16 quantitative traits, with effect sizes ranging from 1.61% to 34.1% in terms of the percentage of variation explained (PVE). The distribution of effect sizes of QTLs followed a negative exponential, suggesting that a few genes of large effect and many genes of small effect were responsible for the phenotypic divergence. We observed 18 clusters of QTLs (QTL hotspots), with each involving multiple adaptive traits, demonstrating the importance of coinheritance of loci/genes in ecological adaptation/speciation. Analysis of effect direction and v-test statistics revealed that interspecific differentiation of most traits was driven by divergent natural selection, supporting the argument that ecological adaptation/speciation would proceed rapidly under coordinated selection on multiple traits.

Fresh leaves of two parental lines and 600 F₂ individuals randomly chosen from the 862 F₂ population were collected and dried by silica-gel. Genomic DNA was extracted using the cetyltrimethylammonium bromide (CTAB) method (Murray & Thompson, 1980). The libraries of mapping parents were constructed following the manufacturer’s recommendations (Illumina) for 500 bp insert size and sequenced by BGIseq500 platform (BGI; Shenzhen, China) with 150 bp paired-end reads. The F₂ individuals were genotyped by a specific-locus amplified fragment-sequencing (SALF-seq) method (Sun et al., 2013). In brief, two restriction enzymes (RsaI and HaeIII) were selected to digest the genomic DNA. The digested fragments (SLAF tags) were ligated to the adapters with T4 DNA ligase. After PCR amplification, purification, sample mixing and electrophoresis with agarose gels, the size of fragments ranging from 264–314 bp were obtained and purified. Subsequently, the products were sequenced using Illumina HiSeq 2500 platform (Illumina; San Diego, U.S.) with 125 paired-end reads according to the manufacturer’s instruction.

The short reads of parental lines and F₂ individuals were filtered by removing low-quality reads with more than 10% of bases missing. Then, short reads were aligned to the reference sequence of Nipponbare genome (IRGSP-1.0) (Kawahara et al., 2013) using BWA (Li & Durbin, 2010) with the MEM algorithm. Furthermore, Samtools (Li et al., 2009) were applied to sort the mapping results and build an index for each BAM file. Variant calling was conducted using the Genome Analysis Toolkit (GATK, version 4.0.2.1) (Van der Auwera et al., 2013). SNPs were filtered with VariantFiltration of GATK “AC < 2 || QD < 2.0 || FS > 60.0 || MQ < 40.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0” as suggested by the manual.