A novel discovery of a long terminal repeat retrotransposon-induced hybrid weakness in rice.

Hybrid weakness is a post-zygotic hybridization barrier frequently observed in plants, including rice. In this study, we describe the genomic variation among three temperate japonica rice (Oryza sativa ssp. japonica) varieties 'Aranghyangchalbyeo' ('CH7'), 'Sanghaehyangheolua' ('CH8') and 'Shinseonchalbyeo' ('CH9'), carrying different hybrid weakness genes. The reciprocal progeny obtained from crossing any two varieties displayed characteristic hybrid weakness traits. We mapped and cloned a new locus, Hwc3 (hybrid weakness 3), on chromosome 4. Sequence analysis identified that a long terminal repeat (LTR) retrotransposon was inserted into the promoter region of the Hwc3 gene in 'CH7'. A 4-kb DNA fragment from 'CH7' containing the Hwc3 gene with the inserted LTR retrotransposon was able to induce hybrid weakness in hybrids with 'CH8' plants carrying the Hwc1 gene by genetic complementation. We investigated the differential gene expression profile of F1 plants exhibiting hybrid weakness and detected that the genes associated with energy metabolism were significantly down-regulated compared with the parents. Based on our results, we propose that LTR retrotransposons could be a potential cause of hybrid weakness in intrasubspecific hybrids in japonica rice. Understanding the molecular mechanisms underlying intrasubspecific hybrid weakness is important for increasing our knowledge on reproductive isolation and could have significant implications for rice improvement and hybrid breeding.

Introduction

Increased heterozygosity through wide hybridizations is generally favored in crops for superior qualities of traits and wider adaptability. Hybridization between diverse genomes generally generates an F1 generation with greater biomass and higher yield—a phenomena termed hybrid vigor or heterosis. However, in some cases hybrids between different species show defective development in terms of viability, fertility, and other phenotypic traits. Such defective hybrid characteristics are collectively called hybrid incompatibility and can be further subdivided into hybrid sterility, hybrid inviability, hybrid weakness, and hybrid breakdown, depending on the hybrid phenotype and underlying molecular mechanisms (Coyne and Orr, 2004; Bomblies and Weigel, 2007; Rieseberg and Willis, 2007; Koide ; Chen ; Yu ). These hybrid incompatibilities are the most important reproductive barriers for promoting the process of speciation (Johnson, 2008). Hybrid incompatibilities usually arise due to deleterious interactions between genes, chromosomal inversions, gene duplication, and transposons (Lynch and Force, 2000; Masly ; Kirkpatrick and Barton, 2006). Transposons are DNA fragments that have the ability to move from one position in the genome to another, causing mutations, and are divided into two categories (Wessler, 2006; Rebollo ), retrotransposons and DNA transposons. Retrotransposons are DNA fragments capable of copying themselves via an RNA intermediate, with the copy then moving into a new genomic position (Wessler, 2006). Retrotransposons are divided into long terminal repeat (LTR) and non-LTR retrotransposons, based on the presence of an LTR in their structure. Retrotransposons have been suggested to be involved in hybrid seed failure in interspecific Arabidopsis hybrids (Josefsson ). However, unlike the well-characterized roles of transposons in inducing hybrid incompatibility in the model organism Drosophila (Masly ; Chambeyron ), the involvement of transposons in causing hybrid incompatibility in plants is not yet well established.

Hybrid weakness, a form of hybrid incompatibility that appears at the post-embryonic stage during plant development, is frequently observed throughout plant taxa (Bomblies and Weigel, 2007). Hybrid weakness manifests itself through characteristic dwarfing of the F1 plants, chlorotic phenotype, stunted growth, necrotic tissues, defective root development, and partial or complete sterility (Bomblies and Weigel, 2007; Hannah ; Hatano ). Expression of hybrid weakness can lead to significant decreases in yield and even lethality and is thus a highly undesirable agronomic trait. Hybrid weakness in Asian rice (Oryza sativa L.) was first described in the 1950s (Oka, 1957). Subsequently, many genetic sources of hybrid weakness have been reported in rice crosses (Amemiya and Akemine, 1963; Chu and Oka, 1972; Okuno and Fukuoka, 1999; Ichitani , 2011; Fu ; Chen , 2014). Hybrid weakness in monocots, especially in rice, is frequently reported (Koide ). However, the detailed mechanisms of hybrid weakness in rice are not yet fully understood. Few genes have been narrowed down to sufficiently small genomic intervals (Ichitani , 2011; Kuboyama ). Two loci, Hwi1 and Hwi2, consisting of three genes that cause hybrid weakness in an interspecific F1 hybrid between O. sativa and a wild rice Oryza rufipogon have been identified and cloned (Chen , 2014). Hwi1 consists of two genes encoding proteins with leucine-rich repeats that are receptor-like kinases, while Hwi2 encodes a secreted putative subtilisin-like protease. The F1 hybrids appeared to have greater resistance response at higher temperatures (Chen ).

In addition to interspecific hybrid weakness, intrasubspecific hybrids may also show some degree of hybrid weakness, which is phenotypically similar to that observed in interspecific hybrids (Zhang, 2012; Fu ). Early studies suggested that Hwc1 and Hwc2 were two dominant complementary genes that, when occurring together in a hybrid background, conferred hybrid weakness (Amemiya and Akemine, 1963; Ichitani ). The Hwc2 gene is found to be prevalent among temperate japonica rice cultivars but not among tropical japonica or indica varieties (Ichitani ). The Hwc1 gene is rare and is reported to be present only in the Jamaica cultivar (Amemiya and Akemine, 1963; Ichitani ). Kuboyama performed high-resolution mapping of Hwc2, narrowing the area of interest down to 19 kb, and identified five cDNAs in this region on the long arm of chromosome 4. Ichitani performed linkage analysis and fine mapped Hwc1 to a 60 kb region on the long arm of chromosome 1.

This study endeavored to provide a comprehensive genetic and functional analysis of hybrid weakness in rice. We describe an intrasubspecific hybrid weakness in F1 plants derived from the japonica rice varieties ‘Aranghyangchalbyeo’ (‘CH7’), ‘Sanghaehyangheolua’ (‘CH8’), and ‘Shinseonchalbyeo’ (‘CH9’) carrying different hybrid weakness genes (Wei, 2013). We observed that F1 hybrids produced from reciprocal crossing between any two of these plants exhibited the hybrid weakness phenotype. The plant material studied showed significant genomic variation. In our previous studies, we had reported that ‘CH8’ contained the rare hybrid weakness gene Hwc1 (Zhang, 2012). In the present study, we identified and cloned a new locus, Hwc3, in ‘CH7’ that induced the expression of hybrid weakness in japonica rice hybrids. We suggest that the intrasubspecific hybrid weakness observed in the reciprocal hybrids in our study might be due to the interaction of three genes, Hwc1, Hwc2, and Hwc3. Taken together, our results suggest the possible involvement of retrotransposons in the expression of hybrid weakness in japonica rice hybrids. These studies will not only improve our understanding of reproductive isolation but will also assist in developing new strategies for crop improvement.

Materials and methods

Plant material

Three temperate japonica rice (O. sativa ssp. japonica) varieties ‘Aranghyangchalbyeo’ (‘CH7’, white pericarp, aromatic, and waxy type), ‘Sanghaehyangheolua’ (‘CH8’, purple pericarp, aromatic, and waxy type) and ‘Shinseonchalbyeo’ (‘CH9’, white pericarp and waxy type) carrying different hybrid weakness genes were used in this study. ‘LiyuB’, a wild-type japonica rice, was used as the control plant. Reciprocal crosses among ‘CH7’, ‘CH8’ and ‘CH9’ were performed to generate six F1 hybrid genotypes (‘CH7/CH8’, ‘CH8/CH7’, ‘CH7/CH9’, ‘CH9/CH7’, ‘CH8/CH9’, and ‘CH9/CH8’).

Phenotypic characterization of hybrid weakness

In order to characterize the phenotype of hybrid weakness, quantitative analysis of important traits, namely plant height, root length, panicle length, and panicle number per plant, was conducted using 10 mature replicate plants of each of ‘CH7’, ‘CH8’, ‘CH9’, and their eight F1 hybrids. Plant height of seedlings was recorded at 3-d intervals starting from 3 d after transplantation to 60 d after transplantation. Plant height was measured as the distance from the basal nodes to the uppermost tips of the leaf blades and/or panicles. For root length measurements, 10 mature plants from each line at the 3-leaf stage were removed from the soil and the root length was measured (cm) for a period of 7 d. Panicle length was measured as the average value from the panicle neck to the panicle tip based on an evaluation of three panicles from each of 10 mature plants from each line at the maturity stage. To check the significance of results, the data were analysed by employing Prism 6.0 using one-way-ANOVA.

Preparation of binary constructs and rice transformation

To prepare the constructs for complementation tests, a 2-kb fragment, designated AT70 and consisting of the LTR retrotransposon, and a 4-kb fragment, designated AT71 and consisting of the LTR retrotransposon and the Hwc3 gene promoter and gene, were amplified using primers Hwc3a F/R and Hwc3b F/R (Supplementary Table S1 at JXB online). PCR was performed with an initial denaturation step of 3 min at 95 °C, with two PCR cycles (1 min at 94 °C, 1 min at 55 °C, and 4 min at 72 °C) followed by 33 PCR cycles (1 min at 94 °C, 1 min at 62 °C, and 4 min at 72 °C) and a final extension step at 72 °C for 10 min. The amplified products were separated by electrophoresis on a 0.8% (w/v) agarose gel stained with ethidium bromide. The amplified fragments were cloned into the HPE-203 binary vector to generate binary constructs. These constructs were transformed into Agrobacterium tumefaciens strain EHA105 and were used for transformation of ‘CH8’ embryogenic callus, with transgenic plants being generated.

RNA extraction and RT-PCR analysis

To confirm the presence of the Hwc3 candidate gene and to examine its expression pattern, RT-PCR analysis was performed on samples from the parental lines and F1 plants. Fresh leaf samples at tillering stage were taken from the plants and the total RNA was extracted using the TRNzol reagent (TRNzol, TianGen Biotech Co. Ltd, Beijing, China). The cDNAs were synthesized (RevertAid First Strand cDNA Synthesis Kit, Thermo Fisher Scientific, MA, USA) from 2 mg total RNA, according to the manufacturer’s protocol. The cDNA primers for the genes were designed using Primer-BLAST software (https://www.ncbi.nlm.nih.gov/tools/primer-blast/). The primers for RT-PCR are listed in Supplementary Table S1. RT-PCR was performed with the initial denaturation step of 3 min at 95 °C, followed by 30 RT-PCR cycles (20 s at 94 °C, 20 s at 58 °C, and 20 s at 72 °C) and a final extension step at 72 °C for 7 min. The amplified products were separated on 1% (w/v) agarose gel stained with ethidium bromide. The rice actin gene (ACT1) was used as an internal control to normalize the expression of the tested genes. Gel imaging and documentation was performed using the Gel Doc™ XR molecular imager (Bio-Rad, CA, USA).

Bioinformatics analysis and characterization of the Hwc3 candidate gene

The Hwc3 gene was annotated according to the full-length cDNA sequence on the NCBI database. Sequence similarity was analysed using the BLAST program (http://www.ncbi.nlm.nih.Gov/ BLAST/) provided by the NCBI website. The promoter sequences were analysed through the Plant CARE database (http: bioinformatics. psb.ugent.be/webtools/plant care/html/). Protein homology was performed using BLASTp databases and a phylogenetic tree was constructed with ClustalW software (https://www.ebi.ac.uk/Tools/msa/clustalo/) using the neighbor-joining method.

Whole genome re-sequencing

Total genomic DNA was extracted from a leaf sample following a modified CTAB method (Doyle and Doyle, 1990). DNA degradation and contamination were monitored on 1% agarose gels. DNA purity was checked using the NanoPhotometer® spectrophotometer (Implen, CA, USA). DNA concentration was measured using the Qubit® DNA Assay Kit in a Qubit® 2.0 Fluorometer (Life Technologies, CA, USA). A total amount of 1.5 µg gDNA per sample was used for library preparation. Libraries were generated using the TruSeq Nano DNA HT Sample Preparation Kit (Illumina, CA, USA) following the manufacturer’s recommendations. Briefly, the DNA sample was fragmented by sonication to a size of 350 bp, then the DNA fragments were end-polished using T4 DNA polymerase, A-tailed, and ligated with the full-length adapter (Illumina, CA, USA) for Illumina sequencing with further PCR amplification. Finally, PCR products were purified and libraries were analysed for size distribution using an Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA) and quantified using real-time PCR. The genomic re-sequencing was conducted by CapitalBio Technology (Beijing, China) to generate filtered paired-end 100-bp-long reads using the Illumina HiSeq 2000 platform. The raw data were processed using FastQC software (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). The filtered paired-end reads were mapped to the Oryza sativa ‘Nipponbare’ reference genome using the Burrows–Wheeler alignment (BWA) tool with default parameters (Li and Durbin, 2009). The aligned reads in the SAM (Sequence Alignment/Map) file were sorted using SortSam (Picard tools v1.118, http://broadinstitute.github.io/picard/), and the sorted SAM file was converted to a BAM file (binary version of SAM file) for variant calling using VarScan software (http://varscan.sourceforge.net). Single nucleotide polymorphisms (SNPs) and insertions/deletions (InDels) (<50 bp) were calculated and identified with snpEff-4.1 software (Cingolani ). Re-sequencing data have been submitted to the SRA repository. The SRA accession number is PRJNA488896.

Microarray analysis

Three biological replicates (each consisting of total RNA from an independent plant) for each line extracted at the tillering stage were selected for microarray analysis. Total RNA extraction, cDNA and labelled cRNA synthesis and Affymetrix expression GeneChip hybridization were conducted by CapitalBio Technology. Arrays were scanned on an Affymetrix GeneChip scanner 3000 (Thermo Fisher Scientific). Data were normalized using the Robust Multi-Array Average (RMA) algorithm, which included quantile normalization and background subtraction. Differentially expressed genes were identified using SAM, at the reliability threshold (q-value)≤5% combined with standard-fold difference for gene screening with a fold-change≥2 (up-regulation) or ≤−2.0 (down-regulation). Gene description was recorded from NCBI BLAST (https://www.ncbi.nlm.nih.gov/) and Uniprot (http://www.uniprot.org/). The microarray records have been submitted to the NCBI database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE119431) and assigned GEO accession numbers as GSE119431.

Results

Characterization and inheritance of hybrid weakness in japonica rice hybrids

F1 progeny obtained from all the reciprocal crosses were inferior compared with their parents with respect to all traits (Fig. 1A–C). The hybrid weakness symptoms started to appear 5 d after transplantation and became more obvious at the tillering stage. Morphological analysis showed that hybrid plants were weak, pale yellow (chlorotic), of short stature, and with wilted leaves. The overall growth of the F1 plants was very slow. At the seedling stage, seedling height in all the studied hybrids was significantly less than that of the parents (Fig. 1D). Compared with the parents, root development in the hybrid plants was very severely affected in terms of root number as well as morphology, and was very clearly weaker than that of the corresponding parents in all the studied hybrids (Fig. 1E). The majority of the plants that developed hybrid weakness symptoms died before reaching the adult stage, while those plants that did reach maturity produced few tillers. Important agronomic traits, namely tiller number, panicle number, and panicle length in the F1 plants exhibiting hybrid weakness, were significantly reduced compared with their parent plants (Fig. 2).

Fig. 1.

Morphological and phenotypic characterization of hybrid weakness. (A–C) Morphological evaluation of hybrid weakness in ‘CH7’, ‘CH8’, and the reciprocal hybrid ‘CH7/8’ and ‘CH8/7’ F1 hybrid plants (A), ‘CH7’, ‘CH9’, and the reciprocal hybrid ‘CH7/9’ and ‘CH9/7’ F1 hybrid plants (B), ‘CH8’, ‘CH9’, and the reciprocal hybrid ‘CH8/9’ and ‘CH9/8’ F1 hybrid plants (C) at 30 d after transplantation. (D) Plant height comparison of the parental plants and their hybrid weakness F1 progeny at the seedling and tillering stages. (E) Root length comparison of the parental plants and their F1 hybrids with hybrid weakness at the seedling stage. Data points in (D, E) represent the mean ±SE. (This figure is available in color at JXB online.)

Fig. 2.

Comparison of yield components in the parents and their F1 progeny exhibiting hybrid weakness at the maturity stage. Comparison of tiller number (A), panicle length (B) and panicle number per plant (PNPP) (C) between parents and their F1 exhibiting hybrid weakness. Significant difference determined by the one-way-ANOVA: **P<0.01, ***P<0.001; ns, P>0.05).

Genome re-sequencing and variant discovery

Because the three japonica rice varieties used in this study were obtained from different geographic regions, it might be expected that the genetic differences among them would be considerable. To test this, we generated a catalog of sequence variants using whole-genome Illumina next-generation sequencing data. The raw data were filtered and mapped to the O. sativa ‘Nipponbare’ reference genome (Supplementary Table S2). In some parallel experiments, the three hybrid weakness plants ‘CH7’, ‘CH8’, and ‘CH9’ were crossed with japonica wild-type rice ‘LiyuB’. The results (Zhang, 2012) showed no hybrid weakness in the F1 plants, which indicated that ‘LiyuB’ does not possess a hybrid weakness gene. Therefore, ‘LiyuB’ was used as a positive control plant. Sequencing results showed that ‘CH7’, ‘CH8’, ‘CH9’, and ‘LiyuB’ produced 123792816, 142658002, 124915310, and 131997874 million reads, respectively, with a mean depth coverage of 30×. Furthermore, in ‘CH7’, ‘CH8’, ‘CH9’, and ‘LiyuB’, 93.82%, 94.06%, 93.88%, and 93.98%, respectively, of the total bases produced high-quality scores (Q30). The genetic variation among the three japonica rice parents, ‘CH7’, ‘CH8’, and ‘CH9’, was individually identified using Samtools. We identified a total of 574 551 polymorphic sites including 487 964 SNPs and 86 587 InDels in ‘CH7’, a total of 1 026 428 polymorphic sites including 874 667 SNPs and 151 761 InDels in ‘CH8’, and a total of 792 465 polymorphic sites including 685 118 SNPs and 10 747 InDels in ‘CH9’. These findings revealed a remarkably high density of molecular polymorphisms among each of the three studied lines. According to their position with respect to genes, SNPs and InDels were divided into different types such as upstream, downstream, intergenic, intron, and exon. The results also indicated that the distribution of SNPs and InDels was mainly in the downstream, exon, intron, intergenic, splice_site_acceptor, splice_site_donor, splice_site_region, transcript, upstream, UTR_3_prime and UTR_5_prime locations (Supplementary Tables S3, S4), and their proportions were very similar in all three parental materials (Supplementary Fig. S1A, B). These regions of the genes are mainly involved in the regulation of gene expression and we speculate that these changes might be involved in changes in gene expression profiles.

Physical linkage map at the Hwc3 locus and prediction of the Hwc3 candidate gene

In this study, a physical linkage map was developed for the locus, based on the previous mapping and sequencing results (Yang, 2005; Wei, 2013). A sequence-tagged site (STS) marker C11112 and a cleaved amplified polymorphic sequence (CAPS) marker C1016 were used for fine mapping around the Hwc3 locus. The Hwc3 locus was located between the two DNA markers, RM3687 and RM5473 (Fig. 3A). The target region was narrowed down to 14 kb, and four candidate genes were identified between the markers KGC4M10 and KGC4M20, approximately 53 kb apart from Hwc2, as previously reported by Kuboyama (Fig. 3A). In order to predict the candidate genes for Hwc3, expression analysis of the genes was studied by RT-PCR. The results indicated that only gene Os.89494 was not expressed in ‘CH7’, ‘CH8’, and ‘CH9’ (Fig. 3B). For further molecular studies of hybrid weakness, we confined our research to ‘CH7’, ‘CH8’, and their F1 hybrid. Sequence analysis of this mapped region in ‘CH7’ and ‘CH8’ was performed. The results revealed that the mapped region showed significant variation between ‘CH7’ and ‘CH8’. The sequence analysis of ‘CH7’ at the Hwc3 locus showed that an approximately 1775 bp region was deleted in the ‘CH8’ genome (Supplementary Fig. S2). A 1740 bp retrotransposon flanked at its ends by a 103 bp LTR was inserted into the promoter part of gene Os.89494 in ‘CH7’ (Fig. 3D), while this insertion was not detected in ‘CH8’. A BLAST search showed that this sequence was found only in the temperate japonica rice varieties and was mostly absent from tropical japonica and indica rice varieties. Since transposon insertions are known to affect transcriptional regulation, the insertion of a LTR retrotransposon may affect the expression of adjacent genes. In order to verify the role of this gene (Os.89494) and the detected polymorphism at the Hwc3 locus in inducing hybrid weakness, we analysed the gene expression pattern of the mapped genes at the Hwc3 locus. Results showed that the two genes, Os.69469 and Os.51535, showed expression in ‘CH7’ and ‘CH8’, as well as in the F1 hybrid (Fig. 3C). Strikingly, the gene Os.89494 was expressed in the F1 hybrid, which showed hybrid weakness, but not in either of the parents (Fig. 3C). Since the parental line that did not show gene expression grew and developed normally, our results suggested that hybrid weakness in the F1 hybrid might be induced by the expression of this gene. Hence, we speculated that Os.89494 could be the candidate gene for the Hwc3 locus.

Fig. 3.

Genetic characterization of hybrid weakness. (A) Physical mapping of Hwc3 on rice chromosome 4. Approximate location of the RM3687 and RM5473 DNA markers used for recombinants around the Hwc2 locus. The square denotes the Hwc2 locus reported by Kuboyama . The circle represents the candidate Hwc3 gene Os.89494. (B) RT-PCR analysis of the candidate Hwc3 gene, where 7 represents the cDNA of ‘CH7’, 8 represents the cDNA of ‘CH8’, 7/8 represents the cDNA of ‘CH7/8’ F1 hybrid and ck represents the genomic DNA of ‘CH7’. The Actin gene was used as the internal control to normalize gene expression. (C) Gene expression analysis at the Hwc3 locus, where 7 represents the cDNA of ‘CH7’, 8 represents the cDNA of ‘CH8’, 7/8 represents the cDNA of ‘CH7/8’ F1 hybrid and ck represents genomic DNA of ‘CH7’. (D) Genomic organization at Hwc3 locus in ‘CH7’ and ‘CH8’. (This figure is available in color at JXB online.)

Complementation test

In order to verify that the gene at the Hwc3 locus is the hybrid weakness gene, we performed a complementation test. A 2-kb fragment (designated AT70), which contained the retrotransposon together with the flanking LTR, and a 4-kb fragment (AT71), which contained the LTR retrotransposon, the promoter part and the gene Os.89494 (Fig. 4A), were amplified from the ‘CH7’ genome. The fragments were cloned into a plant transformation binary vector, HPE203, and the constructs were used to transform ‘CH8’. Overall, 12 and 16 independent T0 transgenic plants were obtained from AT70 and AT71, respectively. We observed that 80% of the transgenic plants carrying AT71 exhibited the hybrid weakness phenotype, while plants carrying AT70 grew normally (Fig. 4B). This indicated that the hybrid weakness in the transgenic plants was associated with only AT71 and not with AT70, suggesting that the gene Os.89494 was able to induce the weakness syndrome in ‘CH8’ only when the LTR retrotransposon was inserted, and that, without the LTR retrotransposon, Os.89494 could not induce hybrid weakness, as depicted in the postulated model of hybrid weakness (Fig. 4C). Previously, Zhang (2012) had detected that the ‘CH8’ genome carried the rare hybrid weakness gene Hwc1. As the parental line carrying only Hwc1 or Hwc3 could grow and develop normally, it is suggested that the hybrid weakness observed in the F1 hybrids derived from ‘CH7/CH8’ resulted from the complementary action of two genes as previously reported (Chen ). We further considered that the insertion of the LTR retrotransposon into the promoter part of the Hwc3 gene was transcribed by the transcriptional regulator, Hwc1, present in the ‘CH8’ genome. Although further experiments are needed to elucidate how genomic DNA from ‘CH8’ (Hwc1) induces the expression of Hwc3, our results preliminarily demonstrate that hybrid weakness in intrasubspecific japonica hybrids may be associated with a two-locus interaction. As the candidate Hwc3 gene is generally not expressed under normal conditions, we speculate that the activation and expression of this gene are responsible for the hybrid weakness phenotype.

Fig. 4.

Genetic complementation and postulated model of hybrid weakness. (A) AT70 and AT71 segments derived from ‘CH7’ used for complementation tests. (B) Independent transformants carrying AT70 and AT71 fragments in ‘CH8’ induced the expression of hybrid weakness. ‘LiyuB’ was used as the wild-type. Each line was represented by five independent transgenic plants. (C) Postulated model for Hwc1–Hwc3 interaction leading to hybrid weakness in ‘CH7/8’ F1 hybrid, where X represents no transcription of the Hwc3 gene and O represents transcription activation of the Hwc3 gene due to the inserted promoter and transcription activator Hwc1 in the ‘CH8’ genome. (This figure is available in color at JXB online.)

The candidate gene of Hwc3 encodes a leucine-rich repeat protein

Hwc3 (Os.89494, AK288361) is 1341 bp. The candidate gene encodes a protein of 446 residues with 10 leucine-rich conserved domains (Supplementary Fig. S3). Apart from the leucine-rich repeats (LRR), no other functional domain was found in the gene. Protein homology analysis search by BLASTp showed that the Hwc3 protein shared high sequence similarity with XA1, the bacterial rice blast-resistance protein from O. sativa (Supplementary Fig. S4). Analysis of the promoter region via Plant CARE software indicated the presence of eight cis-acting regulatory elements related to various physiological and metabolic processes. These include Skn-1, meristem specific activation, TATA box, and CAAT-box (Supplementary Table S5).

Differential gene expression related to hybrid weakness

To obtain a global picture of gene expression changes that occurred in the F1 hybrids exhibiting the characteristic hybrid weakness phenotype, we performed two independent microarray analysis. Three biological replicates were used for each sample. These Affymetrix GeneChip arrays target 44 104 genes in the material studied. Putatively, differentially expressed genes involved in the hybrid weakness phenotype were identified by the following selection criteria: (i) statistically significant differences in gene expression were detected in both datasets; and (ii) the average of the fold-ratio change in expression level was greater than 2.0 or less than −2.0. The number of significantly differentially expressed genes was 443. Comparing the expression of these genes, our dataset indicated that transcript levels of 184 genes exhibited differential expression in the ‘CH7/CH8’ hybrid compared with either parent. Only those genes whose expression level showed a fold-ratio greater than +3.0 or less than −3.0 are listed in Table 1. Supplementary Fig. S5 shows a scatter plot graph and a heat map of expression-distribution patterns of differentially expressed genes after microarray hybridization with labeled cDNA probes obtained from mRNA of the studied materials. With the information from NCBI and Uniport databases, we assigned putative functions to the regulated genes and grouped them into functional categories depending on their biological function. Genes whose products had no known putative function or whose function was not very clear were grouped as unclassified genes. Results indicated that genes with altered expression were involved in various important cellular processes including metabolism, growth and development, signal transduction, and defense-related functions. Moreover, we noticed that among those differentially expressed genes, genes encoding metabolic pathways were over-represented. Contrary to previous studies (Bomblies and Weigel, 2007; Chen ), immune-related genes were not found to be differentially expressed, indicating that the hybrid weakness reported in this study was not related to autoimmunity. Differentially expressed genes involved in photosynthesis showed a high percentage of up-regulated expression, while differentially expressed genes involved in protein modifications and secondary metabolism exhibited a high proportion of down-regulated expression.

Table 1.

List of differentially expressed genes in ‘CH7/8’ hybrid weakness F1 progeny

Annotation	Fold-change	Seq name	Unigene
Down-regulated genes
Metabolism
Thiosulfate sulfurtransferase	−3.0	chr6	Os.10162
Tubulin alpha-4A chain	−3.3	chr12	Os.38240
GDSL esterase/lipase	−3.5	chr5	Os.7727
Tryptophan N-hydroxylase 1	−4.9	chr4	Os.18266
Beta-glucosidase 16	−5.1	chr9	Os.81365
Phosphoribosylanthranilate transferase	−10.5	chr9	Os.26406
Cell growth/ development and other cellular processes
Floral organ regulator 1 (FOR1)	−3.7	chr7	Os.10733
Cellular transportation/transportation
Zinc transporter 5	−5.3	chr5	Os.70563
Phospholipid transfer protein precursor	−3.3	chr11	Os.37890
Cell defense and other cellular responses
RING-H2 finger protein ATL80	−3.0	chr9	Os.87833
Secondary metabolites
(−)-Germacrene D synthase-like	−5.3	chr4	Os.11728
O-Methyltransferase ZRP4-like	−6.1	chr9	Os.53458
Cell homeostasis
Putative glutaredoxin-C12	−6.9	chr11	Os.51468
Unclassified or unknown proteins
Unknown	−3.7	chr6	XM_015787066.1
Unknown	−5.7	chr4	Os.87151
Up-regulated genes
Photosynthesis
Cytochrome P450, family 78, subfamily A, polypeptide 7	19.7	chr10	Os.46789
Chlorophyll a–b binding protein 4	4.0	chr8	Os.7890
Cytosolic glutamine synthetase	7.6	chr3	Os.12728
Cell defense/rescue/responses
Glutathione S-transferase	3.1	chr10	Os.21842
Cellular transportation/translocation
Mitochondrial carnitine/acylcarnitine carrier-like protein	3.1	chr10	Os.62649
Translocon at the inner envelope membrane of chloroplasts 110	3.4	chr10	Os.88107
Transcription factor
F-box protein PP2-B10-like	3.3	chr2	Os.53262
Unclassified or unknown proteins
Unknown	3.8	chr12	XM_015763550.1
Unknown	3.1	chr9	Os.55306

Only those genes which show a fold-ratio greater than +3.0 or less than −3.0 are listed. The annotation identity is from BLAST.

Discussion

Genomic variations and intrasubspecific hybrid weakness in japonica rice

It is well known that the accumulation of genetic differences ultimately leads to reproductive isolation in plants (Widmer ; Baack ). As rice is a self-pollinated species, it is prone to accumulating mutations that can lead to genetic variations (Wright ). In rice, almost all kinds of reproductive isolation mechanisms reported so far in plants have been observed (Nadir ). Hybrid weakness is a common type of reproductive isolation mechanism that is frequently observed in hybrids between wild and cultivated rice or between indica and japonica rice (Ichitani , 2011; Kuboyama ; Chen , 2014). Strong selection during the domestication process, mutation accumulation and adaptation to different environments are the key factors leading to the development of genetic differentiation and subsequently leading to hybrid incompatibilities between rice populations (Nadir ).

Some recent studies, including the current one, have reported that, in addition to interspecific and intersubspecific hybridization in rice, reproductive barriers in the form of hybrid weakness also exist between different lines of the same subspecies, indica or japonica (Zhang, 2012; Fu ). However, the exact understanding of the intrasubspecific or intervarietal hybrid weakness is still not very clear. In this study, the F1 generation clearly exhibited the hybrid weakness phenotype with its characteristic chlorotic phenotype, short stature, aberrant root formation, and markedly decreased tiller number (Figs 1, 2). Consistent with previous studies (Amemiya and Akemine, 1963; Chen ), the hybrid weakness symptoms became apparent at the seedling stage and became more obvious at tillering, indicating that hybrid weakness is expressed during plant development. Our results indicated abundant polymorphic sites (Supplementary Fig. S2; Supplementary Table S4) among the three japonica rice varieties, in the form of SNPs and InDels (Supplementary Fig. S1). We found that there was an uneven distribution of SNPs and InDels on the chromosomes. The DNA polymorphisms were mainly distributed in the intergenic, upstream, and downstream regions of genes. Fewer variants, resulting from InDels and SNPs, could be seen in protein-coding regions, splice junctions, and UTR regions. Non-coding regions that are highly conserved between species contain significantly fewer variants than other non-coding regions (Li ). These regions of the genes are mainly involved in the regulation of gene expression, so these changes might be involved in changes in gene expression. Similar results were previously reported for significant genetic polymorphism in the genomes of indica rice varieties RGD-7S and Taifeng B, whose F1 progeny exhibited the hybrid weakness phenotype (Fu ). Overall, these results indicated that the genetic backgrounds of these japonica rice varieties were very different and that they might have originated and evolved under very different conditions. This genome-level polymorphism might not complement one another when the two genomes are hybridized together, which probably leads to the defective development observed in the F1 progeny. The genome-wide polymorphism results obtained in this study will help further in designing InDel markers for further map-based cloning of hybrid weakness genes.

Intrasubspecific hybrid weakness likely arose via an LTR retrotransposon

The results of the present study suggest that a cis-activation is caused by the insertion of an LTR retrotransposon in the promoter region of the Hwc3 gene in ‘CH7’, which seems to enhance its expression in the F1 hybrid (Figs 3D, 4C). This was further supported by the genetic complementation tests, which showed that a 4-kb sequence of ‘CH7’, including the inserted LTR retrotransposon, the promoter, and the coding sequence of the Hwc3 gene, could induce hybrid weakness in transgenic ‘CH8’ plants (Fig. 4B). The well-known explanation of the genetic basis of the evolution of hybrid incompatibilities comes from the Bateson–Dobzhansky–Muller (BDM) model, which states that hybrid incompatibility is due to negative interactions between at least two genetic loci (Johnson, 2008). However, many cases of hybrid incompatibility have been identified that showed other genetic basis for the evolution of hybrid incompatibilities, including chromosomal inversions, gene duplications, and gene transpositions. LTR retrotransposons are the most abundant class of transposons in plants (Wessler, 2006). Barbara McClintock noted that the activity of transposons can cause restructuring of the genome of hybrids. The involvement of a transposon or retrotransposon in inducing hybrid incompatibility has been previously reported in Arabidopsis and Drosophila (Blumenstiel and Hartl, 2005; Masly ; Josefsson ), although the transposon’s effects on hybrid fitness were not clear (Maheshwari and Barbash, 2011). Our study demonstrates that hybrid weakness in the present case is caused by a promoter on the LTR retrotransposon inducing gene expression at the Hwc3 locus. The candidate Hwc3 gene encoded an LRR protein (Supplementary Fig. S3), a type of protein that is mainly involved in protein–protein interactions. The Hwc3 gene was found not to be expressed under normal conditions and the plants carrying this gene grew normally, which suggests that the gene has no apparent function and is a silenced gene. However, upon hybridization of ‘CH7’ with ‘CH8’, transcriptional activation occurred and the gene started to be expressed in the hybrid background, a finding that suggests that the ‘CH8’ genome contains some transcriptional factor that activates the gene expression in the F1 hybrid (Fig. 4C). Interactions between transcription factors and the promoter region of genes has been found to be an underlying mechanism of heterosis in rice F1 hybrids (Zhang ). Based on our current data, we hypothesize that the transcriptional activation of Hwc3 by a transcriptional regulator, Hwc1, leads to hybrid weakness in japonica rice.

Patterns and mechanisms of gene expression changes during hybrid weakness

To identify the candidate genes that play key roles in hybrid weakness in the intrasubspecific hybridization of japonica rice, we analysed differential gene expression profiling in the ‘CH7/8’ F1 hybrid at tillering by Affymetrix gene chip microarray. We found that the genes involved in the carbohydrate metabolism pathway had a higher proportion of down-regulated genes, while genes related to cellular defense, calcium signaling and secondary metabolism also tended to be down-regulated, as were genes related to protein biosynthesis and peptide transport (Supplementary Fig. S5; Table 1). Plant metabolism is highly coordinated with development as it influences the overall energy status of the plant (Creelman and Mullet, 1997; Ponnu ). Down-regulation of the genes associated with metabolism leads to drastic changes in the nutrient and energy balance of the plant and causes arrested growth and development (Ponnu ). The weakness phenotype of the F1 progeny observed in our study and the down-regulation of energy-metabolism-related genes might suggest that the observed phenotype could be due to the low energy status of the hybrids. Previously, carbohydrate metabolism-related genes were found to be up-regulated in heterotic rice hybrids (Zhang ; Wei ). Based on our data, we suggest that the weakness-associated phenotypic changes may result from quantitative differences in gene expression. The differences in gene expression pattern and gene expression level observed in our study might represent intravarietal expression polymorphisms. To dissect the complex mechanisms underlying heterosis, several researchers have attempted to identify genes leading to heterosis when species from different genetic backgrounds are crossed together (Guo ; Huang ; Zhu ). Changes in gene expression have been documented in hybrids of Arabidopsis (Andorf ; Fujimoto ; Shen ), O. sativa (Wei ; He ; Song ), Zea mays (Guo ; Stupar and Springer, 2006; Swanson-Wagner ; Lai ), and Triticum (Wang ). However, differential gene expression in relation to hybrid weakness had not been explored until this current study. The present study identified differentially expressed genes that represent a wide range of molecular functions. This will help to identify candidate hybrid weakness genes and to study their molecular mechanisms. As LRR-containing proteins are mainly involved in protein–protein interactions that regulate a variety of cellular functions, including cell cycle progression, signal transduction, and metabolic pathways, our results suggest that the expression of Hwc3 might regulate the transcriptional functions of other proteins encoded by the differentially expressed genes in the ‘CH7/CH8’ F1 hybrid.

Supplementary data

Supplementary data are available at JXB online.

Fig. S1. Proportions of InDels and SNPs in different gene regions of the japonica rice varieties ‘CH7’, ‘CH8’, and ‘CH9’.

Fig. S2. Sequence analysis at the Hwc3 locus in ‘Nipponbare’, ‘CH7’, and ‘CH8’.

Fig. S3. The protein sequence of the candidate Hwc3 gene product and its phylogeny.

Fig. S4. Protein homology of the Hwc3 and XA1 genes using the NCBI database.

Fig. S5. The differential gene expression patterns related to hybrid weakness.

Table S1. Primer sequences used in this study.

Table S2. Coverage of the reads mapping to the ‘Nipponbare’ reference genome from re-sequencing of the japonica rice varieties ‘CH7’, ‘CH8’, ‘CH9’, and ‘LiyuB’.

Table S3. Annotation of InDels identified in ‘CH7’, ‘CH8’ and ‘CH9’ and ‘LiyuB’ in comparison with ‘Nipponbare’.

Table S4. Annotation of SNPs identified in ‘CH7’, ‘CH8’ and ‘CH9’ and ‘LiyuB’ in comparison with ‘Nipponbare’.

Table S5. Bioinformatics analysis of the promoter structure of Hwc3.

Click here for additional data file.

Click here for additional data file.