Genetic and tissue-specific RNA-sequencing analysis of self-compatible mutant TSC28 in Brassica rapa L. toward identification of a novel self-incompatibility factor

Self-incompatibility (SI) is a sophisticated system for pollen selectivity to prevent pollination by genetically identical pollen. In Brassica , it is genetically controlled by a single, highly polymorphic S -locus, and the male and female S -determinant factors have been identified as S -locus protein 11 (SP11)/ S -locus cysteine-rich protein (SCR) and S -locus receptor kinase (SRK), respectively. However, the over-all molecular system and identity of factors in the downstream cascade of the SI reaction remain unclear. Previously, we identified a self-compatible B. rapa mutant line, TSC28, which has a disruption in an unidentified novel factor of the SI signaling cascade. Here, in a genetic analysis of TSC28, using an F 2 population from a cross with the reference B. rapa SI line Chiifu-401, the causal gene was mapped to a genetic region of DNA containing markers BrSA64 and ACMP297 in B. rapa chromosome A1. By fine mapping using an F 2 population of 1,034 plants, it was narrowed down to a genetic region between DNA markers ACMP297 and BrgMS4028, with physical length approximately 1.01 Mbp. In this genomic region, 113 genes are known to be located and, among these, we identified 55 genes that were expressed in the papilla cells. These are candidates for the gene responsible for the disruption of SI in TSC28. This list of candidate genes will contribute to the discovery of a novel downstream factor in the SP11–SRK signaling cascade in the Brassica SI system. mapping and sequencing of papilla-expressed genes from the target genomic region. The list of candidate genes will contribute to discovery of the causal gene of SC in TSC28, coding for a novel downstream factor in the SP11-SRK signaling cascade in the Brassica SI system.


INTRODUCTION
To maintain genetic diversity, plant species have estab-lished several pollination systems, such as dioecy and dichogamy (reviewed by de Nettancourt, 2001). Among these, self-incompatibility (SI) is one of the most elegant pollination systems in plant reproduction. From the definition by de Nettancourt (2001), SI is the inability of fertile hermaphrodite plants to produce zygotes after self-pollination. SI in Brassica species is controlled by a single locus, termed the S-locus, with multiple alleles (Bateman, 1955;Nou et al., 1993). From several molecular analyses, three genes at the S-locus were characterized, which are involved in self/non-self pollen recognition and the subsequent downstream reaction (reviewed by Suwabe et al., 2010;Watanabe et al., 2012;Doucet et al., 2016). The female S-determinant, a receptor-type protein kinase (S-receptor kinase; SRK), consisting of an extracellular domain (S-domain, also known as SLG homologous domain), a transmembrane domain and a serine/threonine kinase domain, has been identified and characterized (Stein et al., 1991;Watanabe et al., 1994;Takasaki et al., 2000). SP11/SCR, a small cysteine-rich protein expressed in anther tapetum cells and/or pollen grains, has been identified as the male S-determinant (Schopfer et al., 1999;Suzuki et al., 1999;Takayama et al., 2000;Shiba et al., 2001). SLG (S-locus glycoprotein; Watanabe et al., 1994), encoding a protein secreted onto the papilla cell wall (Kishi-Nishizawa et al., 1990), is required for full manifestation of the self-pollen recognition reaction, in combination with SRK (Takasaki et al., 2000). Physical interaction between SP11 and SRK, in an S-allele-specific manner, induces the SI response to prevent self-pollination ; thus, a signal from the pollen grain is transduced in stigmatic papilla cells, resulting in pollen rejection.
In self-pollination, information from the SP11-SRK interaction transduces to intracellular molecules of the SI reaction that interact with the kinase domain of SRK. Of these molecules, three (MLPK, ARC1 and THL1) have been identified and characterized (Gu et al., 1998;Cabrillac et al., 2001;Murase et al., 2004). From molecular analysis of a self-compatible (SC) mutant variety, yellow sarson, a membrane-anchoring protein kinase, MLPK (M-locus protein kinase), which interacts with the SRK kinase domain (Murase et al., 2004;Kakita et al., 2007), has been identified as the causal gene for selfcompatibility in yellow sarson. ARC1 (ARM repeat containing 1) interacts with the kinase domain of SRK, in a phosphorylation-dependent manner (Gu et al., 1998), and has E3 ubiquitin ligase activity, suggesting proteasome protein degradation in SI signal transduction (Stone et al., 2003;Doucet et al., 2016). ARC1-repressed transgenic plants, produced by an antisense construct, showed an SC phenotype (Stone et al., 1999). MLPK and ARC1 have been shown to be dominant genes in the SI signaling cascade, because their recessive mutant or repressed transgenic plants have an SC phenotype. THL1 (thioredoxin-h-like protein 1) is reported to be a negative regulator of SRK autophosphorylation in vitro (Cabrillac et al., 2001). While more molecules are expected to be involved in the downstream cascade of the SI reaction, further candidates have not yet been identified in Brassica species (reviewed by Watanabe et al., 2008Watanabe et al., , 2012. To identify all components of the SRK downstream pathway, alternative plant material and/or methodology is necessary. Previ-ously, we identified two novel SC mutant lines (TSC4 and TSC28) in B. rapa (Isokawa et al., 2010), from a screening of bulk populations of Japanese traditional and commercial B. rapa vegetables. The putative SC gene of TSC28 was mapped on linkage group A1, while S-locus, MLPK, ARC1 and THL were mapped on A7, A3, A4 and A6, respectively. Thus, the SC gene of TSC28 is independent from the S-locus or known SI-related genes, indicating a novel SC mutant line. In TSC28, the SC phenotype is inherited as a homozygous recessive allele, like MLPK and ARC1 (Stone et al., 1999;Murase et al., 2004). Thus, a comparison of candidate genes on corresponding genomic regions in an SI line and TSC28 is an informative method to identify the causal gene in the downstream cascade of the SI reaction.
For a genome-wide transcriptome analysis of specific tissues by RNA-sequencing (RNA-seq), deep-sequencing technology by next-generation sequencing is effective and informative (Wang et al., 2009). RNA-seq has been performed on various tissues and/or organs of several plant species (Zenoni et al., 2010;Zhang et al., 2010;Tong et al., 2013;Kim et al., 2014). Recently, this methodology was applied in plant reproductive research to understand pollination mechanisms and reproductive organ development and differentiation (Osaka et al., 2013;Zhou et al., 2014;Matsuda et al., 2015;Xu et al., 2016;Zhang et al., 2017;Takeda et al., 2018). Thus, in combination with a genetic approach, RNA-seq data have potential for identification of the causal gene for the SC phenotype in TSC28. Therefore, in this study we characterized TSC28 from the aspect of fine mapping and papilla-expressed genes in the target genomic region. We have previously established the reference RNA-seq data of papilla cells in the SI line Chiifu-401 (Osaka et al., 2013). From the comparison of RNA-seq data in papilla cells between these SI and SC lines, we identified several genes that are potentially involved in self-compatibility of TSC28 in the corresponding quantitative trait locus (QTL) region of the B. rapa chromosome A1.

MATERIALS AND METHODS
Plant materials and observation of pollen tube behavior on the stigma An SC mutant line, TSC28, which was identified by screening Japanese bulk populations of B. rapa vegetable varieties and partially characterized as a novel SC mutant line that is independent from the S-locus or known SI-related genes (Isokawa et al., 2010), was used. As a reference B. rapa SI line, an inbred line, Chiifu-401 (Osaka et al., 2013), was used. F 2 populations (n = 82 in 2012, n = 82 in 2013, n = 1,034 in 2014) were established by selfing of F 1 between TSC28 and Chiifu-401. All plants were grown in a glass house in the 2012 to 2014 season. Pollen tube behavior in the stigma on pollination was observed by the aniline blue staining method, according to Hatakeyama et al. (1998). The degree of incompatible/compatible phenotype was determined by pollen tube behavior in self-pollination, and allocated to the following six categories based on the rate of pollen tube penetration into the stigma: (0) no pollen tube penetrated the stigmatic papilla or grew into a style, (1) 1-5 pollen tubes penetrated into a style, (2) 6-10 pollen tubes penetrated into a style, (3) 11-15 pollen tubes penetrated into a style, (4) 16-20 pollen tubes penetrated into a style, and (5) more than 20 pollen tubes penetrated into a style.
Genetic analysis of the SC genes of TSC28 DNA markers BrgMS (Xu et al., 2010) and ACMP (Ramchiary et al., 2011) were used for fine mapping of SC genes in TSC28. We developed additional simple sequence repeat (SSR) markers (BrSA, SSRA; Supplementary Table S1) in the corresponding genomic region, using the Brassica Database (BRAD; http://brassicadb.org/brad/; Cheng et al., 2011). PCR analysis of DNA markers was conducted according to Isokawa et al. (2010). PCR amplicons were assessed using the Multina system (Shimadzu, Japan).
To construct the genetic map for the SC gene of TSC28, genetic analysis was performed using F 2 populations, according to Suwabe et al. (2006). For segregation analysis of F 2 progenies grown in 2012 and 2013, a chi-squared test was performed. For fine mapping of the SC gene in TSC28, 1,034 F 2 plants grown in 2014 were used. The position of the SC gene on the linkage map was estimated using MapQTL 5 software (van Ooijen, 2004), as described by Isokawa et al. (2010).
Tissue preparation and laser microdissection Tissue fixation, paraffin embedding, cutting sections of paraffin-embedded stigma samples, removing paraffin from samples, and laser microdissection were performed according to Osaka et al. (2013) and Matsuda et al. (2015). Briefly, chemically fixed pistils were embedded in paraffin and 6-to 8-μm-thick sections were cut using a microtome (RV240, Yamato Koki, Japan). After removing paraffin from the tissue sections using Histo-Clear II (National Diagnostics, USA), the Arcturus XT Microdissection System (Life Technologies, USA) was used for laser microdissection.
RNA sample preparation for RNA-seq and sequence data processing Total RNA was extracted from dissected papilla cells using a Pico-Pure RNA isolation kit (Life Technologies). RNA quality and quantity were evaluated using an Agilent 2100 Bioanalyzer (Agilent, USA). Linear amplification of mRNA and synthesis of cDNA were performed according to Osaka et al. (2013) and Matsuda et al. (2015) with slight modification. Briefly, for SOLiD sequencing, a cDNA fragment library was prepared according to the manufacturer's instructions (Life Technologies). Approximately 0.7 to 1.0 pmol of the cDNA fragment library was applied and sequenced using the SOLiD 5500xl system (Life Technologies). Processing of sequence data for assembly was performed according to Osaka et al. (2013) and Matsuda et al. (2015). Read counts per gene locus were calculated from reads mapped to the genome, and expression values were normalized by reads per kilobase of exon per million mapped reads (RPKM).
qRT-PCR Total RNAs isolated from stigma tissues were reverse-transcribed to first-strand cDNA using the Super-Script III First-Strand Synthesis system (Invitrogen, Life Technologies). The quantity and quality of cDNA were confirmed by spectrophotometric analysis. Quantitative real-time PCR (qRT-PCR) was performed using SYBR Green Realtime PCR Master Mix (Toyobo, Japan) by the CFX96 Real-Time Detection System (Bio-Rad, USA). Quantitative validation was calculated using the delta-delta threshold cycle relative quantification method in three replicates. A primer list for qRT-PCR analysis is given in Supplementary Table S2. As a control gene, UBC21 (encoding ubiquitin-conjugating enzyme 21) was used (Osaka et al., 2013). The PCR reaction mixture (20 μl) consisted of 0.2 μM forward and reverse primers, 1 × SYBR Advantage qPCR Premix and approximately 100 ng cDNA. The PCR cycle conditions were 95 °C for 30 sec followed by 40 cycles of 95 °C for 5 sec and 60 °C for 10 sec.

RESULTS AND DISCUSSION
Genetic analysis of the SC trait in line TSC28 Pollination behavior of TSC28 has been examined previously and identified as SC (Isokawa et al., 2010). In that study, the SC trait of TSC28 was shown to be regulated by a single recessive gene which is independent from the S-locus or known SI-related genes. To confirm the SC phenotype with a single recessive gene model in a different genetic background and a different environmental condition, 82 F 2 plants from a cross between Chiifu-401 and TSC28 were investigated in 2012 and 2013. The segregation ratio for pollen tube behavior (SI:SC) of each F 2 plant was 64:18 in 2012; 61:21 in 2013; and in total for both years 125:39. From the chi-squared test for a single gene segregation ratio goodness-of-fit to ratio 3:1, the resulting segregation ratio for total data from the two years was chi-squared = 0.13, P > 0.05, df = 1, confirming that this SC trait in TSC28 is regulated by a single recessive gene, as described by Isokawa et al. (2010). This also indicates that the SC trait in TSC28 is a consistent phenotype, which is not affected by a different genetic background or environmental conditions.
Mapping and genetic analysis of the SC trait DNA markers described in Materials and Methods were used to construct a linkage map, based on the 2012 and 2013 F 2 populations. From the association between pollination behavior phenotype and genotype of DNA markers in each F 2 progeny, a linkage map for the SC gene was constructed (Fig. 1, Table 1). In QTL analysis in each linkage map, the highest logarithm of the odds (LOD) score was observed at BrSA64 and ACMP297 in B. rapa chromosome A1.
Using 1,034 F 2 plants, we performed a precise linkage analysis between genetic markers and SI/SC phenotype by graphical genotyping (Fig. 2). From this, the SI/ SC phenotype and the genetic markers ACMP297 and BrgMS4028 were found to co-segregate. The physical distance between ACMP297 and BrgMS4028 genomic regions was approximately 1.01 Mbp. This region contains 113 genes, according to the BRAD database, indicating that one of these 113 genes is responsible for the SC trait in TSC28 (Fig. 2).
Characterization of the SC trait genetic region in TSC28 To dissect the potential genomic region of the SC trait in TSC28, comparison of expressed genes in papilla cells in SI (Chiifu-401) and SC (TSC28) lines is an informative approach. For Chiifu-401, data on genes expressed in papilla cells have been published (Osaka et al., 2013). Following the same methodology, we established transcriptome data for papilla cells in TSC28 ( Table  2). As shown in Table 2, the trend of the data (mapped read, mapQ ≥ 10, unigenes reads > 0, unigenes reads ≥ 10) for the sequenced tags of TSC28 was similar to that of Chiifu-401. Because B. rapa is a self-incompatible outcrossing species and has genomic diversity among lines and/or cultivars (Brassica rapa Genome Sequencing Project Consortium, 2011; Golicz et al., 2016), the difference in the sequenced tags should reflect the incompatibility systems in Chiifu-401 and TSC28.
To identify the genes expressed in the papilla cells corresponding to the QTL genomic region, the transcriptome data were compared between Chiifu-401 and TSC28. When a tag count was observed only in  Chiifu-401 or TSC28, we designated these as Chiifu-401or TSC28-specific genes, respectively. In contrast, when a tag count was observed both in Chiifu-401 and in TSC28, we designated these as common genes. In the QTL region for the SC trait in TSC28, 113 genes are located and scattered within a 1.01-Mbp genomic region (Fig. 2). Among these, 55 genes were expressed in papilla cells (Table 3) and 43 genes (78%) were commonly expressed in both lines. The remaining 12 genes (22%) were expressed only in either Chiifu-401 or TSC28 (Fig. 3). Because genetic analysis revealed that the SC trait in TSC28 is regulated by a single recessive gene, the nine Chiifu-401-specific genes and 43 common genes, but not the three TSC28-specific genes, are candidates for the gene responsible for the SC trait in TSC28 (Table  2 and Fig. 3). After confirmation of the expressed genes by qPCR analysis (Fig. 4), we evaluated candidate genes for SC in TSC28, based on the gene descriptions from the list in Table 3. In the Chiifu-401-specific genes, one interesting gene, Bra031482, whose description is "Disease resistance protein (CC-NBS-LRR class) family", was found. In Arabidopsis, AT1G63350.1, a gene orthologous to Bra031482, has been characterized as encoding an N-myristoylated protein, which is involved in R-gene function (Boisson et al., 2003). Among the common genes BrSA64 BrgMS4320 13,565,130 14,363,608 14,714,012 16,999,928 18,019,583 20,889,441 21,190,641 21,360,359 21,491,002 21,596,344 23,659,131 23,988,520      ( Table 3), three are involved in plant-pathogen interaction: Bra031456 (disease resistance protein, TIR class), Bra031426 (seven transmembrane MLO family protein) and Bra031408 (lung seven transmembrane receptor family protein). TIR was first identified as the tobacco mosaic virus resistance gene, N, in tobacco (Whitham et al., 1994) and several TIR-type resistance genes have been identified and characterized so far (Meyers et al., 2005). Bra031426 and Br031408 are classified in the seven transmembrane MLO family and are involved in powdery mildew resistance (Büschges et al., 1997). As described by Hodgkin et al. (1988), Brassica self-incompatibility and plant-pathogen interaction have many common features, suggesting that the SI downstream cascade and the plant-pathogen signaling cascade share common molecular mechanisms and factors. Another type of interesting gene, encoding ARM repeat superfamily proteins, was identified in two common genes, Bra031473 and Bra031431. In B. napus, ARC1 (ARM repeat containing 1) has been identified as interacting with the kinase domain of SRK (Gu et al., 1998;Stone et al., 1999). Thus, Bra031473 and Bra031431 gene products also have the potential to function as SRKdownstream molecules, similar to ARC1.
Protein kinases are important for signal transduction in higher plants (Shiu and Bleeker, 2001). In the Brassica SI system, MLPK (M-locus protein kinase) has been identified as a downstream factor of SRK from analysis of the SC mutant variety yellow sarson of B. rapa (Murase et al., 2004). Three genes, Bra031475, Bra031468 and Bra031416, encode protein kinase (super) family proteins that may also function as downstream factors of SRK and/or MLPK for signal transduction.
Thus, there are several candidate genes whose products appear to meet the requirements of a factor that functions downstream of the self/non-self recognition reaction in the Brassica SI signaling/response mechanism.
In conclusion, from this study we confirmed that the SC trait in TSC28 is regulated by a single recessive mutation which is independent from the S-locus or known SI-related genes, and we have narrowed down the genetic region containing the SC gene to a physical genomic length of approximately 1.01 Mbp. This region contains 113 genes and, among these, 52 genes (nine Chiifu-401-specific and 43 common genes) are candidates for the gene determining SC in TSC28, as identified by a combination of fine mapping and sequencing of papilla-expressed genes from the target genomic region. The list of candidate genes will contribute to discovery of the causal gene of SC in TSC28, coding for a novel downstream factor in the SP11-SRK signaling cascade in the Brassica SI system.