Comparative analyses of Brassica accessions have demonstrated that both genes relocation and sequence polymorphisms between species are common in the Brassica genome (Hu et al., 1998; Sebastian et al., 2000). Here, we show that the genome of cabbage and cauliflower are highly collinear, with macro collinearity punctuated by rearrangements mainly involving translocations and inversions. The majority of the observed rearrangements involved very short distances, but an apparent inversion spanning distance of more than 20 cM was noted at the central section of linkage group 04 between Marker13078 and Marker17566 (Figure S3; Table S3). Similar violation of genetic collinearity have also been identified within both intraspecific and interspecific level in other Brassica species, which may be responsible for the high degree of morphological polymorphism among these species (Sebastian et al., 2000; Bancroft et al., 2011). Since the cabbage and cauliflower

# REFERENCES

subspecies shared a polyphyletic origin in primitive B. oleracea populations (Truco et al., 1996), our results imply that the inverted region seem to have been involved in the divergence of cauliflower and cabbage. Although additional genomic data from other varieties/lines are necessary to trace the genome evolution, comparative information about the collinearity between these two closely related subspecies has important implications for further marker and gene identification in cauliflower through the use of sequence data from cabbage as a genomic model.

In conclusion, we demonstrated the application of SLAF-seq strategy in cauliflower for large-scale SNP discovery and highdensity genetic map construction. The parents used in this study are not only representatives of cultivated cauliflower, but also elite lines in our breeding program. Therefore, the genome-wide SNPs identified and the linkage map developed in this study will provide an important foundation not only for QTL identification and map-based cloning, especially of curd development-related traits, but also for MAS within cauliflower breeding programs.

#### AUTHOR CONTRIBUTIONS

ZZ and HG conceived and performed experiments and wrote the manuscript; XS, HY, and JW were involved in data analysis; LH and DW performed the SLAF library construction and sequencing. All the authors have commented, read and approved the final manuscript.

# FUNDING

This project is supported by the National Natural Science Foundation of China (grant number, 31501768); the Research project in Zhejiang Province Science and Technology department (grant number, 2012C12903-3-4, 2016C32102); the Research project in Zhejiang Academy of Agricultural Sciences (grant number, 2014CX031, 2015R23R08E05).

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2016. 00334

based on large scale marker development using next-generation double-digest restriction-site-associated DNA sequencing (ddRADseq). BMC Genomics 15:351. doi: 10.1186/1471-2164-15-351

**Conflict of Interest Statement:** 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.

Copyright © 2016 Zhao, Gu, Sheng, Yu, Wang, Huang and Wang. This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Genetic and Epigenetic Alterations of *Brassica nigra* Introgression Lines from Somatic Hybridization: A Resource for Cauliflower Improvement

Gui-xiang Wang1 , Jing Lv 1, 2, 3, Jie Zhang1 , Shuo Han1 , Mei Zong1 , Ning Guo1 , Xing-ying Zeng1 , Yue-yun Zhang1 , You-ping Wang2 and Fan Liu1 \*

*1 Beijing Vegetable Research Center, Beijing Academy of Agriculture and Forestry Sciences, Key Laboratory of Biology and Genetic Improvement of Horticultural Crops (North China), Ministry of Agriculture, Beijing, China, 2 Yangzhou University, Yangzhou, China, 3 Zhalute No.1 High School, Tongliao, China*

#### *Edited by:*

*Narendra Tuteja, International Centre for Genetic Engineering and Biotechnology, India*

#### *Reviewed by:*

*Anca Macovei, University of Pavia, Italy Carl Gunnar Fossdal, Norwegian Institute of Bioeconomy Research, Norway*

> *\*Correspondence: Fan Liu<EMAIL_ADDRESS>
#### *Specialty section:*

*This article was submitted to Crop Science and Horticulture, a section of the journal Frontiers in Plant Science*

*Received: 13 May 2016 Accepted: 08 August 2016 Published: 30 August 2016*

#### *Citation:*

*Wang G-x, Lv J, Zhang J, Han S, Zong M, Guo N, Zeng X-y, Zhang Y-y, Wang Y-p and Liu F (2016) Genetic and Epigenetic Alterations of Brassica nigra Introgression Lines from Somatic Hybridization: A Resource for Cauliflower Improvement. Front. Plant Sci. 7:1258. doi: 10.3389/fpls.2016.01258* Broad phenotypic variations were obtained previously in derivatives from the asymmetric somatic hybridization of cauliflower "Korso" (*Brassica oleracea* var. *botrytis*, 2*n* = 18, CC genome) and black mustard "G1/1" (*Brassica nigra*, 2*n* = 16, BB genome). However, the mechanisms underlying these variations were unknown. In this study, 28 putative introgression lines (ILs) were pre-selected according to a series of morphological (leaf shape and color, plant height and branching, curd features, and flower traits) and physiological (black rot/club root resistance) characters. Multi-color fluorescence *in situ* hybridization revealed that these plants contained 18 chromosomes derived from "Korso." Molecular marker (65 simple sequence repeats and 77 amplified fragment length polymorphisms) analysis identified the presence of "G1/1" DNA segments (average 7.5%). Additionally, DNA profiling revealed many genetic and epigenetic differences among the ILs, including sequence alterations, deletions, and variation in patterns of cytosine methylation. The frequency of fragments lost (5.1%) was higher than presence of novel bands (1.4%), and the presence of fragments specific to *Brassica carinata* (BBCC 2*n* = 34) were common (average 15.5%). Methylation-sensitive amplified polymorphism analysis indicated that methylation changes were common and that hypermethylation (12.4%) was more frequent than hypomethylation (4.8%). Our results suggested that asymmetric somatic hybridization and alien DNA introgression induced genetic and epigenetic alterations. Thus, these ILs represent an important, novel germplasm resource for cauliflower improvement that can be mined for diverse traits of interest to breeders and researchers.

Keywords: introgression lines, genetic diversity, *Brassica oleracea*, *Brassica nigra*, somatic hybridization, epigenetic variation, cauliflower

# INTRODUCTION

Cauliflower (Brassica oleracea var. botrytis, 2n = 18, CC genome) is a major vegetable crop valued worldwide for its nutrition and flavor. As with other modern crop species, intense selection for preferred traits, and a tendency for inbreeding have resulted in low genetic diversity among cauliflower breeding resources. To address this problem, related species such as black mustard (Brassica nigra, 2n = 16, BB genome), with a large reservoir of genes conferring desirable characteristics, have been proposed as a valuable source of genetic diversity for Brassica crop improvement. In principle, this diversity can be transferred into crops via sexual hybridization and subsequent backcrossing. However, in practice, genetic manipulation using distantly related plants has been severely restricted by difficulties in creating the initial sexual hybrid and by sterility issues in the early backcross generations (Glimelius, 1999; Jauhar, 2006). Asymmetric somatic hybridization is a potential alternative for gene transfer from relatives to cultivated crops, especially when wide crosses are not applicable (Gerdemann et al., 1994; Xia, 2009). The Brassicaceae is a model plant family commonly used for somatic hybridization (Glimelius, 1999; Navrátilová, 2004), resulting in interspecies or even intertribal hybrids. Moreover, agronomically important traits, such as disease resistance and specific fatty acid compositions, have been successfully integrated into the crops (Gerdemann et al., 1994; Hansen and Earle, 1994, 1995, 1997; Wang et al., 2003; Tu et al., 2008; Scholze et al., 2010).

Brassica is also an excellent model for the study of allotetraploid speciation. Brassica is closely related to the classic plant model Arabidopsis and a wealth of germplasm exists from species in this genus. Three diploid species are widely cultivated: Brassica rapa (AA, x = 10), B. nigra (BB, x = 8), and B. oleracea (CC, x = 9). Recent interspecific hybridization among these species has resulted in polyploidization and the production of three allotetraploid species: Brassica napus (AACC, x = 19), Brassica juncea (AABB, x = 18), and Brassica carinata (BBCC, x = 17) (Nagaharu, 1935).

The formation of allotetraploids is an influential mode of speciation in flowering plants, often accompanied by rapid, extensive genomic and epigenetic changes that globally alter gene expression, termed "genomic shock." These include fragment gain and loss through chromosome rearrangement or the activation of transposable elements (Kashkush et al., 2003; Kraitshtein et al., 2010), extensive alteration of DNA cytosine methylation (Zhao et al., 2011; Gautam et al., 2016), histone modification (Ha et al., 2011), and changes to small RNA (Ha et al., 2009). However, genomic-shock-induced changes during sexual polyploid synthesis should be distinct from the changes that occur from genomic shock during somatic hybridization. Somatic hybrids combine both the nuclear and cytoplasmic genomes within a single cell. Therefore, the introgression of chromatin segments via asymmetric somatic hybridization likely occurs via non-homologous end-joining of fragmented genome pieces, rather than by the homologous recombination that occurs through sexual reproduction (Liu et al., 2015). Moreover, the epigenetic states of somatic cells and gametal cells tend to differ, given that gametal cells are more conserved to ensure genetic fidelity (Bird, 1997, 2002; Liu et al., 2015). Thus, the variations induced by "somatic genomic shock" likely have unique genetic and epigenetic characteristics compared with "sexual genomic shock." However, little data exist on the exact nature of these differences, or if the types of genomic shock indeed differ in their effects.

A number of hybrid progenies have been regenerated from asymmetric somatic hybrids between cauliflower "Korso" and black mustard "G1/1." The hetero-cytoplasmic nature of these hybrids was confirmed through the finding that most progenies have chloroplast genomic components from "G1/1" and mitochondrial DNA from "Korso," while mitochondrial genome recombination occurred in a few hybrids (Wang G. X. et al., 2011). Among the progenies, dozens of lines containing 18 chromosomes showed some obvious characters of "G1/1" origin, such as waxless leaves with lobes and ears, as well as green petioles; however, they were also densely covered with trichomes and exhibited purple stems and leaf veins, characteristic of "Korso." This combination of characters indicated that they were putative introgression lines (ILs). Although hybrid synthesis primarily focused on the transfer of disease-resistance genes from "G1/1" into cultivated cauliflower (Wang G. X. et al., 2011) because the former is resistant to several Brassica pathogens (Westman et al., 1999), the putative ILs contained considerable genetic diversity for many traits beyond disease resistance. Therefore, these lines should possess many characters that differ from near-isogenic lines and feature discrete portions of "G1/1" chromatin in a "Korso" genetic background. Additionally, asymmetric somatic introgression should induce further diverse genomic variations involving change in DNA sequences and epigenetic modifications. To verify this, we used 12 putative ILs shown to be phenotypically stable for several generations to characterize genetic and epigenetic alterations from "somatic genomic shock."

### MATERIALS AND METHODS

#### Plant Materials

Asymmetric somatic hybridization between cauliflower "Korso" and black mustard "G1/1" resulted in the establishment of 117 individual hybrids in the field; of these, 13 fertile plants with preferred traits were selected for continued selfing and backcrossing (Wang G. X. et al., 2011). Hundreds of derivatives were obtained following year-by-year selection (since 2006) that combined phenotypic observation, cytological and molecular analysis, as well as pathogen-resistance assays (**Figure 1**). In this study, 28 putative ILs were chosen for phenotypic diversity and chromosomal constitution analysis, based on possession of cauliflower-like morphology and phenotypic stability over three generations. Next, 12 putative ILs (IL1-12) were used for genetic an epigenetic analysis. Chinese cabbage "Asko" (B. rapa, AA, 2n = 20) and tetraploid B. carinata (BBCC, 2n = 34) were chosen as reference materials for molecular and chromosomal analysis. The selfed seeds from ILs, parents, and B. carinata were sown in a greenhouse and are available upon request.

# Multi-Color FISH Analysis

Multi-color FISH was conducted on meiosis metaphase I chromosomes of parent lines and putative ILs to confirm their chromosomal composition. B. nigra and CentBr (centromericspecific tandem repeats of Brassica; accession numbers: CW978699 and CW978837, respectively; Wang G. et al., 2011) DNA were selected as probes and labeled with biotin (dig)-14-dUTP (Roche, Indianapolis, IN) using nick translation with an average length of 500 bp. Slides were prepared following a previously published protocol (Zhong et al., 1996) with minor

modifications. To decompose the cell walls of pollen mother cells, anthers approximately 1–3 mm in length were digested at 37◦C for 3 h in an enzyme mixture containing 6% cellulose R-10 and 6.5% pectinase (Sigma, solution in 40% glycerol). Previously described methods (Wang G. et al., 2011) were followed for FISH analysis. Hybridized probes were visualized using an FITC-conjugated avidin antibody or rhodamine-conjugated antidigoxin antibody (Roche, Indianapolis, IN). Chromosomes were counterstained using 0.1 mg/mL DAPI (Vector Laboratories, Burlingame, CA). Images of the signals and chromosomes were captured using a CCD camera (QIMAGING, RETIGA-SRV, FAST1394) attached to a Nikon Eclipse 80i epifluorescence microscope (Tokyo, Japan). Image contrast and brightness were adjusted in Adobe Photoshop (8.0).

#### SSR, AFLP, and MSAP Fingerprinting

Simple sequences repeat (SSR), amplified fragment length polymorphism (AFLP), and methylation-sensitive amplification polymorphism (MSAP) fingerprinting were performed following previously described methods (Shaked et al., 2001; Liu et al., 2015). Sixty-five SSR, 77 AFLP, and 16 MSAP primer combinations were used; for SSR and AFLP, markers were collected from the A genome of B. rapa due to the limited number of cytological and genetic studies on the B genome of B. nigra (Tables S1–S3). For AFLP and MSAP, two independent technical replicates were performed for each of the three biological replicates. Only clear and reproducible bands larger than 100 bp were scored.

# RESULTS

#### Phenotypic Diversity of Introgression Lines

The 28 ILs were extremely diverse in multiple morphological and physiological traits (**Table 1**), including distinct variations in leaf or flower morphology, fertility, and resistance to pathogens (i.e., black rot or clubroot). Here we describe diagnostic differences in flower heads as an example. Four phenotypic groups were apparent based on flower curd characters (**Figure 2**). In the first group (type I), flower heads were white, compact, and hemispheric, with a wheel-like radial arrangement; flower buds were small, granular, and close. The second group (type II) possessed flat, loose flower heads with no obvious wheel-like arrangement; flower buds were light yellow, tiny, and soft, while the pedicel was green and slightly longer than Type I pedicels. The third group (type III) possessed relatively flat flower heads and pale yellow-green color; flower buds were small, loose, soft, and distributed in clusters, but with no obvious ball flower formation. The pedicel was green and long. The flower heads of the fourth group (type IV) were relatively flat, yellow green, with small flower buds that were fine and soft.

# Chromosomal Constitution of the Putative Introgression Lines

**Figure 3** shows the results of multi-color FISH on chromosomes of "G1/1," "Korso," and 28 putative ILs. The two parent lines exhibited different hybridization patterns with both probes. As expected, when using "G1/1" genomic DNA, all 16 "G1/1" chromosome pairs exhibited red signals (**Figure 3A**). In "Korso," red signals (2 strong and 2 weak; **Figure 3B**) were detected on only two pairs of satellite chromosomes that likely correspond to the hybridization signals of 45SrDNA (Fukui et al., 1998).

No signals were detected on "G1/1" chromosomes using the CentBr probe. In contrast, seven pairs of strong signals and one pair of weak signals were detected on "Korso" chromosomes (**Figure 3B**); the latter was hybridized to chromosomes that also exhibited strong red signals. Only 1 "Korso" chromosome pair exhibited no signals (**Figure 3B**).

Most cells in the putative ILs contained 18 chromosomes and exhibited hybridization patterns identical to "Korso." The results of genomic in situ hybridization (GISH) found no "G1/1" chromosomal segments present in these putative ILs, likely because GISH cannot detect introgressions from Brassica chromosomes due to their small size and compact structure.

# Genetic Sequence Variation Analyzed Using SSR and AFLP Markers

We analyzed the sequence introgression and variation in 12 ILs (IL1–12) using 65 SSR and 77 AFLP markers that were polymorphic between the parents. The results of SSR and AFLP profiling revealed 1799 loci from both parents, including 682 (37.9%) those were polymorphic. The IL profiles were similar to the "Korso" profile, but loss of 5–13 "Korso"-specific bands were common. Additionally, all putative ILs were confirmed to possess anywhere from 3 to 33 "G1/1"-derived fragments. However, we noted that in 3 ILs, SSR markers failed to detect "G1/1"-derived bands but AFLP profiling revealed 20–25 "G1/1"-specific loci, indicating that the latter method is more effective when the introgression size is small. Next, all 12 ILs contained B. carinataspecific loci (35–49), and we were able to amplify 1–7 new bands from them.

TABLE 1 | Phenotypic traits associated with introgression lines (IL1-12) and their parents.

*(Continued)*

#### TABLE 1 | Continued

*(Continued)*

TABLE 1 | Continued

FIGURE 2 | Four types (I–IV) of representative flower heads in introgression lines. C: "Korso," B: "G1/1."

These results clearly showed that all 12 lines were ILs and that genomic-sequence variations were common among them. Notably, IL genomes were not simply similar to the "Korso" genome but with some "G1/1"-derived fragments added. For example, we noted that some ILs lacked bands present in "Korso," or possessed fragments lacking in "Korso." Moreover, several B. carinata-specific bands were present in all ILs, hinting at the occurrence of common changes in the early stages of both somatic hybridization and naturally occurring hybridization.

DNA profiling of the ILs revealed 1493 "Korso," 1423 "G1/1," and 1559 B. carinata fragments; 1117 of these were shared by both parents (**Figure 4**). The number of fragments present in the profiles of the 12 representative ILs ranged from 3 to 33, with an overall frequency of 7.5%, whereas the number of fragments lost ranged from 7 to 32, with an overall frequency of 5.1% (**Table 2**). Note that the percentages of "G1/1"-fragment presence and "Korso"-fragment loss were respectively based on the 306 "G1/1" fragments that did not co-migrate with "Korso" and the 376 "Korso" fragments that did not co-migrate with "G1/1." In addition, the 12 ILs possessed 1–7 new fragments (1.4%) that were not present in the "Korso" profile, and 4–49 (15.5%) B. carinata-specific fragments from the 226 that did not co-migrate with "Korso" and "G1/1."

# Epigenetic Changes in Methylation Patterns Analyzed with MSAP

The MSAP analysis used isoschizomers MspI and HpaII, which differ in their sensitivity to cytosine methylation (Shaked et al., 2001; Liu et al., 2015). Our comparison of the profiles generated from EcoRI-MspI and EcoRI-HpaII digestion revealed three kinds of band classifications (**Figure 5**). First, A type was characterized by fragments from both the EcoRI-MspI and EcoRI-HpaII digestions, representing non-methylated sites (or methylation within a single strand). Second, B type was characterized by half-methylation sites, present in EcoRI-HpaII digestion but not in EcoRI-MspI digestion. Third, C-type bands stemmed from complete methylation sites present in EcoRI-MspI digestion but not in EcoRI-HpaII digestion. Using 16 MSAP primer combinations amplifying clear and stable bands, we detected 401 fragments in parents, with 290 present in "Korso" but not in "G1/1."

**Table 3** summarizes variation in methylation patterns and their heredity in parents and the ILs. The average frequency of hypermethylated loci in different ILs was 12.4%, while the average frequency of hypomethylated loci was 4.8%.

# DISCUSSION

#### Somatic Hybridization and Introgression Can Potentially Increase Genetic Diversity

The low genetic diversity in modern Brassica varieties is concerning because it reduces potential genetic gains in breeding programs. To combat this problem, ILs containing fragments from related species can be used to generate improved cultivars.

In this study, we generated a set of ILs using somatic hybridization between cauliflower and black mustard. The original objective was to select disease- and drought-resistant ILs for use in a Brassica breeding program. In fact, these ILs showed considerable genetic diversity in many morphological and physiological traits beyond resistance. Some of these traits are useful for increasing cauliflower adaptation to stressors. For instance, IL2 and IL6 respectively exhibited strong resistance to black rot and clubroot, major Brassica diseases with worldwide distribution. Other ILs exhibited commercially desirable traits such as early maturation and high nutritional quality, especially high phosphorus and potassium content (IL3 and IL9), as well as less desirable traits such as loose flower heads, Early inflorescent, and long pedicel. This diversity indicates that all of the ILs could be a source of new alleles for developing Brassica cultivars. Thus, somatic hybridization provides an effective and efficient means of achieving introgression into a crop species from its wild relatives and should continue to be important in producing novel germplasm for breeding programs.

#### Combining Cytological and Molecular Methods Effectively Identified Brassica Introgression Lines

Although GISH has been successfully used to investigate genomic relationships in several agriculturally important genera, including Brassica (e.g., Lysak et al., 2005), and has detected introgressions in crops such as wheat (e.g., Liu et al., 2015), GISH was unable to detect black mustard chromosomal segments in our ILs. We believe this outcome was likely due to the small and compact Brassica chromosomes or else very small introgression sizes. We dealt with this issue by first using multi-color FISH to identify chromosome composition and then combining SSR and AFLP markers to detect "G1/1" sequence introgression. Primer sources and amplification results from the genomes of B. rapa, B. nigra, B. oleracea, and B. carinata are shown in Table S1. As those data indicate, the primers are usable in multiple Brassica species and can detect abundant polymorphism. Moreover, compared with the SSR primers, AFLP markers were more effective when working with small introgressions, such as those in our ILs (**Table 2**). In sum, our study demonstrated that Brassica ILs can be clearly identified with a combination of cytological and molecular methods.

#### TABLE 2 | DNA profiles and phenotypic summary of partial introgression lines (IL1-12).

#### TABLE 2 | Continued

### Genetic Changes from Somatic Hybridization in the Brassica Introgression Lines

The extensive phenotypic variation of the ILs was caused by genetic and epigenetic changes in the ILs relative to their parents. However, black mustard introgression was actually more limited (7.5%; **Table 2**) than other ILs, such as the somatic hybrids between wheat and tall wheatgrass (Liu et al., 2015). We believe the low introgression was probably due to suboptimal UV treatment. Thus, treatment factors such as irradiation time and intensity will require further testing and validation to achieve higher introgression rates. At the same time, the amount of introgression observed may be enough to achieve considerable genomic diversity. For example, in stable rice introgression lines derived from intergeneric hybridization (followed by successive selfing) between rice (Oryza sativa L.) and a wild relative (Zizania latifolia Griseb.), very low introgression (0.1%) by foreign DNA was enough to trigger up to 30% of the genomic changes (Wang et al., 2005). The mechanism underlying these changes is unclear, but some researchers have suggested that they happen very early, through a cryptic pathway that differs from conventional or unorthodox meiotic recombination of homeoalleles between rice and Zizania (Zhang et al., 2013).

Similar to previous findings (Wang et al., 2005), we also observed sequence loss and new bands (**Table 2**) following somatic hybridization, backcrossing, and self-pollination (**Figure 1**). The frequency of fragments lost (5.1%) was higher than presence of novel bands (1.4%). It is not entirely clear why the percentage of novel bands was so much lower, but one reason may be that some of them were B. carinata-specific, TABLE 3 | Variation in cytosine methylation across the partial introgression lines (IL1-12).

and B. carinata-specific fragments occurred at a relatively high frequency (15.5%). Thus, future investigations may need to confirm the exact nature of the lost and novel bands, as well as exclude any B. carinata-specific loci.

Our results were consistent with previous reports in artificial hybrid B. carinata (BBCC, x = 17), produced via crossing and polyploidization between B. nigra (BB, x = 8) and B. oleracea (CC, x = 9). Around 47% of the hybrid B. carinata (BBCC) genome possessed isoenzyme sites specific to natural allotetraploid B. carinata (Jourdan and Salazar, 1993). Given

FIGURE 4 | Amplified fragment length polymorphism (AFLP) fingerprinting of the parents and ILs, with primer pair E45–M57. C: "Korso," B: "G1/1," BC: *B. carinata*, 1–11: introgression lines. Arrows indicate "Korso"-specific bands, diamond-headed pointers indicate "G1/1"-specific bands, and circle-headed pointers indicate *B. carinata*-specific bands.

FIGURE 5 | Methylation-sensitive amplified polymorphism (MSAP) profiling based on primer combinations ETTC + HMTCC. 1–7: introgression lines, B: "G1/1," C: "Korso," H: DNA cleaved with EcoRI + HpaII, M: DNA cleaved with EcoRI + MspI. A, B, and C indicate three separate band types.

this similarity, we propose that the early generations of artificial synthesis and natural polyploid evolution of both genomes may have experienced the same change events. Supporting this idea, previous analyses of recombination in other genera, such as the homoploid hybrid sunflower (Helianthus anomalus), showed that newly synthesized hybrids converged on the linkage pattern of wild hybrids within five generations (Rieseberg et al., 1996).

Based on our results and those of previous studies, we propose that genetic changes in Brassica ILs occur via the following mechanisms. First, alien genetic elements are incorporated into cauliflower early in hybridization, through the activation of transposable elements or DNA methylation, as suggested in rice introgression lines (Wang et al., 2005). Next, chromosomal rearrangements occur through intergenomic translocations or transpositions (Udall et al., 2005) and homoeologous pairing (Leflon et al., 2006; Nicolas et al., 2007). During Brassica evolution, ancestral group rearrangements contributed to the high homology found especially in the three linkage groups of the B genome (B5, B6, and B4) and the A genome (A5, A6, and A4) (Panjabi et al., 2008). Finally, consecutive selection results in rapid sequence elimination, as reported in synthetic hybrids and hybrids less related to the original parents (Osborn et al., 2003; Pires et al., 2004; Gaeta et al., 2007; Zhang et al., 2013).

#### Changes to DNA Methylation in the Brassica Introgression Lines

Cytosine methylation plays an important role in epigenetic gene regulation at both the transcriptional and the posttranscriptional levels (Paszkowski and Whitham, 2001). Our ILs exhibited fairly high proportions of change to methylation patterns (17.2%; **Table 3**), compared with several other plant hybrid/allopolyploid systems analyzed to date using MSAP. For instance, in resynthesized allotetraploid Arabidopsis suecica, 8.3% of the fragments experienced methylation changes (Madlung et al., 2002), whereas in newly synthesized allohexaploid wheat, ∼13% of the loci saw alterations to cytosine methylation (Shaked et al., 2001). In contrast, a 23.6% change to cytosine methylation was observed in asymmetric somatic wheat introgressions (Liu et al., 2015). Taken together, these data allow us to conclude that somatic hybridization induced a broad spectrum of cytosine-methylation changes that perturbed gene expression to a larger extent than allopolyploidization.

Consistent with previous findings (Zhang et al., 2013), changes to hypermethylation (9.7–13.5%; **Table 3**) occurred more frequently than changes to hypomethylation (2.4–7.6%). Zhang et al. (2013) proposed that this pattern may be caused by the activation and subsequent silencing of some transposable elements after hybridization, but this hypothesis requires more data for confirmation. Additionally, the methylation patterns in ILs may be related to changes in enzymatic machinery, as well as the expression of siRNAs and long noncoding RNAs (Goll and Bestor, 2005; Wierzbicki et al., 2008; Zhang et al., 2013).

In conclusion, somatic hybridization mimics many of the genetic alterations induced by polyploidization or sexual wide hybridization, but to a stronger extent and with considerably less time. Therefore, somatic hybridization is both effective and efficient in achieving introgression of a crop species and its wild relatives. Moreover, this approach provides a potential means to explore the genetic and epigenetic events induced by "somatic genomic shock" (Liu et al., 2015).

#### REFERENCES

# AUTHOR CONTRIBUTIONS

GW conceived the research work and wrote the paper. FL provided plant materials and intellectual advice on the project. JL performed molecular and FISH experiments. JZ and SH performed the field management and character survey. XZ performed the methylation-sensitive amplified polymorphism analysis. NG participated in the technical guidance of relative experiment. MZ and YZ performed materials classification and preservation. YW provided technical guidance for pathogenresistance assays.

# ACKNOWLEDGMENTS

This work was supported by the Natural Science Foundation of China (No. 31000538), the Natural Science Foundation of Beijing, China (No. 6142012), the Project of Technology Innovation Ability from Beijing Academy of Agriculture and Forestry Sciences (KJCX20140111), and the Youth Science Research Foundation of Beijing Academy of Agriculture and Forestry Sciences (No. QNJJ201601). National Key Research and Development Project (2016YFD010 0204-14).

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2016. 01258

(Festuca pratensis) into Italian rye grass (Lolium multiflorum) and physical mapping of the locus. Heredity 91, 396–400. doi: 10.1038/sj.hdy.6800344

**Conflict of Interest Statement:** 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.

Copyright © 2016 Wang, Lv, Zhang, Han, Zong, Guo, Zeng, Zhang, Wang and Liu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Genome-specific differential gene expressions in resynthesized *Brassica* allotetraploids from pair-wise crosses of three cultivated diploids revealed by RNA-seq

Dawei Zhang1 , Qi Pan1 , Cheng Cui 2 , Chen Tan1 , Xianhong Ge1 , Yujiao Shao3 \* and Zaiyun Li 1 \*

#### *Edited by:*

*Naser A. Anjum, University of Aveiro, Portugal*

#### *Reviewed by:*

*Genlou Sun, Saint Mary's University, Canada Maoteng Li, Huazhong University of Science and Technology, China Naghabushana K. Nayidu, University of Saskatchewan, Canada*

#### *\*Correspondence:*

*Yujiao Shao<EMAIL_ADDRESS>Zaiyun Li<EMAIL_ADDRESS>
#### *Specialty section:*

*This article was submitted to Crop Science and Horticulture, a section of the journal Frontiers in Plant Science*

*Received: 31 July 2015 Accepted: 20 October 2015 Published: 04 November 2015*

#### *Citation:*

*Zhang D, Pan Q, Cui C, Tan C, Ge X, Shao Y and Li Z (2015) Genome-specific differential gene expressions in resynthesized Brassica allotetraploids from pair-wise crosses of three cultivated diploids revealed by RNA-seq. Front. Plant Sci. 6:957. doi: 10.3389/fpls.2015.00957* *1 National Key Lab of Crop Genetic Improvement, National Center of Crop Molecular Breeding Technology, National Center of Oil Crop Improvement (Wuhan), College of Plant Science and Technology, Huazhong Agricultural University, Wuhan, China, 2 Crop Research Institute, Sichuan Academy of Agricultural Sciences, Chengdu, China, 3 College of Chemistry and Life Science, Hubei University of Education, Wuhan, China*

Polyploidy is popular for the speciation of angiosperms but the initial stage of allopolyploidization resulting from interspecific hybridization and genome duplication is associated with different extents of changes in genome structure and gene expressions. Herein, the transcriptomes detected by RNA-seq in resynthesized *Brassica* allotetraploids (*Brassica juncea*, AABB; *B. napus*, AACC; *B. carinata*, BBCC) from the pair-wise crosses of the same three diploids (*B. rapa*, AA; *B. nigra*, BB; *B. oleracea*, CC) were compared to reveal the patterns of gene expressions from progenitor genomes and the effects of different types of genome combinations and cytoplasm, upon the genome merger and duplication. From transcriptomic analyses for leaves and silique walls, extensive expression alterations were revealed in these resynthesized allotetraploids relative to their diploid progenitors, as well as during the transition from vegetative to reproductive development, for differential and transgressive gene expressions were variable in numbers and functions. Genes involved in glucosinolates and DNA methylation were transgressively up-regulated among most samples, suggesting that gene expression regulation was immediately established after allopolyploidization. The expression of ribosomal protein genes was also tissue-specific and showed a similar expression hierarchy of rRNA genes. The balance between the co-up and co-down regulation was observed between reciprocal *B. napus* with different types of the cytoplasm. Our results suggested that gene expression changes occurred after initial genome merger and such profound alterations might enhance the growth vigor and adaptability of *Brassica* allotetraploids.

Keywords: *Brassica* species, allopolyploidization, differential gene expressions, transgressive gene expression, transcriptome

# INTRODUCTION

Allopolyploidization which is realized through the merger and duplication of distinct parental genomes following interspecific hybridizations of two or more related species results in the origin of new allopolyploid species (Otto, 2007; Doyle et al., 2008). This pattern of speciation occurs widely in angiosperms, largely because the allopolyploids relative to their progenitors show the enhanced growth vigor (the phenomenon of heterosis) and the advantage in ecological adaptation (Chen, 2007; Leitch and Leitch, 2008). The obvious success of allopolyploidy in nature has allured the extensive investigations of genetic effects caused by genome merger at levels of the chromosomes, DNA sequences, gene expression, proteins, small RNA, and DNA methylation during the last two decades, by utilizing the continuously improved approaches of molecular biology and genome sequencing (Song et al., 1995; Soltis and Soltis, 2012; Li et al., 2014; Soltis et al., 2014). The results from recent and synthetic allopolyploids of different taxa (Arabidopsis, Brassica, cotton, wheat, triticale, etc.) demonstrate that the initial stage of allopolyploidization is accompanied by the various genetic, epigenetic and transcriptome changes, but the degrees of variations are obviously different between allopolyploids (Soltis and Soltis, 2012). So the new allopolyploids respond to the genome merger by rapid alterations in genomic structure and function, in order to coordinate the divergent genome at different aspects and to establish the novel plants for further evolution. As to the crucial gene expression in hybrids and new allopolyploids, widespread changes are revealed by transcriptomic analysis and extensive gene expression changes are non-additive in allopolyploids relative to their parents, including expression level dominance and transgression (outside the range of either parent), such profound changes to gene expression may enable new hybrids to survive in novel environments not accessible to their parent species (Wang et al., 2006; Chen, 2007; Rapp et al., 2009; Flagel and Wendel, 2010; Yoo et al., 2013; Li et al., 2014).

The six cultivated Brassica species offer a text book system of allopolyploidization through the pair-wise crosses of three diploids, which is illustrated as U-triangle (U 1935). Brassica carinata (2n = 34, BBCC), B. juncea (2n = 36, AABB), and B. napus (2n = 38, AACC) are three allotetraploids, which are derived from various two-way combinations of the three diploids B. nigra (2n = 16, BB), B. oleracea (2n = 18, CC), and B. rapa (2n = 20, AA). Resynthesized Brassica allotetraploids through the interspecific hybridizations between three extant Brassica diploid progenitors have been widely investigated to elucidate the genetic alterations during the initial stage of allopolyploid formation, since the seminal work of Song et al. (1995). The reciprocal synthetics of allotetraploids would remedy the drawback of unprecise progenitors for natural counterparts and also uniparental cytoplasm background. The results obtained provide many new insights into the genome regulations following genome merger, though the studies focused mainly on the younger B. napus with much more agricultural and economic importance (Gaeta et al., 2007; Nicolas et al., 2007, 2012; Xu et al., 2009; Szadkowski et al., 2010, 2011; Xiong et al., 2011; Cui et al., 2012, 2013; Ge et al., 2013). The genome of B. rapa has been sequenced, which should enhance the evolutionary analysis of these Brassica allotetraploids (Wang et al., 2011). By using the B. rapa as reference genome, the analysis on resynthesized B. napus across four generations revealed that the gene expression was more complicated than the simple combination of two genomes, and non-additive gene regulation was also detected (Jiang et al., 2013). Another transcriptomic study on synthesized Brassica allohexaploid and its parents showed that genome-wide changes in gene expression were involved in adaptation and evolution processes, and non-additive genes associated with important biological processes were identified (Zhao et al., 2013). An intriguing genetic interaction in three Brassica allopolyploids was the hierarchy of nucleolar dominance (genomes BB> AA> CC) in three allotetraploids, in which both B. juncea (BB> AA) and B. carinata (BB> CC) expressed the rRNA genes from B. nigra, and B. napus (AA> CC) expressed those from B. rapa, but the genes from another parent were silenced (Chen and Pikaard, 1997; Ge et al., 2013). Meanwhile, the rRNA genes silenced in vegetative tissues were expressed in reproductive tissues, indicated that the expression of rRNA genes could be tissue-specific (Chen and Pikaard, 1997). However, whether the parent-specific expression of the genes for ribosomal protein also occurs is still an open question.

In this study, the gene expressions detected by RNA-seq in three resynthesized Brassica allotetraploids prior to meiosis from the pair-wise crosses of the same three diploids are compared, with the aims: (1) to reveal the differential gene expressions from progenitor genomes after genome merger and duplication and with the exclusion of meiotic homoeologous exchanges; (2) to assess transgressive gene expression and its effect on adaption environment changes in Brassica allotetraploids; (3) to identify the expression of ribosomal protein genes after genome merger and the relationship with the nucleolar dominance. The results should give new insights into the genetic contribution of three progenitors and their interactions at the beginning of Brassica allopolyploidization.

#### MATERIALS AND METHODS

#### Plant Materials and Sample Preparations

Four Brassica allotetraploids (AACC/CCAA, AABB, and BBCC) used in this study were previously synthesized from pair-wise crosses of three cultivated diploids, B. rapa L. (2n = 20, AA genome, genotype 3H120), B. nigra (L.) Koch cv. Giebra (2n = 16, BB), B. oleracea var. alboglabra L. (2n = 18, CC, genotype Chi Jie Lan), with the aid of embryo rescue and colchicine inducing chromosome doubling in vitro (Cui et al., 2012, 2013). Three allotetraploids (AABB, AACC, and BBCC) were synthesized by doubling the chromosome numbers of the respective digenomic hybrids, and CCAA was directly obtained from the cultured embryo derived plantlets, probably by the spontaneous chromosome doubling in vitro. The plants of each synthetic were derived from single hybrid embryo by successive

**Abbreviations:** DEGs, differentially expressed genes; FC, fold change; GO, gene ontology.

subculturing on MS medium, and those of each parent were derived from one seed. The rooted plantlets of four allotetraploids and three parents cultured on the medium were transferred to the pots and kept in the unheated greenhouse. Leaves at the same stage were collected from these seven materials, while the silique walls after 21 days of pollination were collected. All samples were stored in liquid nitrogen and kept at −80 until RNA extraction.

### RNA-seq Library Preparation and Sequencing

Total RNAs were extracted by TRIzol reagent (invitrogen) and were then treated with RQ1 DNase (promega) to remove DNA. The quality and quantity of the purified RNA were determined by measuring the absorbance at 260/280 nm (A260/A280) using smartspec plus (BioRad). RNA integrity was further verified by 1.5% agrose gel electrophoresis.

For each sample, 10µg of total RNA was used for RNA-seq library preparation. Polyadenylated mRNAs were purified and concentrated with oligo(dT)-conjugated magnetic beads (invitrogen) before used for directional RNA-seq library preparation. Purified mRNAs were iron fragmented at 95◦C followed by end repair and 5′ adaptor ligation. Then reverse transcription was performed with RT primer harboring 3′ adaptor sequence and randomized hexamer. The cDNAs were purified and amplified and PCR products corresponding to 200– 500 bps were purified, quantified, and stored at −80◦C until used for sequencing.

For high-throughput sequencing, the libraries were prepared following the manufacturer's instructions and applied to illumina GAIIx system for 80 nt single-end sequencing by ABlife. Inc (Wuhan, China). Owing to the high cost for sequencing 2 years ago, the total of 14 (2 tissues × 7 plant samples) libraries without replicates were made for RNA-seq.

#### Reads Filtering and Alignment

After the high-throughput sequencing, the raw data which contained adapter and low quality cycle of reads were filtered. Due to lack of reference genome from B. nigra and the genome sequence of B. oleracea was still unpublished yet 2 years ago (Liu et al., 2014), the clean reads were aligned to the B.rapa var. pekinensis Chiifu-401 genome as described in previse study (Jiang et al., 2013; Zhao et al., 2013) using the default parameters for TopHat, allowing no more than two mismatched bases. Subsequently, only unique mapped reads were used in further study to provide sensitive and accurate alignment results, even for highly repetitive genomes. The gene expression level was calculated by using RPKM method (Reads Per kb per Million reads). If there were more than one transcript for a gene, the longest one was used to calculate its expression level and coverage.

### Analysis of Differentially Expressed Genes (DEGs)

Program edgerR was used to analyze the differentially expressed genes between samples. Fold change ≥ 2 and P ≤ 0.01 between compared samples were considered as simulated biological variation in DEG analysis. GO enrichment was performed using blast2Go (http://www.blast2go.com/b2ghome). In addition, the web-based Brassica database (http://brassicadb.org/brad/index. php) and WEGO (Ye et al., 2006) server for gene ontology analysis was also used in this study (P < 0.05). We also used DAVID to investigate the transgressively expressed genes, GO terms with enrichment score ≥ 0.5 and P < 0.05 were considered significantly enriched (Huang et al., 2008).

Pearson correlation between biological samples was calculated using IBM SPSS V19.0 software. Additionally, the Venn diagrams in this study were performed using online tool (http://bioinfogp. cnb.csic.es/tools/venny/).

#### Analysis of r-protein Genes

All the mapped genes were blasted to the Arabidopsis genome, those genes matching the Arabidopsis ribosomal protein genes were considered to be r-protein genes in Brassica. To group rprotein genes with similar expression patterns, a hierarchical clustering was generated using the normalized expression values (log2 RPKM) from each library. The analysis was conducted using HemI software with Pearson correlation as the distance measure (Deng et al., 2014).

### Real-time Quantitative RT-PCR (qRT-PCR) Analysis

The RNA samples used for the qRT-PCR assays were the same as for the RNA-seq experiments. First-strand cDNA synthesis was performed with 1500 ng of total RNA using Thermo Scientific RevertAid First Strand cDNA Synthesis Kit, total RNA (0.5µg) was reverse-transcribed with oligo (dT)18 primer (0.5µg/µl) according to the described protocol. Gene-specific primers were designed according to the reference unigene sequences using the Primer 3.0, all primer sequences are listed in **Supplementary Table 5**. A primer was also designed for B. napus actin gene to normalize the amplification efficiency. QRT-PCR assays in triplicate were performed using Kapa Probe Fast qPCR Kit with a Bio-Rad CFX96 Real-Time Detection System. The actin gene was used as an internal control for data normalization, and quantitative variation in the different replicates was calculated using the delta-delta threshold cycle relative quantification method.

# Availability of Supporting Data

The datasets supporting the results of this article are available in GenBank SRA under accession ID PRJNA281555. Other supporting data are included within the article and its additional files.

# RESULTS

#### Gene Expressions between Nascent Allotetraploids and Parents

To investigate and compare the mRNA expression levels in the resynthesized Brassica allotetraploids relative to their diploid progenitors (**Figure 1**), 14 RNA-seq libraries were constructed for two types of tissue: leaves (L. for short) and silique walls (S.),

respectively. As a result, we obtained an average of 12,748,770 (84.02%) high-quality and clean reads from the raw reads (**Supplementary Table 1**). Among the clean reads, the average 41.76% reads were matched either to unique or multiple genomic positions using the B. rapa reference genome. As most of uniquely mapped reads were mature mRNA or ncRNA and multiple mapped reads were mainly rRNA and tRNA, only uniquely and perfectly mapped reads (3,696,642) were used to measure the transcriptional activity of each gene. Finally, our RNA-Seq data revealed an average of 28944 genes in Brassica allotetraploids and its progenitors, accounting for 70.56% of the total genes in B. rapa reference genome. For further comparative analysis, the gene expression level was calculated using RPKM method (Reads Per kb per Million reads).

A correlation dendrogram showed that the gene expression patterns were distinct between leaves and silique walls, for the samples of the two tissues were clustered separately (**Figure 2**). However, the global relationships of gene expression among the synthetics and their parents were the same in two tissues. In the dendrogram, AACC and CCAA which were most closely correlated were clustered firstly with AA and then with AABB, and BBCC was clustered with CC, while BB appeared as the outlier.

#### Differentially Expressed Genes (DEGs) in Nascent Allotetraploids

To study the gene expression patterns during allopolyploidization process, we first performed pair-wise comparisons between allotetraploids to their parents to identify differentially expressed genes (DEGs) using edgerR (fold change ≥ 2 and P ≤ 0.01 as criteria; **Supplementary Figure 1**). As a result, average 3791 (13.1% of expressed genes) DEGs were up-regulated and 3534(12.2% of expressed genes) DEGs were down-regulated in allotetraploids, respectively (**Figure 1**). There was no significant difference between the average number of up and down-regulated genes in leaves (3593 vs. 3721; t-test, P > 0.05), however, the difference was statistically significant (3990 vs. 3347; t-test, P < 0.05) in silique walls, suggesting that the direction of DEGs was affected by tissue type. The maximum number (8898) of DEGs was observed between L.CCAA and L.CC among all comparisons, including 4780 up-regulated and 4118 down-regulated. But the minimum number (5368) of DEGs was between S.CCAA and S.CC, only 2774 genes were up-regulated and 2594 were down-regulated.

Notably, by comparing the total number of genes showing differential expression (both up and down regulation) between the allotetraploids and their parents, there was a bias in the direction of differential expression relative to the parents (**Table 1**). For example, more expressed genes remained statistically unchanged (less differential expression) between L.AACC/CCAA and L.AA than between L.AACC/CCAA and L.CC (Chi square test, P < 0.01). This asymmetric gene expression was also observed in L.AABB and L.BBCC, the global expression patterns were closer to either L.AA in L.AABB or L.CC in L.BBCC (Chi square test, P < 0.01). Whereas, the expression bias in silique walls was not as tangible as in leaves, the global expression patterns were closer to either S.BB in S.AABB or S.CC in L.AACC.

#### Identification of Transgressively Expressed Genes

Among the differentially expressed genes, we then filtered the transgressively expressed genes in allotetraploids which showed more than two-fold changes in expression relative to both parents (**Table 2**). Briefly, on average 878 genes (3.0% of total expressed genes) were transgressively up-regulated while 652 genes (2.3% of total expressed genes) were transgressively down-regulated. When we compared the number of genes for up and downregulation in each sample, no significant bias was observed in leaves, but there were more genes exhibiting transgressive upregulation in silique walls.

Moreover, the comparisons of transgressively expressed genes in each Brassica allotetraploids revealed that transgressive expression varied between tissues (**Figure 3**). For instance, more genes showed transgressive up-regulation in silique walls than leaves in AABB (1172 vs. 566). Among those genes, a multitude of them were specifically expressed in different direction and tissues, suggesting that the majority of transgressively expressed genes were tissue-specific. In addition, certain number of genes showed co-transgressive up- or down-regulation (105 and 40) regardless of tissue type, likely due to genome merger. However, 28 genes showing transgressive down-regulation in leaves displayed up-regulation in silique walls, while only nine genes showed expression changes in the opposite direction (From up in leaves to down in silique walls; **Figure 3**). The similar tendency of gene expression changes was observed in other Brassica allotetraploids (28 vs. 6 in BBCC; 21 vs. 5 in CCAA; 28 vs. 9 in AACC), suggesting that genes in silique walls tended to be transgressively up-regulated.

Transgressively expressed genes were further functionally classified according to Gene Ontology (GO) terms using

DAVID (Huang et al., 2008). Ten functional clusters with the highest enrichment score were selected with criteria that enrichment score > 0.5 and P < 0.05(see more details in **Supplementary Table 2**). Although remarkable difference was observed among samples, it was interesting to note that a certain functional cluster encoding glucosinolates was up-regulated among a group of allotetraploids (L.BBCC, L.CCAA, L.AACC, S.AABB, S.AACC), suggesting that the increased expression of resistance-related genes, especially glucosinolates, was common in Brassica allotetraploids (**Supplementary Table 2**).

#### Co-transgressive Gene Expression in *Brassica* Allotetraploids

To further ascertain which genes were specially expressed in tissues regardless of genome composition and what the function of those genes was, we investigated genes showing cotransgressive expression in leaves and silique walls, respectively. To avoid confusion, three representative allotetraploids (AABB, BBCC, and CCAA) with different cytoplasm were selected for further comparative analysis.

As illustrated in the Venn diagrams, 17 genes exhibited cotransgressive up-regulation and 18 genes gave co-transgressive TABLE 1 | Summary of differential expression between the *Brassica* allotetraploids and their parents.

*aThe total number of genes showing differential expression between the allotetraploids and their maternal parent, including both up and down-regulated genes.* \**Chi square test with expected ration of 1:1, P* < *0.01.*

down-regulation among the allotetraploids in leaves, respectively (**Figures 4A,B**). Interestingly, two key methylase genes (Bra002610 and Bra022537) which were involved in DNA methylation on cytosine and methyltransferase activity reflected co-transgressive up-regulation in the allotetraploids, suggesting that methylation-related genes were up-regulated during Brassica allopolyploidization in leaves (**Supplementary Table 4**). In silique walls, there were more genes displaying co-transgressive up-regulation than down-regulation (23 vs. 5, **Figures 4C,D**).

#### TABLE 2 | Transgressive expressions in *Brassica* allotetraploids.

*aCalculated by dividing the number of expressed genes in each Brassica allotetraploid.* \**Chi square test with expected ration of 1:1, P* < *0.01.*

#### Novel Gene Expression and Silencing in Allotetraploids

To explore novel gene expression and gene silencing in the allotetraploids, a strict criterion was set: novel expression was defined when both parental lines had no reads but in allotetraploids and the genes displayed more than 10 reads and RPKM ≥ 2. On the contrary, if both parental lines contained genes with more than 10 reads and RPKM ≥ 2, yet allotetraploids had no reads, we considered this situation as gene silencing.

Overall, a total of 160 genes were found to be novel expression, while 102 genes exhibited silencing in allotetraploids, suggesting that these changes were tightly controlled (**Table 3**). There was no significant difference between the total number of genes showing novel expression and silencing in both tissues, but considerable variation between the two different expression patterns was found in silique walls: more genes exhibited silencing than novel expression in S.CCAA, whereas more genes showed novel expression than silencing in S.BBCC and S.AABB. Furthermore, the number of genes showing novel expression in each sample was negatively correlated with those genes being silenced (Pearson correlation, r = −0.92, P = 0.008 < 0.01; **Figure 5**).

The biological functions and processes of the two expression patterns were analyzed using gene ontology annotations according to function annotation convention (**Supplementary Figure 2**). The enriched GO terms were similar between tissues, while the number of genes was various in the enriched GO terms. The genes associated with cell part, cellular, as well as metabolic process and response to stimulus were over-represented in both tissues. However, there were more genes showing novel expression than silencing in those GO terms in leaves, but more genes exhibited silencing than novel expression in silique walls.

#### Differentially Expressed Genes between Reciprocal *B. napus*

To evaluate how gene expression was affected by cytoplasm types, we examined gene expression changes between the reciprocal synthetics of B. napus (AACC/CCAA) which kept the same nuclear genomes but the cytoplasm from B. rapa or B. oleracea. Though a large proportion of genes showed differential expression, the total number of up-regulated genes was similar to that of down-regulated genes (6210 vs. 6270), suggesting that there was no bias in the direction of differential expression between AACC and CCAA (Chi square text,χ 2 = 0.29, P > 0.01; **Table 4**). Among these differentially expressed genes, only those showing co-regulation were selected for subsequent function analysis to avoid the tissues specific effects, the difference between co-up (581) and co-down (530) regulation was also insignificant (Chi square text, χ 2 = 2.25, P > 0.01; **Table 4**).

To have an overview of the major differences between the reciprocal synthetics, we performed GO analysis of the differential expressed genes. We found that the number of co-up regulated genes was similar to that of co-down regulated in the most enriched classes (**Figure 6**). Among the cellular component categories, cell and cell parts were the largest enrichment, followed by organelle. As to the molecular function categories, the most over-represented GO terms were binding and catalytic, far more than others. In the biological process categories, there were more GO terms compared to the other two categories, and slightly more genes involved in the biological regulation, cellular process, metabolic process as well as response to stimulus were found to be up-regulated in these GO terms.

genes in silique walls.

TABLE 3 | Genes showing novel expression and silencing in three *Brassica* allotetraploids.

*a*,*bCalculated by dividing the number of expressed genes in each Brassica allotetraploid.*

#### Analysis of r-protein Gene Expression

Besides the afore-mentioned expression patterns, we further identified ribosomal protein genes among samples to evaluate whether the genome merger has similar effect on housekeeping genes. We compared the Arabidopsis r-proteins against the above mapped genes based on the sequencing data. 363 genes were identified in our libraries which matched a total of 79 of the 80 r-proteins identified in Arabidopsis, suggesting that there was a high number of homologous r-proteins genes between Brassica and Arabidopsis (**Supplementary Table 4**; Whittle and Krochko, 2009). However, the expression of r-protein genes varied extensively across 14 libraries, the average expression level of r-protein genes in leaves was generally higher than that in silique walls (**Supplementary Figure 3**). The highest expression level of r-proteins genes were detected in L.AABB and L.BBCC, however, it was notable that the lowest expression level were also found in their silique walls (S.AABB and S.BBCC).

To further reveal the relationships of r-proteins gene expression among Brassica allotetraploids and diploids in different tissue, hierarchical clustering of expression level (standardized by log2 RPKM) was conducted for the 363 expressed r-proteins genes (encoding 79 r-proteins) using HemI (Deng et al., 2014). Those Brassica allotetraploids and diploids with similar expression profiles of r-proteins genes were more closely clustered in leaves and silique walls, respectively (**Figures 7**, **8**). In leaves, two allotetraploids with B-genome (L.AABB and L.BBCC) were closer to each other and clustered with L.BB, while reciprocal B. napus (L.AACC/L.CCAA) exhibited similar transcription profiles with L.AA and diverged markedly from L.CC. In silique walls, r-protein genes retained similar expression pattern to that in leaves on whole genome scale. Three allotetraploids with A-genome (S.AACC, S.CCAA, and S.AABB) were closer to S.AA and separated from S.CC and S.BBCC, while all of them were quite distinct from S.BB.

#### Verification of DEGs by qRT-PCR

To confirm that the above results of differentially expressed genes, a set of gene-specific primers were designed for quantitative RT-PCR assays (**Supplementary Table 5**). The relative transcript levels in the allotetraploids and their parents were compared with those of RNA-seq data (RPKM value). For 24 out of the 30 comparisons, qRT-PCR analysis revealed

\**Chi square test with expected ration of 1:1, P* < *0.01.*

the same expression trends as the RNA-seq data, despite some quantitative differences (**Supplementary Figure 4**). Among 25 genes showing transgressive expression in the allotetraploids relative to their parents using RNA-Seq, 12 were up-regulated, eight were down-regulated, and five were differentially expressed to only one of their parents, respectively. Moreover, four of five genes showing differential expression between the reciprocal synthetics of B. napus also showed the same expression trends revealed by qRT-PCR and RNA-seq, confirming the reliability of RNA-seq data.

#### DISCUSSION

Widespread changes of gene expression during allopolyploidization have been revealed in different aspects, including homoeolog expression bias, novel gene expression, or silencing, transgressive up or down-regulation, expression level dominance, and altered expression times and locations (Doyle et al., 2008). In this study, global transcriptome analyses through RNA-Seq were made for serial synthetic Brassica allotetraploids derived from the same three diploid parents (**Figure 1**; Cui et al., 2012). We focused on general analysis of transgressive gene expression, novel gene expression, and gene silencing, r-protein genes expression upon Brassica allopolyploidization. In addition, qRT-PCR analysis was also performed to confirm the reliability of RNA-seq data.

# Effects of Transgressive Gene Expression during Allopolyploidization

In consistence with other findings of gene expressions in cotton (Yoo et al., 2013), wheat (Li et al., 2014), and Senecio (Hegarty et al., 2008), transgressive gene expression in our Brassica allotetraploids was genome-wide and temporal. A lot of genes exhibited opposite expression direction between tissues, suggesting that transgressively expressed genes had different roles in development (**Table 2**). Additionally, transgressively up-regulated gene expression increased over time, more genes were or reversed to be transgressively up-regulated in silique walls (**Figure 3**). Moreover, wide ranges of alterations occurred for the functional clusters of transgressively expressed genes between tissues (**Supplementary Table 2**). All these changes indicated that transgressive gene expression was sophisticated in different genome background and different growth phases. In particular, transgressive up-regulation of resistance-related genes including for the synthesis of glucosinolates could be responsible for immediate physiological pre-adaptation of Brassica allotetraploids, since glucosinolates and their breakdown products played a role in the defense of plants against pathogen (Brader et al., 2001), fungi (Hiruma et al., 2010), and insects (Ahuja et al., 2015).

Interestingly, the genes which involved in DNA methylation were transgressively up-regulated in leaves across three types of Brassica allotetraploid (**Figure 4A**; **Supplementary Table 3**). DNA methylation was an important and best-studied epigenetic phenomenon in allopolyploids. Allopolyploids (e.g., wheat, Senecio, Spartina) underwent rapid and widespread metylation state changes and significant proportion of these changes showed non-additivity in previous studies (Dong et al., 2005; Salmon et al., 2005; Hegarty et al., 2011). Recent investigations in synthesized Brassica hybrids and allotetraploids also provided evidence for rapid changes in DNA methylation by transcriptome analysis or genetic analysis (Xu et al., 2012; Cui et al., 2013; Jiang et al., 2013; Zhao et al., 2013). Consistent with these studies, our data further demonstrated that the expression of methylation-related genes was non-additive and indeed transgressively up-regulated. DNA methylation was vitally important in silencing transposons and regulating gene expression, typically reducing expression (Zilberman et al., 2007; Zhang, 2008), as a consequence, the reduced expression and transposon silencing could prove capacity for adaptation.

Together, these data suggested that transgressive gene expression could be various in number and function between tissues. Such variation (including DNA methylation) might play a causative role in regulation of gene expressions in allopolyploids and increase their relative fitness over their parents in novel environment.

# Negative Correlation between Novel Gene Expression and Silencing

Evidence for rapid gene silencing and novel expression came from the detection of missing parental fragments by cDNA - AFLP screens and verification by RT-PCR in Arabidopsis and cotton, and these expression changes could be variable in different parts of the plants (Adams et al., 2003; Wang et al., 2004). Recently, reciprocal gene silencing and novel expression were also found both in natural and synthetic cottons by using RNA-Seq, and specifically increased in natural allopolyploids (Yoo et al., 2013).

To explore patterns of novel gene expression and silencing during the Brassica allopolyploidization, synthetic Brassica allotetraploids (U 1935) were examined using RNA-Seq data, since RNA-Seq technology could provide more data than cDNA -AFLP and RT-PCR. By restrictive criteria, the percentage of genes showing novel expression and silencing was significantly less than that of Arabidopsis and cotton using AFLP-cDNA and RT-PCR method (Adams et al., 2003; Wang et al., 2004), while similar to that of cotton from RNA-Seq (Yoo et al., 2013). This difference may be explained by technical consideration and various restrictive criteria. As a check in present study, we have also validated of five genes using qRT-PCR, four genes and one gene, respectively, were confirmed as novel expression and silencing (**Supplementary Figure 4**). Our results also found a wide range of variations between the two different expression patterns and different tissues of the allotetraploids (**Table 4**). More genes displayed novel expression or silencing in silique walls than in leaves, suggesting that tissue-specific expression partitioning could arise quickly after the onset of allotetraploids formation, which was consistent with the previous study (Adams et al., 2003; Wang et al., 2004; Adams and Wendel, 2005; Buggs

et al., 2009). Interestingly, there was no significant difference between the total number of genes, but inverse correlation between the two expression patterns was observed among samples, indicating that there was a tradeoff between novel gene expression and silencing during plant development (**Figure 5**).

#### Cytoplasmic Effects on Genes Expression in Reciprocal *B. napus*

The cytoplasm has exhibited considerable influence on the evolution of nuclear genomes of allopolyploids (Prakash et al., 2009), because the presence of the paternal nuclear genome in the maternal cytoplasm could result in nuclear-cytoplasmic incompatibilities. Extensive and rapid genomic changes as well as the variations in chromosome meiotic pairings were observed in the reciprocal Brassica allopolyploids, but the gene expression changes were not as obvious as the genomic changes (Song et al., 1995; Cui et al., 2012, 2013). Although large-scale differences in gene expression existed between the reciprocal synthetics of B. napus (AACC and CCAA), the numbers of genes showing co-up and co-down regulation gave insignificant difference (**Table 4**). It was interesting to see that these co-regulated genes (co-up and co-down) were enriched in the similar functional classes with similar numbers, revealing that certain gene networks may be particularly susceptible by cytoplasm (**Figure 6**). Overall, the equivalence in either number or function of co-regulated genes indicated that there might be a balance of differentially expressed genes in nuclear-cytoplasmic interactions.

#### Tissue-specific and Differential Expression of r-protein Genes

Although the ribosome was well-known to be an immensely complex "molecular machine" in translating the genetic code, recent studies indicated that the r-protein genes were likely involved in tissue-specific process and had a regulatory role in the development of plant (Byrne, 2009; Whittle and Krochko, 2009). Our data indicated that those sets of r-protein genes which represented a wide range of r-proteins were differently expressed between tissues (**Supplementary Table 4**). In general, the average expression level of r-protein genes in leaves was higher than that in silique walls (**Supplementary Figure 3**). Nevertheless, an early finding showed that r-protein gene expression was greatest among highly differentiating reproductive tissues (not including silique walls), for example, microspores, embryos, and seeds, using available EST data. However, the lowest levels were also detected in some reproductive tissues, such as anthers, pollen, seed coats (Whittle and Krochko, 2009). Our results demonstrated that the r-protein gene expression levels in silique walls were lower than that in leaves and such difference could be attributable to the fact that silique walls were more mature and the translation activity was not as high as in young leaves. Additionally, in two Brassica allotetraploids with B-genome (AABB and BBCC), the expression level of r-protein genes which were highest among samples in leaves turned to be lowest in silique walls. An recent study in the same synthesized allotetraploids showed higher number and percentage of absent bands of B genome than that of A or C genome in AABB and BBCC using AFLP-cDNA, suggesting that gene expression was particularly susceptible to perturbation when B-genome was present in Brassica allotetraploids during plant development (Cui et al., 2013).

Hierarchical clustering of these genes indicated that the expression relationships among Brassica allotetraploids and diploids were diverse in different tissues (**Figures 7**, **8**). The hierarchy of r-protein genes in leaves (B > A > C) was consistent with that of the expression of rRNA genes (Chen and Pikaard, 1997). However, the expression of rRNA genes could be developmentally regulated and tissue-specific, because the rRNA genes silenced in vegetative tissues were found to be expressed in reproductive tissues, including sepals and petals (Chen and Pikaard, 1997). Thus, the expression divergence of r-protein

#### REFERENCES

Adams, K. L., Cronn, R., Percifield, R., and Wendel, J. F. (2003). Genes duplicated by polyploidy show unequal contributions to the transcriptome and organ-specific reciprocal silencing. Proc. Natl. Acad. Sci. U.S.A. 100, 4649–4654. doi: 10.1073/pnas.0630618100

genes could be associated with the expression change of parental rRNA genes.

In summary, widespread changes of gene expression were observed among the serial resynthesized Brassica allotetraploids relative to their diploid progenitors. There were considerable alterations on temporal and spatial expressions, and the range of variations in silique walls was much wider than in leaves. The expression of r-protein genes was tissue-specific and associated with nucleolar dominance. Furthermore, novel gene expression was negatively correlated with silencing during the transition from vegetative to reproductive development. The balance between the transgressive up- and down-expressions was also observed in leaves, as well as between the co-up and co-down regulation in different cytoplasm, such profound changes might enhance the fitness and adaptability of Brassica allopolyploids.

#### ACKNOWLEDGMENTS

This work was funded by National Natural Science Foundation of China (Grant No. 31375656) and by Special Fund for Agro-Scientific Research in the Public Interest (201203026-6).

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2015. 00957

Supplementary Figure 1 | Analysis of the differentially expressed genes between each allotetraploid and its parents. (A) Leaves. (B) Slique walls.

Supplementary Figure 2 | GO functional categories of genes showing novel expression and silencing. (A) Leaves. (B) Slique walls.

Supplementary Figure 3 | The average expression level of total r-protein genes among samples.

Supplementary Figure 4 | RT-PCR confirmation of the differentially expressed genes. Columns and bars represent the means and standard error (*n* = 3), respectively. The gene expression levels from RNA-Seq data are added on the top of each gene.

Supplementary Table 1 | Summary of alignment statistics of RNA-Seq in the 14 samples.

Supplementary Table 2 | Top 10 GO items of transgressively regulated genes.

Supplementary Table 3 | Annotation of co-transgressively regulated genes among allotetraploids.

Supplementary Table 4 | R-protein genes matched each *Arabidopsis* r-proteins.

Supplementary Table 5 | The corresponding primers of qRT-PCR.

**Conflict of Interest Statement:** The reviewer Genlou Sun declares that, despite having previously collaborated with the authors Xianhong Ge, Yujiao Shao, and Zaiyun Li, the review process was handled objectively. 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.

Copyright © 2015 Zhang, Pan, Cui, Tan, Ge, Shao and Li. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# The miRNAs and their regulatory networks responsible for pollen abortion in Ogura-CMS Chinese cabbage revealed by high-throughput sequencing of miRNAs, degradomes, and transcriptomes

Xiaochun Wei 1 † , Xiaohui Zhang2 †, Qiuju Yao1 , Yuxiang Yuan1 , Xixiang Li 2 , Fang Wei 3 , Yanyan Zhao1 , Qiang Zhang1 , Zhiyong Wang1 , Wusheng Jiang1 and Xiaowei Zhang1 \*

#### Edited by:

*Naser A. Anjum, University of Aveiro, Portugal*

#### Reviewed by:

*Maoteng Li, Huazhong University of Science and Technology, China Hao Peng, Washington State University, USA*

#### \*Correspondence:

*Xiaowei Zhang<EMAIL_ADDRESS>
*† These authors have contributed equally to this work.*

#### Specialty section:

*This article was submitted to Crop Science and Horticulture, a section of the journal Frontiers in Plant Science*

Received: *29 July 2015* Accepted: *08 October 2015* Published: *21 October 2015*

#### Citation:

*Wei X, Zhang X, Yao Q, Yuan Y, Li X, Wei F, Zhao Y, Zhang Q, Wang Z, Jiang W and Zhang X (2015) The miRNAs and their regulatory networks responsible for pollen abortion in Ogura-CMS Chinese cabbage revealed by high-throughput sequencing of miRNAs, degradomes, and transcriptomes. Front. Plant Sci. 6:894. doi: 10.3389/fpls.2015.00894* *1 Institute of Horticulture, Henan Academy of Agricultural Sciences, Zhengzhou, China, 2 Key Laboratory of Biology and Genetic Improvement of Horticultural Crops, Ministry of Agriculture, Institute of Vegetables and Flowers, Chinese Academy of Agricultural Sciences, Beijing, China, 3 College of Life Science, Zhengzhou University, Zhengzhou, China*

Chinese cabbage (*Brassica rapa* ssp. *pekinensis*) is one of the most important vegetables in Asia and is cultivated across the world. Ogura-type cytoplasmic male sterility (Ogura-CMS) has been widely used in the hybrid breeding industry for Chinese cabbage and many other cruciferous vegetables. Although, the cause of Ogura-CMS has been localized to the orf138 locus in the mitochondrial genome, however, the mechanism by which nuclear genes respond to the mutation of the mitochondrial orf138 locus is unclear. In this study, a series of whole genome small RNA, degradome and transcriptome analyses were performed on both Ogura-CMS and its maintainer Chinese cabbage buds using deep sequencing technology. A total of 289 known miRNAs derived from 69 families (including 23 new families first reported in *B. rapa*) and 426 novel miRNAs were identified. Among these novel miRNAs, both 3-p and 5-p miRNAs were detected on the hairpin arms of 138 precursors. Ten known and 49 novel miRNAs were down-regulated, while one known and 27 novel miRNAs were up-regulated in Ogura-CMS buds compared to the fertile plants. Using degradome analysis, a total of 376 mRNAs were identified as targets of 30 known miRNA families and 100 novel miRNAs. A large fraction of the targets were annotated as reproductive development related. Our transcriptome profiling revealed that the expression of the targets was finely tuned by the miRNAs. Two novel miRNAs were identified that were specifically highly expressed in Ogura-CMS buds and sufficiently suppressed two pollen development essential genes: *sucrose transporter SUC1* and *H* +*-ATPase 6*. These findings provide clues for the contribution of a potential miRNA regulatory network to bud development and pollen engenderation. This study contributes new insights to the communication between the mitochondria and chromosome and takes one step toward filling the gap in the regulatory network from the orf138 locus to pollen abortion in Ogura-CMS plants from a miRNA perspective.

Keywords: miRNAs, Brassica rapa ssp. pekinensis, Ogura-CMS, bud, pollen, deep sequencing

# INTRODUCTION

Cytoplasmic male sterility (CMS) is a maternally inherited trait that results in the loss of the ability to produce fertile pollen; CMS has been used extensively for hybrid crop breeding (Woodson and Chory, 2008; Luo et al., 2013). Ogura-CMS was originally discovered in the wild radish (Raphanus sativus) (Ogura, 1968) and has been widely introduced into Cruciferous vegetables and oil crops and successfully applied in the heterosis breeding industry due to its complete male sterility and stability (Pelletier et al., 1983). Ogura-CMS is controlled by a mitochondrial orf138 locus that consists of two co-transcribed open reading frames: orf138 and atp8 (ATP synthase, subunit 8) (Krishnasamy and Makaroff, 1993; Grelon et al., 1994). However, the mechanism by which the orf138 locus results in male sterility is unclear. Many researchers have focused on the crosstalk between the mitochondria and nucleus, and the CMS is believed to be controlled by mitochondrial-nuclear interactions (Chase, 2007). Several CMS-associated nuclear genes, such as those related to programmed cell death (PCD) and reactive oxygen species (ROS), have been identified and proposed to play roles in the CMS pathway (Balk and Leaver, 2001). Other regulatory signaling pathways, such as the biogenesis of jasmonic acid, are also impaired in the CMS line (Liu et al., 2012).

MicroRNAs (miRNAs) are a class of 21–24 nt small non-coding RNAs that regulate gene expression by posttranscriptional repression (Carrington and Ambros, 2003; Bartel, 2004). In plants, miRNAs are generated from primary miRNA transcripts (pri-miRNAs) that are generally transcribed by RNA polymerase II (Bartel, 2004; Khraiwesh et al., 2010). Primary transcripts containing a distinctive hairpin structure are first trimmed by DICER-LIKE1 (DCL1) to generate miRNA precursors (pre-miRNAs) in the nucleus. Then, these premiRNAs are transported to the cytoplasm and further processed by DCL1 to generate ∼21 nt mature miRNAs (Kurihara and Watanabe, 2004). Subsequently, the mature miRNAs are loaded onto the RNA-induced silencing complex (RISC) and guide RISC recognition of complementary sites on target mRNAs, thereby inducing transcript cleavage (Bartel, 2004; Baulcombe, 2004) or translational repression (Aukerman and Sakai, 2003; Chen, 2004). Plant miRNAs mainly function by directing the cleavage of their highly complementary target transcripts. Thus, it is easy to predict and verify their targets using bioinformatics and experimental methods.

miRNAs play important roles in the regulation of a wide range of plant developmental processes, including plant architecture (Zhang et al., 2011b), leaf development (Floyd and Bowman, 2004; Juarez et al., 2004), root development (Guo et al., 2005), tuberization (Bhogale et al., 2014), the vegetative-reproductive phase change (Wang et al., 2011a), flowering time (Zhou et al., 2013), floral organ identity (Aukerman and Sakai, 2003; Bartel, 2004; Chen, 2004; Nagpal et al., 2005), self-incompatibility (Tarutani et al., 2010), plant nutrient homeostasis (Yamasaki et al., 2007; He et al., 2014), and response to environmental biotic and abiotic stresses (Juarez et al., 2004; Navarro et al., 2006; Sunkar et al., 2007; Katiyar-Agarwal and Jin, 2010; Zhang et al., 2011a). However, few miRNAs have been confirmed to control pollen development or CMS. Efforts have been made to screen for pollen-specific or enriched miRNAs in rice (Wei et al., 2011), Arabidopsis (Chambers and Shuai, 2009; Grant-Downton et al., 2009a,b), and pakchoi (Jiang et al., 2014).

Chinese cabbage (Brassica rapa ssp. pekinensis) is an important leafy vegetable that is a cross-pollinated crop with significant heterosis. Ogura-CMS was transferred to the Chinese cabbage in the 1980s and is still widely used in the breeding industry (Yamagishi and Bhat, 2014). The miRNAs from the Ogura-CMS and the maintainer line buds of a closely related subspecies (B. rapa ssp. chinensis) have been profiled (Jiang et al., 2014). In that study, 54 new conserved miRNAs and 25 pairs of novel miRNAs were identified. Among them, 18 miRNAs were differentially expressed between the male sterile and fertile lines. A genome-wide analysis of miRNAs in B. rapa by deep sequencing has been reported (Kim et al., 2012). However, to the best of our knowledge the miRNAs that underlie flower bud development and respond to the Ogura-CMS in Chinese cabbage have not been examined.

In present study, the miRNAs from the buds of a Chinese cabbage Ogura-CMS line (Tyms) and its maintainer line (231- 330) were profiled by small RNA deep sequencing. The targets of the miRNAs were identified by degradome sequencing analysis. The expression levels of the corresponding targets were monitored by transcriptome sequencing. These results provide a full set of miRNA-target messages that underlie flower bud development in the Chinese cabbage and provide insights into miRNA-related regulatory functions that may play roles in the mitochondrial–nuclear interactions in Ogura-CMS.

#### MATERIALS AND METHODS

#### Plant Materials and RNA Extraction

The Chinese cabbage Ogura-CMS sterile line (Tyms) and its maintainer line (231-330) used in this study were grown in the Henan Academy of Agricultural Sciences, Yuanyang, Henan Province, China. The lines possessed isogenic chromosomes with different cytoplasmic genes. Flower buds <6 mm length were stripped from 10 individual of each genotype of plants. These buds contain the entire developmental progress from anther generation to pollen abortion. Tyms and 231-330 were sampled separately and then snap-frozen in liquid nitrogen and kept at - 80◦C for further use. Total RNA was extracted using the TRIzol reagent (Invitrogen, USA). DNase (Promega, USA) was used to remove potential DNA contamination.

#### Small RNA Library Construction and Sequencing

The RNA samples from 231-330 (Control) and Tyms (CMS) were quantified and equalized so that equivalent amounts of RNA were analyzed. A total of 30µg of RNA was resolved on denatured polyacrylamide gels. Gel fragments with the size range of 18–30 nt were excised and recovered. These small RNAs were ligated with 5′ and 3′RNA adapters with the T4 RNA ligase. The adapterligated small RNAs were subsequently transcribed into cDNA by Super-Script II Reverse Transcriptase (Invitrogen) and amplified using primers specific for the ends of the adapters. The amplified cDNA products were purified and finally sequenced using Solexa sequencing technology (BGI, Shenzhen, China).

# Identification of Known and Novel miRNAs

The adapter sequences, impurities, and sequences with less or more than 18–30 nt were filtered out from the raw sequence reads. The remaining sequences that ranged from 18 to 30 nt in length were used for known miRNA prediction. First, we aligned the tags to the B. rapa miRNA precursors in miRBase (version 20.0, http://www.mirbase.org/index.shtml) with no mismatches allowed. Then, the obtained tags were aligned to the mature miRNAs of B. rapa with at least a 16 nt overlap to allow offsets. The miRNAs that satisfied both of the above criteria were counted to obtain the expression number of identified miRNAs. The rest of the small RNA tags were aligned to the miRNA precursors/mature miRNAs of all plants in miRBase, allowing two mismatches and free gaps. The miRNAs with the highest expression levels for each mature miRNA family were chosen as a temporary miRNA database. Then, the precursors of the temporary miRNAs in the B. rapa genome were predicted. Those that failed to fold into hairpin structures were regarded as pseudo-miRNAs and discarded, while those that fulfilled the miRNA criteria were adopted as newly identified known miRNAs.

By comparing our sequences with those in the databases and the Chinese cabbage genome, the sRNAs can be annotated into different categories, including siRNA, piRNA, rRNA, tRNA, snRNA, snoRNA, repeat associated sRNA, degraded tags of exons or introns, and sRNAs that could not be annotated. The tags annotated as intron, exon antisense, and unknown were used to predict novel miRNAs using the software Mireap. The key conditions are as follow: hairpin miRNAs can fold into secondary structures and mature miRNAs are present in one arm of the hairpin precursors; the 5-p and 3-p mature miRNAs present 2-nucleotide 3′ overhangs; hairpin precursors lack large internal loops or bulges; the secondary structures of the hairpins are steady, with the free energy of hybridization lower than or equal to −18 kcal/mol; and the copy number of mature miRNAs with predicted hairpins must be greater than five in the alignment result. The expression of novel miRNAs was produced by summing the count of miRNAs with no more than three mismatches on the 5′ and 3′ ends and no mismatches in the middle from the alignment result.

The differentially expressed miRNAs were calculated with the following procedures. First, the expression of miRNAs in two samples was normalized to obtain the expression of transcript per million (TPM). The normalization formula was as follows: Normalized expression = Actual miRNA count/Total count of clean reads × 1,000,000; Second, the fold-change and P-value were calculated from the normalized expression.

# Degradome Library Construction and Target Identification

Total RNA extracted from the 231-330 (Control) and Tyms (CMS) lines. Approximately 200µg of the total RNA was polyadenylated using the Oligotex mRNA mini kit (Qiagen). A 5′ RNA adapter was added to the cleavage products (which possessed a free 5′ -monophosphate at their 3′ termini) using the T4 RNA ligase (Takara). Then, the ligated products were purified using the Oligotex mRNA mini kit (Qiagen) for reverse transcription to generate the first strand of cDNA using an oligo dT primer via SuperScript II RT (Invitrogen). After the cDNA library was amplified for 6 cycles (94◦C for 30 s, 60◦C for 20 s, and 72◦C for 3 min) using Phusion Taq (NEB), the PCR products were digested with the restriction enzyme Mme I (NEB). A double-stranded DNA adapter was ligated to the digested products using T4 DNA ligase (NEB). The ligated products were selected based on size in a 10% polyacrylamide gel and purified for the final PCR amplification (94◦C for 30 s, 60◦C for 20 s, and 72◦C for 20 s) for 20 cycles (German et al., 2009). The PCR products were gel purified and used for high-throughput sequencing with the Illumina HiSeq 2000.

Low quality sequences and adapters were removed, and the unique sequence signatures were aligned to the database of Chinese cabbage transcript assemblies in the Brassica Gene Index (http://brassicadb.org/brad/downloadOverview.php) using the SOAP software (Li et al., 2008) (http://soap.genomics.org.cn/). CleaveLand was used to detect potentially cleaved targets based on degradome sequences. The 20 and 21 nt distinct reads were subjected to the CleaveL and pipeline for small RNA target identification as previously described (Addo-Quaye et al., 2009). The tags mapped to cDNA sense strands were used to predict cleavage sites. The miRNA-mRNA pairs were searched and pvalues were calculated using PAREsnip (http://srna-workbench. cmp.uea.ac.uk/tools/paresnip/). Only p-values less than 0.05 were adopted for the t-plot figure (Folkes et al., 2012). All alignments with scores not exceeding 4 that possessed the 5′ end of the degradome sequence coincident with the 10th and 11th nucleotides that were complementary to the small RNA were retained. To evaluate the potential functions of miRNAtargeted genes, gene ontology (GO) categories (http://www. geneontology.org/) were used to assign the identified target genes according to the previously described method (Du et al., 2010).

#### Transcriptome Sequencing and Target Gene Profiling

Total RNA (10µg) was subjected to poly-A selection, fragmentation, random priming, first and second strand cDNA synthesis with the Illumina Gene Expression Sample Prep kit (CA, USA). The cDNA fragments were subjected to an end repair process and then ligated to adapters. The products were enriched with PCR, and the 200-bp fragments were purified with 6% TBE PAGE gel electrophoresis. After denaturation, the single-chain fragments were fixed onto the Solexa Sequencing Chip (Flowcell) and consequently grown into single-molecule cluster sequencing templates through in situ amplification on the Illumina Cluster Station. Double-end pyrosequencing was performed on the Illumina Genome Analyzer platform with read lengths of 101 bp for each end. The clean reads were aligned to the Chinese cabbage reference transcript assemblies (http://brassicadb.org/ brad/downloadOverview.php). Gene expression levels were calculated using the RPKM method (Mortazavi et al., 2008) using the following formula: RPKM = (1,000,000 ∗ C)/(N ∗ L ∗ 1000), where RPKM(A) is the expression of gene A, C is the number of reads that uniquely align to gene A, N is the total number of reads that uniquely align to all genes, and L is the number of bases in gene A. Statistical comparison between Tyms and 231–330 was performed using the IDEG6 software (Romualdi et al., 2003). The General Chi-squared method was used and the FDR (false discovery rate) was applied to determine the Q-value threshold. Unigenes were considered to be differentially expressed when the RPKM between Tyms and 231–330 displayed a more than two-fold change with an FDR less than 10−2 .

#### Quantitative Real-time PCR

For analysis of miRNAs, 2.5µg of total RNA was polyadenylated using a miRNA cDNA synthesis kit (Takara, Inc., Dalian, China). The poly(A)-tail-amended total RNA was reversetranscribed by PrimeScript RTase using a universal adapter primer containing oligo-dT. The qPCR was performed on a LightCycler 96 System (ROCHE, USA) using SYBR Premix Ex TaqTMII (TaKaRa, Dalian China). The miRNA-specific forward primer for each miRNA was designed based on the entire miRNA sequence (Table S13), and the universal reverse primer was provided by the miRNA cDNA synthesis kit (Takara, Dalian, China). All reactions were performed with three biological and three technical replicates for each sample, and the U6 snRNA (Forward: GGGGACATCCGATAAAATT, Reverse: TGTGCGTGTCATCCTTGC) was used as the internal control. The reaction volume was 20µL, including 10µL of SYBR Premix Ex Taq II, 0.8µL of 10 mM Forward primer, 0.8µL of 10 mM Uni-miR qPCR primer, 2.0µL of the cDNA sample and 6.4µL of dH2O. The following qPCR program was used: denaturation at 95◦C for 30 s, followed by 40 cycles of 95◦C for 5 s, 55◦C for 30 s, and 72◦C for 60 s. Melting curve analysis with 95◦C for 10 s, 65◦C for 60 s, and 97◦C for 1 s was performed to produce a dissociation curve for verification of the amplification specificity. Relative expression levels of miRNAs were quantified using the 2−11Ct method (Livak and Schmittgen, 2001).

The primers of the selected genes subjected to target analysis are listed in Table S14. β-actin was used as an internal control. Experiments were performed on a similar system as described above. The reaction volume was 20µL, including 10µL of SYBR R Premix Ex Taq™ (Tli RNaseH Plus), 0.8µL of 10 mM Forward primer, 0.8µL of 10 mM Reverse primer, 2.0µL of the cDNA sample and 6.4µL of dH2O. Three independent biological and three technical replicates were performed. The fold change was estimated using the 2−11CT method (Livak and Schmittgen, 2001).

# Paraffin Sectioning and Microscopic Observation

Sterile and fertile flower buds were fixed and embedded in Paraffin (Beeswax, China). Thin (0.8µm) sections were prepared with an Ultracut Eultra microtome (Leica, Germany), stained with hematoxylin, and photographed under a LEICA DMI3000B microscope (Leica, Germany).

# RESULTS AND DISCUSSION

# sRNA Sequencing and miRNA Identification

To analyze the roles of miRNAs in Ogura-CMS in Chinese cabbage, small RNAs were pyrosequenced from the buds of the Ogura-CMS Tyms line and its maintainer line 231– 330 (the morphology and microscopy of the samples are shown in **Figure 1**). A total of 22.2 and 24.4 M reads were produced from the two lines, respectively. After filtering out the low quality reads, 3′ adapter null, insert null, 5′ adapter contaminants, reads smaller than 18 nt and poly A, a total of 21.9 and 24.1 M high quality clean reads were obtained, respectively (Table S1). The majority of the tags ranged in size between 21–24 nt, with the 24 and 21 nt lengths dominant (**Figure 2**). Out of the 46.0 M total tags, 33.8 M (73.5%) were shared by the two samples, while 5.7 M (12.5%) and 6.4 M (14.0%) were Tyms- and 231–330-specific, respectively. Then, the tags were divided into 12.2 M non-redundant unique reads composed of 4.8 M (39.6%) Tyms-specific, 5.6 M (45.7%) 231– 330-specific and 1.8 M (14.8%) reads shared by both lines (Table S2).The reads were aligned to B. rapa miRNAs (B. rapa 1.1) in miRBase (http://www.mirbase.org/); among the 157 mature miRNAs encoded by 96 precursors distributed into 63 families included in the database, 69 mature miRNAs from 59 precursors of 41 families were expressed in our samples (Table S3). The remainder of the reads were aligned to all known plant miRNAs in miRBase 21.0 and then aligned to the B. rapa genome for precursor identification. From this analysis, another 220 mature miRNAs derived from 163 precursors were identified. Among them, 57 precursors contained both 5p and 3p miRNAs, while only one mature miRNA was detected for the remaining 106 precursors. The 139 (39 pairs and 61 singles) mature miRNAs encoded by 100 precursors were new members of 36 existing families in the B. rapa miRBase. Thus, only 17 families (bra-MIR9552—bra-MIR9557 and bra-MIR9559—bra-MIR9569) in the miRBase B. rapa collection were not detected in our present study. By expending the members to the existing B. rapa miRNA families, bra-miR156 was now the largest family (harboring 21 members), followed by bra-miR171 (11 members), bra-miR167 (8 members), bramiR172 (8 members), bra-miR164 (7 members), bra-miR168 (7 members), bra-miR2111 (7 members), bra-miR157 (6 members), bra-miR160 (6 members), bra-miR390 (6 members), and bramiR395 (6 members). The other families contained less than five members (Table S4). In addition to these existing families, we also identified 81 miRNAs (18 pairs and 45 singles from 63 precursors) belonging to 23 families that had not been previously reported in B.rapa. These included 13 members of bra-miR1439, seven members of bra-miR169, five members each of bra-miR166, bra-miR393, bra-miR394, and bra-miR399, three bra-miR165, two bra-miR170, two bra-miR5376, and one member each of bra-miR397, bra-miR827, bra-miR828, bramiR838, bra-miR845, bra-miR858, bra-miR5298, bra-miR5575, bra-miR5641, bra-miR6029, bra-miR6030, bra-miR6033, bramiR6034, and bra-miR6284 (Table S5). The remaining reads

were aligned to Genbank, Rfam and the cabbage genome and annotated as rRNA, scRNA, snRNA, snoRNA, tRNA, repeat, exon sense, exon antisense, intron sense, intron antisense, and unannotated sequences (Table S6). The unannotated reads and those derived from the intron region and exon antisense region were used for novel miRNA analysis. A total of 426 novel miRNAs were identified based on the criteria defined in the Methods section. The sequences and precursor information were listed in Table S7. Among them, 83 novel miRNAs were generated from more than one precursor in the Chinese cabbage. A total of 138 precursors (84 types of mature miRNA) harbored both the 3-p and 5-p miRNAs on their arms (Supplementary file 2).

Of these miRNAs, members of 32 known families were also identified in buds of a closely related vegetable plant B. rapa ssp. chinensis (Jiang et al., 2014). These miRNAs may be common bud development contributors in both plants. Our present study identified more miRNAs of either known or novel type than that of the previous report. One important reason for this is that nearly fourfold of reads were generated in present study compared to the previous one. This facilitated the identification of low abundant miRNAs. Another reason can be attributed to the genetic differences between the two subspecies. One of our sequencing projects showed only half of reads from B. rapa ssp. chinensis could map to the B. rapa ssp. pekinensis reference genome (unpublished data). The publicly released B. rapa ssp. pekinensis genome (Wang et al., 2011b) is not an ideal reference for B. rapa ssp. chinensis. Jiang et al. (2014) used B. rapa ssp. pekinensis reference genome for their miRNA annotation and novel-miRNA prediction, there for many miRNAs were failed to be identified due to the sequence diversity, especially for the novel miRNAs which rely on the hairpins prediction.

#### Comparative Expression Patterns of miRNAs in the Buds of Ogura-CMS and its Maintainer

The expression levels of known and novel miRNAs were profiled based on tag counts. For multiple precursors sharing the same mature miRNA, only one member of the mature miRNA was used for expression profiling. The miRNA with the highest expression levels was bra-miR157, which accumulated 579,691 and 253,845 copies in the 231-330 and Tyms buds, respectively, followed by bra-miR168 and bra-miR156 (Table S3). bra-miR156 and bra-miR157 belong to the same family and regulate the SQUAMOSA promoter binding protein-like (SPL) genes that are important regulators for plant vegetative and reproductive development as well as gynoecium differential patterning and male fertility (Wang et al., 2009; Wu et al., 2009; Xing et al., 2010, 2013; Zhang et al., 2011b; Yu et al., 2012). bra-miR168 targets ARGONAUTE 1 (AGO1), which is an important part of the RISC. The balance between miR168 and AGO1 plays important roles in plant development, including flowering and fruiting (Vaucheret et al., 2004; Xian et al., 2014). The top 20 most highly expressed known miRNAs are listed in **Table 1**. The expression levels of novel miRNAs were much lower than the known miRNAs, with the majority (∼80%) barely accumulating several to less than 100 copies. The top 15 most highly expressed novel miRNAs are listed in **Table 2**.

The expression level between the Ogura-CMS and maintainer buds was compared using normalized tag counts. Ten known miRNAs were down-regulated two-fold (p < 10−3 ) in Tyms compared to 231–330 (**Table 3**). Among them, bra-miR157, bramiR158-3p, and bra-miR5718 were expressed at relatively high levels and could play important roles in pollen development. However, bra-miR6030 was the only known miRNA that was up-regulated in the Tyms buds. In addition to these examples, a total of 49 and 27 novel miRNAs were down- and up-regulated in Tyms, accounting for 11.5 and 6.3% of the total 426 novel miRNAs, respectively (Table S8).

To test the accuracy of the RNA-Seq-based expressional profiling, a set of Q-PCR analyses were performed. As shown in **Figure 3**, the relative expression levels were similar between the Q-PCR and RNA-Seq technologies for seven known and

#### TABLE 2 | The most highly expressed novel miRNAs.

one novel miRNAs. These results indicate that the RNA-Seq expression profiles are reliable.

#### Target Prediction and Validation by Degradome Analysis

miRNAs function by regulating their target genes and especially by degrading their target mRNAs in plants. Thus, we performed a degradome sequencing analysis to validate the miRNA

targets. A total of 36,541,512 and 36,595,548 clean tags were generated from the buds of Tyms and 231-330, respectively (Table S9). Among them, 23,253,666 (63.64%) and 24,768,149 (67.68%) tags were mapped to the reference B. rapa genome. A total of 15,154,675 (41.47%) and 17,009,287 (46.48%) tags were mapped to the cDNA sense chains and were used for target prediction (Table S10). A total of 376 mRNAs were identified as miRNA targets, including 136 targeted by 30 known miRNA families and 248 targeted by 100 novel miRNAs. Among which 26 mRNAs were targeted by more than one miRNA (Table S11). A larger proportion known (50.8%) miRNAs compared to novel (23.5%) miRNAs have targets been detected, indicating that the conserved miRNAs have comparatively more valid cutting functions than the newly formed miRNAs. One possible reason for not detecting the targets could be that the target mRNAs were expressed at low levels; thus, the cut ends were not detected by our pipeline. Another possibility is that some miRNAs function by translation inhibition rather than mRNA digestion. Based on GO annotation, the targets were enriched in the "binding," "catalytic activity," and "nucleic acid binding transcription factor activity" terms in the "molecular function" cluster. "Cell," "membrane," and "organelle" were the three most abundant "cellular component" targets. In the "biological process" cluster, the top three terms were "cellular process," "metabolic process," and "response to stimulus" (**Figure 4**).

Most of the known miRNAs targeted transcription factorencoding genes, such as the targeting of SPL, MYB, auxin response factor (ARF), NAC, scarecrow, APETALA 2 (AP2), GROWTH-REGULATING FACTOR (GR), and C3HC4-type RING finger family transcription factors by miR156/157, miR159, miR160, miR164, miR171, miR172, miR396, and miR5716, respectively. This finding indicates that the transcription factors perform primary important roles in bud development. Another large group of miRNAs, including miR158, miR161, miR400, and miR5654, target the pentatricopeptide repeat (PPR) containing protein-encoding genes. PPR-containing proteins are a large family in plants that are mostly located in the mitochondria and chloroplast and play important roles in RNA processing in the two organelles (Fujii and Small, 2011). Some PPR-containing proteins have been identified as fertility restoration (Rf) genes for CMS (Desloire et al., 2003; Wang et al., 2008; Yasumoto et al., 2009). miR1885, miR5719, and miR6030 target disease resistance (CC-NBS-LRR class) genes. miR167 targets IAA-amino acid hydrolase 3 (IAR3), which encodes an auxin conjugate hydrolase that contributes to the jasmonate pathway (Kinoshita et al., 2012; Widemann et al., 2013). MiR393 targets AUXIN signaling F-BOX (AFB) genes, which are auxin receptors that have been widely identified as stress response regulators and control auxin-related development in plants (Si-Ammour et al., 2011; Terrile et al., 2012). These genes and the above mentioned miR160-ARFs indicate that the miRNA-hormone signaling cascades play important roles in flower bud development. Interestingly, three Vacuolar ATP synthase subunit A (VHA-A) tags were identified as targets of miR5712. VHA-A have been demonstrated to be essential for male gametophyte development due to its important role in Golgi organization

(Dettmer et al., 2005). miR5724 targets embryo defective 1473 (emb1473), which encodes a structural constituent of the ribosome. miR158, which is down-regulated in Tyms buds, targets a glutathione S-transferase gene that has been annotated as a restorer-of-fertility (Krishnasamy and Makaroff, 1994; Woo et al., 2008). These miRNAs could play important roles in pollen development.

For the novel miRNAs, 153 out of the 264 targets (58%) were annotated as "catalytic activity," which accounted for the largest fraction of the targets. Only 16 transcription factors were identified as targets of novel miRNAs: four NAC (including three cup-shaped cotyledon) targeted by novel mir410, 446, and 465, four zinc finger family proteins, two homeobox proteins, one MYB, one WRKY, one SPL, one HSFA1E, and one PIL1 (Phytochrome Interacting Factor 3-like 1). Several of the 19 targets possess "transporter activity," including four genes encoding H+-transporting ATPases, two genes encoding Ca2+:H+ antiporters and two genes encoding sucrose transporter SUC1. These genes could play important roles in gametophyte development by transporting protons, ions, mineral elements and energy. A total of 28 targets were annotated as reproductive process-related, including the aforementioned cupshaped cotyledon 2 transcription factors and sucrose transporter SUC1, auxin-responsive protein, and ubiquitin-activating enzyme E1 (Table S12). Among them, an oxidoreductase, a HSFA1E, a C2H2 zinc finger protein 1 and a beta-glucosidase were targets of novel-miR-6, −93, and −321, which were downregulated in Tyms. An F-box/LRR-repeat protein, a heat shock 70 kDa protein 1/8, an auxin-responsive protein IAA and two sucrose transporter SUC1 were targeted by novel-miR-175, −383, −385, and −448, which were up-regulated in Tyms buds.

#### Target Expression Profiles Based on Transcriptome Sequencing

To profile the expression of the targets, a transcriptome sequencing project was performed separately on buds of 231– 330 and Tyms. A total of 25,015,186 and 24,295,576 clean reads were generated from the 231–330 and Tyms buds, respectively. Approximately 76.4 and 72.6% of the reads were mapped to the reference genome and 59.3 and 57.9% of the reads were mapped to the reference gene sets of B. rapa, respectively. The expression of the target genes was called from the transcriptome data for further analysis. A total of 367 and 352 out of the 376 target genes were expressed in the 231–330 and Tyms buds, respectively. Among them, 44 were up-regulated and 32 were down-regulated (two-fold change and FDR < 10−3 ) in Tyms compared with 231–330. Most of the differentially expressed targets were pollen development-related genes according to the annotation (**Tables 4, 5**).

#### miRNA-target Network Underlying Ogura-CMS

The differentially expressed target genes were targeted by 10 known miRNAs (miR159, miR164, miR171, miR172, miR393, miR396, miR397, miR1885, miR5654, and miR6034) and 33 novel miRNAs. However, none of these known miRNAs were significantly differentially expressed according to the two-fold criteria. Based on the combined analysis of the expression profiles of miRNAs and their targets without the significant test filters, the expression patterns of the miRNA-targets could be classified into four clusters. As shown in **Figure 5A**, cluster I contained more than half of the miRNA-target pairs, in which

miRNAs were down-regulated and the targets were up-regulated in Tyms. Approximately 16 and 24% of the miRNA-target pairs were synergistically down-regulated (cluster II) or up-regulated (cluster III), respectively. Less than 8% of the miRNAs were up-regulated while their targets were down-regulated. These results indicate the miRNAs provide an efficient buffer system for fine-tuning the expression of genes involved in bud and pollen development. When a cut-off of >1.5-fold change was used, seven miRNAs were down-regulated and released the expression of 18 targets in Tyms. In contrast, four miRNAs were up-regulated and suppressed the expression of five targets (**Figure 5B**). miR156 and SPLs represent a cascade of controls for many important developmental processes, including the proper development of sporogenic tissues of the anther (Xing et al., 2010). miR156 and SPLs are expressed in high abundance in many tissues, including buds (Zhang et al., 2011b). Our present study found that miR156 and SPLs were down- or up-regulated in Tyms buds by more than 1.5-fold but less than two-fold. It is possible that the expression of these genes was changed more significantly in sporogenic tissues but was submerged by the expression in other tissues of the buds. Using a two-fold threshold, three and two novel miRNAs were up-regulated in 231–330 and Tyms, resulting in the down-regulation of five and three of their target genes in the corresponding tissues, respectively (**Figure 5C**). Among them, a cytochrome c biogenesis protein was targeted by novel-miR-180, a Ca2+/H+ antiporter was targeted by novel-miR-191, a P450 and a glycogenin glucosyltransferase (GGT) were targeted by novel-miR-6; all of these genes have been implicated in pollen development by previous studies (Balk and Leaver, 2001; Welchen and Gonzalez, 2005; Morant et al., 2007; Song et al., 2009; Li et al., 2010, 2012; Rennie et al., 2014) and were upregulated in Tyms in our analysis. Interestingly, novel-miR-448 and novel-miR-335 were specifically expressed in Tyms; this is especially worth noticing because the novel-miR-335 was expressed at a relatively high level. The expression of its

target genes (sucrose transporter SUC1 and H+-ATPase 6) was high in 231–330 but suppressed by approximately 100-fold in Tyms. SUC1 is an enzyme that catalyzes the degradation of sucrose and performs important roles for pollen germination in Arabidopsis (Sivitz et al., 2008). Thus, the down-regulation of SUC1 may result in energy deficiency in Tyms buds and thus abort pollen development. A study in Nicotiana plumbaginifolia reported that co-suppression of an H+-ATPase impaired sucrose translocation and male fertility (Zhao et al., 2000). In our present study, both a SUC1 and an H+-ATPase were indicated as targets of two novel miRNAs and were sufficiently suppressed by these two miRNAs in the Ogura-CMS buds (**Figure 5C**). This finding indicated that the novel-miR-335/H+-ATPase and novel-miR-448/SUC1 cascade could play important roles in male sterility in Ogura-CMS. The secondary structure of the novelmiR-335 and novel-miR-448 and the degradation map of H+- ATPase and SUC1 were shown in **Figure 6**. The suppression of H+-ATPase and SUC1 expression were validated by Q-PCR analysis (**Figure 7**). This finding filled a gap within the crosstalk network between the orf138 locus in mitochondria and the effecter genes in the chromosome. However, more studies are still needed to linkup the orf138 locus and miRNA regulatory networks.

# CONCLUSION

The present study used deep sequencing technology and performed a series of whole genome small RNA, degradome,

areas indicate the dominant mature miRNAs, violet shaded areas indicates the reverse complementary mature miRNAs. (B) Cutting plots of miRNA targets confirmed using degradome sequencing. The corresponding miRNA:mRNA alignments are shown on the top. The red arrows indicate the miRNA-directed cleavage positions. The y-axis shows the nucleotide position in the target gene. The x-axis indicates the number of cleaved ends detected in the degradome analysis.

and transcriptome analyses on Chinese cabbage buds from both Ogura-CMS and its maintainer. A total of 289 (69 families) known and 426 novel miRNAs were identified, which was much higher than the number of miRNAs found in the buds of the closely relative subspecies B. campestris ssp. chinensis (Jiang et al., 2014). These miRNAs not only validated the finding in the

previous study on B. campestris ssp. chinensis but also contained many novel miRNAs that were not previously reported. A large number of targets were firstly validated in this study. The combinational profiling of miRNAs and targets revealed a regulatory network contributing to bud development, especially for pollen engenderation. The finding of two novel miRNA/target cascades (novel-miR-335/H+-ATPase and novel-miR-448/SUC1) will provide new insights into the communication between the mitochondria and chromosome and take one step toward filling in the gap in the regulatory network mechanism from the orf138 locus to the pollen abortion in the Ogura-CMS plants. The true functions of these two miRNAs warrant more solid experiments, such as the transgenic complementary test.

#### REFERENCES

#### AUTHOR CONTRIBUTIONS

XW and XWZ designed the study. XW performed the experiments. XW and XHZ analyzed the data and drafted the manuscript. YY and QZ assisted with the bioinformatics analysis and aided in writing the manuscript. QY, ZW, YZ, and WJ aided in performing the experiments. XL and FW modified the manuscript. All of the authors carefully checked and approved this version of the manuscript.

#### DATA ACCESS

RNAseq are submitted to EMBL/NCBI/SRA with the accession numbers SRR2132359, SRR2136647, SRR2149955, SRR2132463, SRR2136646, SRR2149956.

#### FUNDING

This work was supported by funding from North Henan Station, the National Vegetable Industry Technology System (CARS-25-27), Excellent technology innovation team of Henan province, the National High Technology Research and Development Program of China (863 Program) (2012AA100202- 7), The National Key Technology R&D Program of China (2012BAD02B01-3).

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2015. 00894

**Conflict of Interest Statement:** 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.

Copyright © 2015 Wei, Zhang, Yao, Yuan, Li, Wei, Zhao, Zhang, Wang, Jiang and Zhang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Cytogenetic Diversity of Simple Sequences Repeats in Morphotypes of Brassica rapa ssp. chinensis

Jin-shuang Zheng1,2, Cheng-zhen Sun1,2, Shu-ning Zhang1 , Xi-lin Hou1 \* and Guusje Bonnema3

1 State Key Laboratory of Crop Genetics and Germplasm Enhancement, Key Laboratory of Biology and Germplasm Enhancement of Horticultural Crops in East China, Ministry of Horticulture, Nanjing Agricultural University, Nanjing, China, 2 Hebei Normal University of Science and Technology, Qinhuangdao, China, 3 Wageningen UR Plant Breeding, Wageningen University and Research Centre, Wageningen, Netherlands

A significant fraction of the nuclear DNA of all eukaryotes is comprised of simple sequence repeats (SSRs). Although these sequences are widely used for studying genetic variation, linkage mapping and evolution, little attention had been paid to the chromosomal distribution and cytogenetic diversity of these sequences. In this paper, we report the distribution characterization of mono-, di-, and tri-nucleotide SSRs in Brassica rapa ssp. chinensis. Fluorescence in situ hybridization was used to characterize the cytogenetic diversity of SSRs among morphotypes of B. rapa ssp. chinensis. The proportion of different SSR motifs varied among morphotypes of B. rapa ssp. chinensis, with tri-nucleotide SSRs being more prevalent in the genome of B. rapa ssp. chinensis. We determined the chromosomal locations of mono-, di-, and tri-nucleotide repeat loci. The results showed that the chromosomal distribution of SSRs in the different morphotypes is non-random and motif-dependent, and allowed us to characterize the relative variability in terms of SSR numbers and similar chromosomal distributions in centromeric/peri-centromeric heterochromatin. The differences between SSR repeats with respect to abundance and distribution indicate that SSRs are a driving force in the genomic evolution of B. rapa species. Our results provide a comprehensive view of the SSR sequence distribution and evolution for comparison among morphotypes B. rapa ssp. chinensis.

Keywords: Brassica rapa, simple sequence repeats, fluorescence in situ hybridization, cytogenetic diversity, heterochromatin

# INTRODUCTION

Simple sequence repeats (SSRs), also known as microsatellites, are composed of 1–6 nucleotide motifs that are repeated in tandem and are widely and non-randomly distributed in 100–1000s of copies in the genomes of both monocots and dicots (Tautz and Renz, 1984; Toth et al., 2000; Mortimer et al., 2005; Lawson and Zhang, 2006; Hong et al., 2007). Microsatellites are found predominantly in heterochromatin regions, such as centromeric, peri-centromeric, and sub-distal regions of eukaryotic chromosomes (Yang et al., 2005; Li et al., 2010), and sex chromosomes in animals (Beckmann and Weber, 1992; Cuadrado and Jouve, 2011), usually associated with

#### Edited by:

Sarvajeet Singh Gill, Maharshi Dayanand University, India

#### Reviewed by:

Maoteng Li, Huazhong University of Science and Technology, China Santosh Tiwari, Maharshi Dayanand University, India

> \*Correspondence: Xi-lin Hou<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Crop Science and Horticulture, a section of the journal Frontiers in Plant Science

Received: 27 February 2016 Accepted: 04 July 2016 Published: 26 July 2016

#### Citation:

Zheng J-s, Sun C-z, Zhang S-n, Hou X-l and Bonnema G (2016) Cytogenetic Diversity of Simple Sequences Repeats in Morphotypes of Brassica rapa ssp. chinensis. Front. Plant Sci. 7:1049. doi: 10.3389/fpls.2016.01049

**202**

the constitutive heterochromatin (Lohe et al., 1993; Pedersen et al., 1996). Variously non-random, taxon-specific patterns of SSR occurrence call for functional interpretations. Morgante et al. (2002) hypothesized that the relative frequency of microsatellites is higher in the single- or low-copy regions of the genome than in the repetitive regions. Mortimer et al. (2005) showed that SSRs are associated with and around transcribed sequences in Arabidopsis. Li et al. (2004) reviewed that SSRs in different positions in a gene can play important roles in regulating its expression and determining the function of its products. The accumulated evidence indicates that SSRs play an important role in chromatin organization, regulation of gene activity (Nagaki et al., 2004), recombination, DNA replication, the cell cycle, the mismatch DNA repair system (Li et al., 2002), and protein coding regions (Sonah et al., 2011).

The fluorescence in situ hybridization (FISH) technique has been used to localize one or more SSR loci on chromosomes. FISH has become a strategy for chromosome diagnosis and for investigating plant genome organization (Cuadrado and Schwarzacher, 1998; Cuadrado and Jouve, 2002). SSRs change rapidly during evolution, and thus display polymorphism at homologous sites between closely related species. Begum et al. (2009) showed that the rapid evolution of repetitive DNA sequences has resulted in species-specific repeat variants and the generation of novel repeat families. This characteristic has made SSRs useful as markers in comparative diversity analysis (Zhebentyayeva et al., 2003; Yu et al., 2010; Fang et al., 2013; Zhu et al., 2013) and genetics research (Nanda et al., 1991). SSRs are highly abundant within genomes, they can be widely dispersed or be confined to certain chromosomal regions, and they display a high degree of length polymorphism (Katti et al., 2001; Lawson and Zhang, 2006). Carmona et al. (2013) defined cytogenetic diversity in terms of the differences in abundance and distribution of microsatellites, and also found some specific and motif-dependent hybridization patterns. The repeats AG, AAG, ACT, and ATC presented different in situ hybridization patterns that provided cytogenetic landmarks for chromosome identification in barley, Hordeum vulgare ssp. vulgare (Carmona et al., 2013). Altogether, such variation could be used to determine evolutionary relationships between related species.

Brassica rapa belongs to the A genome species group in the Brassicaceae with 2n = 20 chromosomes (Nagaharu, 1935), which had a monophyletic origin (Lysak et al., 2005; Cheng et al., 2013). B. rapa comprises several sub-species, such as non-heading Chinese cabbage (B. rapa ssp.chinensis), Chinese cabbage (B. rapa ssp. pekinensis; Koo et al., 2004), and turnip (B. rapa L. ssp. Rapifera; Snowdon, 2007). B. rapa ssp. chinensis is one of the most important leafy vegetable forms in B. rapa, and consists of five morphotypes (Pak-choi, Wu ta cai, Cai xin, Fen nie cai, and Tai cai; Viehoever et al., 19201 ; Gladis and Hammer, 1992; Zheng et al., 2015). The phylogenetic relationships between some B. rapa ssp. chinensis morphotypes have been determined from morphological, ecological, and molecular data (Yu et al., 2010). In this paper, systematic research was performed to investigate the cytogenetic diversity between intra-specific forms of B. rapa.

The objective of the work presented here is to characterize the cytogenetic diversity of SSRs between morphotypes of B. rapa ssp. chinensis that represent a broad range of cytogenetic diversity. The available genome sequence of Chinese cabbage is an important and fundamental resource for understanding this species2 (Wang et al., 2011), but the vast majority of heterochromatic regions remain essentially uncharacterized. The current estimates of SSR frequencies in many organisms differs from reality after comparisons with sequence databases (Hong et al., 2007; Cavagnaro et al., 2010; Gao et al., 2011). The distribution of SSRs in databases has been reported for the Chinese cabbage genome (Hong et al., 2007). In view of their ubiquity and functional importance, detailed information will be necessary to explore the comparative cytogenetics of SSRs in B. rapa.

Simple sequence repeats appear to be more abundant in noncoding regions than in coding regions of plant genomes (Hong et al., 2007; Cavagnaro et al., 2010). Tri-nucleotide repeats, the most abundant SSR types in many species (Gao et al., 2003; Shi et al., 2013), were found to be the most frequent in protein coding regions (Sonah et al., 2011). In addition to tri-nucleotide repeats, mono- and di-nucleotide repeats were also predominant in the B. rapa genome (Hong et al., 2007; Gao et al., 2011). In this study, we selected mono-, di-, and tri-nucleotide repeats for physical mapping on the chromosomes of B. rapa. FISH was performed with mono-, di-, and tri-nucleotides to detect the distributional profile of SSRs and to enhance our understanding of genome organization. The distributional characterization of SSRs revealed a range of cytogenetic diversity that could relate to genome organization, function, and evolutionary trends. We demonstrated that: (1) not all of the SSR-based probes produced FISH signals on all B. rapa ssp. chinensis chromosomes; (2) some SSR signal intensity did not show a relationship to the abundance in the genome database; (3) the distributional patterns of SSR signals depended on the SSR motif used and the species analyzed; and (4) differences in SSR abundance and density were shown within and between genomes.

#### MATERIALS AND METHODS

#### Plant Materials