Comparative transcriptome analysis of the effect of different heat shock periods on the unfertilized ovule in watermelon (Citrullus lanatus)

Journal of Integrative Agriculture 2020, 19(2): 528–540 Available online at www.sciencedirect.com

ScienceDirect

RESEARCH ARTICLE

Comparative transcriptome analysis of the effect of different heat shock periods on the unfertilized ovule in watermelon (Citrullus lanatus) ZHU Ying-chun1, 2, SUN De-xi2, DENG Yun2, AN Guo-lin2, LI Wei-hua2, SI Wen-jing2, LIU Jun-pu2, SUN Xiao-wu1 1

College of Horticulture, Hunan Agricultural University, Changsha 410128, P.R.China

2

Zhengzhou Fruit Research Institute, Chinese Academy of Agricultural Sciences, Zhengzhou 450009, P.R.China

Abstract In vitro gynogenesis is an important tool used in haploid or homozygous double-haploid plant breeding. However, because of low repeatability, embryoid induction rate and quality, the molecular mechanisms remain poorly understood. Heat shock treatment can promote the transformation of the gametophytic pathway into the sporophyte pathway, which induces the occurrence of haploid. In this study, unfertilized ovaries were heat shocked for 0 h (A0) before flowering and for 0 h (A1), 4 h (A3), 8 h (A5), 12 h (A7), and 24 h (A8), respectively, at 37°C at the first day of the flowering stage. The ovule enlargement rate was increased from 0% at 25°C to 96.8% at 37°C (24 h treatment). Thus, we aimed to investigate the gene expression patterns in unfertilized ovules of watermelon after different periods of heat shock by using RNA-Seq technology. The results showed that compared with A3, A5, A7, and A8, the biosynthesis of amino acid, glycine, serine and threonine metabolic pathways in A1 has changed significantly. This indicated that heat shock treatment affected the synthesis and transformation of amino acids during ovule expansion. The transcriptome data suggested gene expressions of ovule growth were significantly changed by heat-specific influences. The results provide new information on the complex relationship between in vitro gynogenesis and temperature. This provides a basis for further study of the mechanism of heat shock affecting the expansion of watermelon ovule. Keywords: watermelon, heat shock, unfertilized ovule, ovule enlargement, transcription factors

(Citrullus lanatus Thunb.) is widely grown in the world

1. Introduction As an economically important fruit crop, watermelon

(Vinoth and Ravindhran 2016) and the yield was 114 million tons in 2017 (Yang et al. 2019). Watermelon contains a variety of nutrients and is a well-known component of the daily nutrition for the world’s population (Zhu et al. 2017). Although distinct advantages may be obtained from cross breeding, continuous self-fertilization by manual selection,

Received 11 February, 2019 Accepted 10 July, 2019 ZHU Ying-chun, E-mail: [email protected]; Correspondence LIU Jun-pu, Tel/Fax: +86-371-60247661, E-mail: liujunpu@caas. cn; SUN Xiao-wu, E-mail: [email protected] © 2020 CAAS. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http:// creativecommons.org/licenses/by-nc-nd/4.0/). doi: 10.1016/S2095-3119(19)62777-2

which is a long and tedious process, is necessary to obtain genetically stable and pure inbred lines. In contrast, haploidy culture shortens the duration of plant breeding and more easily provides homozygous pure lines. San (1979) was the first to use unfertilized ovaries for breeding when he achieved the haploid of Hordeum vulgare in a process termed in vitro gynogenesis. Since then, in vitro

ZHU Ying-chun et al. Journal of Integrative Agriculture 2020, 19(2): 528–540

gynogenesis has become an important method for obtaining haploid or homozygous double-haploid plant production in plant breeding (Shalaby 2007), especially for Cucurbitaceae plants. For instance, haploid and diploid breeding of Cucurbita pepo was obtained by in vitro culture of unfertilized ovules (Chamboner and Dumas de Vaulx 1985). Haploid plants from C. pepo (cv. Eskandarani) may also be produced by unpollinated ovules in ovaries of squash plants picked one day before anthesis and exposed to cold temperature and 2,4-dichlorophenoxy acetic acid (2,4-D) solid MS medium (Metwally et al. 1998a, b). Subsequently, the haploid of summer squash was successfully induced by ovules using in vitro gynogenesis. This work showed that genotype, temperature and sucrose concentration had different effects on the induction of squash (Shalaby 2007). In addition, both homozygous plants of cucumber and melon were produced by in vitro gynogenesis induction (Gémesjuhász et al. 2002; Lotfi et al. 2003; Malik et al. 2011). Although in vitro gynogenesis has been receiving more attention as an important tool for breeding and genetic research of Cucurbitaceae, various factors (e.g., material status, genotype and temperature pretreating) have resulted in low repeatability, and low embryo development and quality (Zou et al. 2017). Thus, related methodologies are in urgent need of exploration. Inductive heat shock treatment is an important method for increasing the efficiency of gynogenesis in the developmental stage of the female gametophyte for Cucurbitaceae, such as C. sativus (Gémesjuhász et al. 2002), C. pepo (Metwally et al. 1998b; Shalaby 2007), C. moschata (Kwack and Fujieda 1988). Heat shock treatment induces the development of the gametophytic pathway to the sporophytic pathway, and leads to haploid production (Shalaby 2007). Although the coordinated expression of genes is the major mechanism of in vitro gynogenesis in plants, the details about heat shock stress and haploid culture in watermelon are still lacking. Since a heat shock study on the ovule enlargement rate of unfertilized ovaries showed that heat shock treatment could affect the ovule enlargement rate in watermelon, we selected cultivated watermelon Zhongke 6 for a comparative transcriptome sequencing (RNA-Seq) study. In this study, we investigated the effect of heat shock stress on unfertilized ovaries from the cultivated watermelon by using transcriptome profiling. The differential expression genes (DEGs) in each of six tissues from the unfertilized ovaries of watermelon in response to heat treatment were identified. The acquired data provide insight into a series of signal pathways and regulators involved in embryogenesis induction in response to heat shock applied for different time periods.

529

2. Materials and methods 2.1. Plant treatments Plants of watermelon (Zhongke 6) is grown in Xinxiang Base of Chinese Academy of Agricultural Sciences, Xinxiang City, Henan Province. Unfertilized ovaries were collected by wiping off the tomentum and washed with running water for 30 min. Subsequently, the samples were fixed with 75% ethyl for 30 s, transversely cut to 1 mm sheet, placed in 8% sodium hypochlorite for 20 min, and washed with sterile water three times. Afterwards, ovary sections were seeded on MS medium containing, 0.02 mg L–1 TDZ (N-phenyl-N´-1,2,3thiadiazol-5-ylurea), 0.25 mg L–1 NAA (a-naphthaleneacetic acid), 0.50 mg L –1 6-BA (6-benzylaminopurine), 30 g saccharose and 6 g agar and cultured in the dark for 1 day at 37°C. RNA was extracted from unfertilized ovaries heat shocked for 0 h (A0) before flowering and for 0 h (A1), 4 h (A3), 8 h (A5), 12 h (A7), and 24 h (A8), respectively, at 37°C at the first day of the flowering stage by Trizol (TaKaRa, Japan) according to the manufacturer’s protocols. The samples were wrapped with silver paper and treated with liquid nitrogen immediately and stored at –80°C for use.

2.2. Bioinformatics analysis Comparative transcriptome analysis was performed on extracted RNA as described by Du et al. (2016) and Song et al. (2017). Total RNA of each sample was isolated using a Quick RNA Isolation Kit (Bioteke Corporation, Beijing, China) and then characterized on a 1% agarose gel and examined with a NanoDrop ND1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA). The RIN (RNA integrity number) values (>8.0) of these samples were assessed using an Agilent 2100 Bioanalyzer (Santa Clara, CA, USA). The construction of the libraries and the RNA-Seq were performed by the Anoroad Corporation (Beijing, China). mRNA was enriched and purified with oligo(dT)-rich magnetic beads and then broken into short fragments. Taking these cleaved mRNA fragments as templates, first- and secondstrand cDNAs were synthesized. The resulting cDNAs were then subjected to end-repair and phosphorylation using T4 DNA polymerase and Klenow DNA polymerase. After that, an ‘A’ base was inserted as an overhang at the 3´ ends of the repaired cDNA fragments and Illumina paired-end solexa adaptors were subsequently ligated to these cDNA fragments to distinguish the different sequencing samples. To select a size range of templates for downstream enrichment, the products of the ligation reaction were purified and selected on a 2% agarose gel. Next, PCR amplification was performed to enrich the purified cDNA template. Finally, the four libraries

2.4. Transcript analysis by quantitative realtime PCR (qRT-PCR)

To verify the reliability of the RNA-Seq analyses, 11 candidate genes were selected for qRT-PCR analysis. A total of 500 ng μL–1 RNA were extracted with Trizol (TaKaRa, Japan) and reverse transcribed to cDNA with a PrimeScript Reverse Transcriptase (RT) Reagent Kit (Perfect Real Time) (TaKaRa, Japan). The templates were performed with CFX96 Real-Time PCR Detection System (Bio-Rad, USA) by using SYBR ® Green Kit (TaKaRa) with the following parameters: 95°C for 15 s, 60°C for 30 s and finally 72°C for 20 s. The primers used are listed in Table 1. The relative DEGs expression was

1)

The GO functional classification for the DEGs with functional significance was analyzed by Gene Ontology Annotation Plot Software (http://www. geneontology.org/). GO terms in each comparison group with corrected P-value≤0.05 were regarded as significant enrichment. In order to identify the pathways that the DEGs participated in, the DEGs were enriched in the KEGG database. Sense (5´→3´) CCCTAATTACTCCGGGTGGG CTGAAAAGGTGCCATTCCGAC GACCGGATGCTTTGACGAAC CTAGCCAGTAGCTTGGAGCC CGGAAGCTTTTCTGGGGCAA AATCGCGGAAGCTTTTCTGG CTACCCTGAACTTCTTAAAGCC CATCGGAGAATTCATCTGGGGTT TGACGCTACCAATCCAAAGC TCGACTTTGACTCCTCTGCG TCGCCAAAAGCTGAGAAGGA CCCTTTGCTGCCCAAGTTTC GGACGAATTCTCCTGGACGG AGAACCCAGAGAAACAATGGC GTTTCCGCCGGTTGATTGTC AATGCGGCTTCACAACTTGG AAGCGAACGATCTCGCTACT TAACCCGCCACGAACTACG GGTTCCTGCAGACGAGGATG AACAGAAGAGAGTTCGCCCC

Anti-sense (5´→3´) CTTCGGAATCACCAGCCTGT GGGTCGTGAACAAAAGCCTC TCTTCTGGAGACGACGTTCTG GGCACTTGCACGATACACAC ACGAGACGCCTTTTTGTCCA CCTGCATCTTTCGTCGCTTG CCCATGGAACATCACCCACA CACCCAACAACCTGAGCCTT GAAGCCCCATCTGTAGTGGTT CCGTTCCTCATCCTCGTCTG TACAAGCTGATGGTGCTGGG TGACAGGACAGAGGTGGACA GCGCCAGTGATTGTGATGTC ATCACCATTGTCAGCCCCTA CAGCCTGACCAACAGGGAG AGCTCCTGACCTTTCTGGGA ACGACGTCCTCAACATTCGG TTGGTAAGAGCGTTGGGACC TGTCTTTGGCTTCTCAGTGTC ATTCAGGTTGGGCTTCGGAC

NR: Seq-id, optimal alignment results of genes in non-redundant protein sequences database; NT: Seq-id, optimal alignment results of genes in nucleotide collection database.

NT: Seq_id gi|778730759|ref|XM_004151012.2| gi|659110120|ref|XM_008456837.1| gi|733577844|emb|LN713260.1| gi|778698981|ref|XM_011656331.1| gi|659069538|ref|XM_008452314.1| gi|659119703|ref|XM_008461576.1| gi|659096791|ref|XM_008451067.1| gi|778705302|ref|XM_004135152.2| gi|659102710|ref|XM_008454052.1| gi|778660005|ref|XM_004150158.2| gi|659127727|ref|XM_008465635.1| gi|659110994|ref|XM_008457299.1| gi|659123959|ref|XM_008463702.1| gi|778669307|ref|XM_011650929.1| gi|778682724|ref|XM_011653461.1| gi|778717308|ref|XM_004143321.2| gi|733578671|emb|LN713266.1| gi|733577844|emb|LN713260.1| gi|659088134|ref|XM_008446599.1| gi|659129399|ref|XM_008466444.1|

2.3. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses

NR: Seq_id gi|449466693|ref|XP_004151060.1| gi|778724794|ref|XP_011658863.1| gi|449463464|ref|XP_004149454.1| gi|778698982|ref|XP_011654633.1| gi|659069539|ref|XP_008450536.1| gi|659119704|ref|XP_008459798.1| gi|449465437|ref|XP_004150434.1| gi|449434833|ref|XP_004135200.1| gi|659102711|ref|XP_008452274.1| gi|659086500|ref|XP_008443966.1| gi|659127728|ref|XP_008463857.1| gi|700188633|gb|KGN43866.1| gi|659123960|ref|XP_008461924.1| gi|449462248|ref|XP_004148853.1| gi|659123761|ref|XP_008461826.1| gi|449451239|ref|XP_004143369.1| gi|659068378|ref|XP_008444113.1| gi|307136351|gb|ADN34165.1| gi|659088135|ref|XP_008444821.1| gi|659129400|ref|XP_008464666.1|

were sequenced using an Illumina HiSeq™ 2000. After removing those reads with only adaptor and unknown nucleotides >5%, or those that were of low quality, the clean reads were filtered from the raw reads. The clean reads were considered for further analyses. The raw RNA-Seq reads were trimmed from the adapter sequences and mapped against the watermelon genome reference (http:// cucurbitgenomics.org/organism/1), to determine the expression levels of the identified genes. According to the data of uniquely mapped reads, FPKM (fragments per kilobase per million mapped reads) values were calculated with Cufflinks software (2.2.1) (Trapnell et al. 2010). The sequence alignments were performed by TopHat v2.0.12 (Trapnell et al. 2009), and DEGs between two sample pairs were analyzed using the DESeq R package, which were carried out and identified by applying an adjusted P-value, in which false discovery rate (FDR, Bonferroni corrected) was taken into account in multiple tests. DEGs were defined with the absolute value of Log2FoldChange>1 (FoldChange>2 or <0.5) and adjusted P-value<0.05 (Wang et al. 2010).

Gene name Cla000855 Cla001617 Cla004223 Cla004923 Cla005777 Cla006725 Cla007544 Cla010153 Cla011760 Cla012214 Cla013275 Cla016077 Cla016401 Cla016552 Cla017291 Cla018349 Cla019599 Cla021189 Cla022413 Cla023261

Table 1 Primer sequences used for qRT-PCR1)

530 ZHU Ying-chun et al. Journal of Integrative Agriculture 2020, 19(2): 528–540

531

ZHU Ying-chun et al. Journal of Integrative Agriculture 2020, 19(2): 528–540

(Fig.  2-B). In total, we generated 726.6 million clean reads, and of which 91.0% on average were mapped to the reference genome (Table 3).

calculated by using the 2–ΔΔCt method (Livak and Schmittgen 2001).

3. Results

3.2. DEGs from ovaries under heat stress

3.1. Mapping transcriptome sequencing

The DEGs identified from each paired group A1 vs. A3 (A3_A1), A1 vs. A5 (A5_A1), A1 vs. A7 (A7_A1) and A1 vs. A8 (A8_A1) were counted in a Venn diagram (Fig. 3-A). The results revealed that 2 363 genes in the A1_A0 pair, 6 332 genes in the A3_A1 pair, 7 235 genes in the A5_A1 pair, 7 444 genes in the A7_A1 pair, and 8 219 genes in the A8_A1 pair were differentially expressed, whereas only 1 127, 91 and 868 genes were differentially expressed in the A5_A3, A7_A5 and A8_A7 paired groups, respectively (Fig.  3-B). Remarkably, more than 3 000 DEGs were evident in groups A3, A5, A7 and A8 compared with group A1 (Fig. 3-C). The number of DEGs was greater at the A8

The effect of culturing at 25 and 37°C on the ovule enlargement rate is shown in Fig.  1. At 25°C, the enlargement rate was 0% for all time periods. At 37°C, however, the enlargement rate was 96.8% in 24 h, showing a positive effect of heat stress on ovule enlargement. The libraries and Illumina-based sequenced reads were performed on the HiSeq™ 2500 platform with 6 sets of unfertilized watermelon ovaries. The output summary of all samples is listed in Table 2. The differences in expression of each sample were evaluated by the overall distribution (Fig.  2-A) and overall dispersion of expression quantity

25°C

37°C

B

0h

0h

4h

4h

8h

8h

12 h

12 h

24 h

24 h

Enlargement rate (%)

A

Zhongke 6

150 100 50 0

25

37

Exposure temperature (°C)

Fig. 1 Ovule enlargement at 25 and 37°C, respectively. A, the phenotype of ovule enlargement. B, the ovule enlargement rate. Table 2 Summary of the Illumina HiSeq™ 2500 sequencing output for all samples Sample1) A0-1 A0-2 A0-3 A1-1 A1-2 A1-3 A3-1 A3-2 A3-3 A5-1 A5-2 A5-3 A7-1 A7-2 A7-3 A8-1 A8-2 A8-3 1)

Total raw reads 47 081 634 49 044 448 42 139 890 48 017 454 48 507 182 50 387 124 46 992 236 41 909 192 50 339 082 47 972 580 46 802 272 43 866 962 47 203 956 48 087 730 47 353 418 42 736 910 46 004 012 49 467 498

Total raw bases 7 062 245 100 7 356 667 200 6 320 983 500 7 202 618 100 7 276 077 300 7 558 068 600 7 048 835 400 6 286 378 800 7 550 862 300 7 195 887 000 7 020 340 800 6 580 044 300 7 080 593 400 7 213 159 500 7 103 012 700 6 410 536 500 6 900 601 800 7 420 124 700

Total clean reads 45 031 430 46 570 226 39 195 412 45 627 734 44 681 186 47 530 602 44 254 722 39 568 706 47 996 666 46 082 170 44 731 414 40 394 200 45 219 312 45 549 304 45 215 288 40 361 216 43 902 902 46 590 014

Total clean bases 6 754 714 500 6 985 533 900 5 879 311 800 6 844 160 100 6 702 177 900 7 129 590 300 6 638 208 300 5 935 305 900 7 199 499 900 6 912 325 500 6 709 712 100 6 059 130 000 6 782 896 800 6 832 395 600 6 782 293 200 6 054 182 400 6 585 435 300 6 988 502 100

Q30 (%)2) 92.95 93.45 93.25 92.20 93.07 93.17 91.79 93.12 93.18 92.36 92.02 90.10 92.95 89.76 92.86 92.95 89.71 89.23

Unfertilized ovaries heat shocked for 0 h (A0) before flowering and for 0 h (A1), 4 h (A3), 8 h (A5), 12 h (A7), and 24 h (A8), respectively, at 37°C at the first day of the flowering stage. Each treatment was repeated for 3 times. 2) Percentage of bases with Phred greater than 30 in total.

532

ZHU Ying-chun et al. Journal of Integrative Agriculture 2020, 19(2): 528–540

Distrubution of sample expression A0-1 A0-2 A0-3 A1-1 A1-2 A1-3 A3-1 A3-2 A3-3

0.20

0.10

0.05

A5-1 A5-2 A5-3 A7-1 A7-2 A7-3 A8-1 A8-2 A8-3

120 100 80 60 40 20

0.00

0 −5

0

5 log2FPKM

10

15

A0-1 A0-2 A0-3 A1-1 A1-2 A1-3 A3-1 A3-2 A3-3 A5-1 A5-2 A5-3 A7-1 A7-2 A7-3 A8-1 A8-2 A8-3

Density

0.15

FPKM distribution

B

FPKM value

A

Sample

Fig. 2 The gene expression quantity (A) and boxplot of FPKM (fragments per kilobase per million mapped reads) (B) under different experimental conditions. A0, unfertilized ovaries heat shocked for 0 h before flowering; A1, A3, A5, A7, and A8, unfertilized ovaries heat shocked for 0, 4, 8, 12, and 24 h, respectively, at 37°C at the first day of the flowering stage. Each treatment was repeated for 3 times. Table 3 Summary of the transcriptome data Sample1) A0-1 A0-2 A0-3 A1-1 A1-2 A1-3 A3-1 A3-2 A3-3 A5-1 A5-2 A5-3 A7-1 A7-2 A7-3 A8-1 A8-2 A8-3 1)

Total reads 45 031 430 46 570 226 39 195 412 45 627 734 44 681 186 47 530 602 44 254 722 39 568 706 47 996 666 46 082 170 44 731 414 40 394 200 45 219 312 45 549 304 45 215 288 40 361 216 43 902 902 46 590 014

Mapped reads 41 582 535 43 038 630 36 256 357 41 390 412 41 213 258 43 930 911 40 016 028 36 505 275 44 364 079 42 209 317 40 561 913 35 666 544 41 680 160 40 095 650 41 600 239 37 023 366 38 608 309 40 818 312

Mapping rate (%) 92.34 92.42 92.50 90.71 92.24 92.43 90.42 92.26 92.43 91.60 90.68 88.30 92.17 88.03 92.00 91.73 87.94 87.61

Exon 15 367 764 (62.60%) 15 199 156 (60.66%) 12 918 972 (60.60%) 16 026 077 (63.57%) 15 534 913 (62.67%) 17 059 879 (63.99%) 14 641 023 (62.82%) 13 283 083 (62.24%) 16 208 673 (62.70%) 15 370 360 (62.74%) 14 821 726 (62.56%) 13 189 759 (63.12%) 15 079 077 (62.70%) 14 579 454 (62.57%) 15 370 502 (62.55%) 13 519 881 (63.50%) 14 402 835 (64.43%) 15 302 297 (64.67%)

Intron 1 702 959 (6.94%) 2 032 361 (8.11%) 1 433 776 (6.73%) 1 606 493 (6.37%) 1 630 798 (6.58%) 1 686 790 (6.33%) 1 589 220 (6.82%) 1 542 656 (7.23%) 1 830 259 (7.08%) 1 689 983 (6.90%) 1 766 038 (7.45%) 1 482 953 (7.10%) 1 624 692 (6.76%) 1 659 195 (7.12%) 1 817 588 (7.40%) 1 286 796 (6.04%) 1 398 032 (6.25%) 1 527 403 (6.46%)

Intergenic 7 479 662 (30.47%) 7 823 656 (31.23%) 6 964 948 (32.67%) 7 576 995 (30.06%) 7 620 797 (30.75%) 7 912 706 (29.68%) 7 075 292 (30.36%) 6 514 694 (30.53%) 7 812 870 (30.22%) 7 439 035 (30.36%) 7 104 896 (29.99%) 6 223 418 (29.78%) 7 346 894 (30.55%) 7 063 450 (30.31%) 7 383 251 (30.05%) 6 485 716 (30.46%) 6 551 954 (29.31%) 6 831 211 (28.87%)

Unfertilized ovaries heat shocked for 0 h (A0) before flowering and for 0 h (A1), 4 h (A3), 8 h (A5), 12 h (A7), and 24 h (A8), respectively, at 37°C at the first day of the flowering stage. Each treatment was repeated for 3 times.

exposure period than at the other periods, indicating that the longer heat time had an influence on the relative gene expression of ovaries during embryo production. DEGs that were identified in different treatments are presented in a heat-map graph (Fig. 4-A). The top 5% of the DEGs were clustered in groups A1_A0, A3_A1, A5_A1, A7_A1, and A8_A1 (Fig. 4-B–F, Appendices A–E).

3.3. Functional annotation of the DEGs To annotate the functions of the DEGs involved in heat

treatment at different times, all DEGs were mapped to terms in the GO and KEGG analyses (Figs. 5 and 6). There were a total of 30 930 DEGs in the A3_A1 group, 35 490 in the A5_A1 group, 36 751 in the A7_A1 group and 40 601 in the A8_A1 group annotated in GO databases. In the group A8_A1 with the most DEGs (Fig. 5-E), 9 090 counts of upregulation and 8 399 counts of down-regulation were involved in biological processes, 8 029 counts of up-regulation and 6 730 counts of down-regulation were involved in cellular component, and 4 459 counts of up-regulation and 3 894 counts of down-regulation were involved in molecular

533

ZHU Ying-chun et al. Journal of Integrative Agriculture 2020, 19(2): 528–540

A A3_A1

A7_A1

A8_A1

167

1 225 439

87 632

A5_A1

170

139

1 348

219

4 645 169

84

193

426 150

150 −lgP 100

50 40 −lgP 20 30 10

1 362 1 064

_A

5

566 302

A8

3

_A A8

1

_A A8

_A A8

3

_A A7

1

_A A7

3

_A A7

1

_A A5

1

_A A5

_A

0

0

5

53 38

7

639

472

382

_A

−5

0 5 Log2FoldChange

3 812

886

745

A3

−5 0 5 10 Log2FoldChange

3 477

1 416

A1

−10

4 407 3 967

2 000 947

6

Down Up None

10 0 5 Log2FoldChange

2 104

1 000

−2 0 2 4 Log2FoldChange

0

Down 3 818

−lgP 100 150 200 250

10 12 14 −lgP 6 8 2 0

200 −5

3 417 3 311 3 021

3 000

−4

−lgP 100 150 0

0

0 5 Log2FoldChange

4 000

A8_A1

A8_A7

50

50

−lgP 100

150

200

8

A8_A5

Up

−4 −2 0 2 4 6 Log2FoldChange

4

40 20 −6 −4 −2 0 2 4 6 Log2FoldChange

A8_A3

5 000

0 5 10 Log2FoldChange A7_A5

0 −5 0 5 10 Log2FoldChange

−5

−5

A7_A3

−lgP 60 80 100 120

80 100 −lgP 60 40 20 0 −10

Number of DEGs

0 5 10 Log2FoldChange

0

50

0

−5

A7_A1

A5_A3

0

20 −lgP 10 15 5 0

−5 0 5 Log2FoldChange

C

A5_A1 −lgP 10 20 30 40 50 60

A3_A1 −lgP 20 40 60 80 100 120

A1_A0

0

B

Fig. 3 Venn diagram from four comparison groups (A), volcano plots from 11 comparison groups (B), and genes distribution from 11 comparison groups (C) of differentially expressed genes (DEGs) in watermelon ovaries under heat stress. Significantly DEGs were treated with red dots (upregulated) or green dots (down-regulated), others indicated with gray dots in Fig. 3-B. A0, unfertilized ovaries heat shocked for 0 h before flowering; A1, A3, A5, A7, and A8, unfertilized ovaries heat shocked for 0, 4, 8, 12, and 24 h, respectively, at 37°C at the first day of the flowering stage.

534

ZHU Ying-chun et al. Journal of Integrative Agriculture 2020, 19(2): 528–540

A Color key

−2 2 Value

6

A1-3 A1-1 A1-2 A0-1 A0-2 A0-3 A3-1 A3-2 A3-3 A8-3 A8-2 A8-1 A7-1 A7-2 A7-3 A5-2 A5-1 A5-3

−6

C

D

E

F

2

B

A1_A0

A3_A1

A5_A1

A7_A1

A8_A1

Fig. 4 Heat map from 18 samples (A) and top 5% of the differentially expressed genes (DEGs) in five comparison groups: A1_A0 (B), A3_A1 (C), A5_A1 (D), A7_A1 (E), and A8_A1 (F), respecitvely. The bar color relects the gene expression levels. A0, unfertilized ovaries heat shocked for 0 h before flowering; A1, A3, A5, A7, and A8, unfertilized ovaries heat shocked for 0, 4, 8, 12, and 24 h, respectively, at 37°C at the first day of the flowering stage.

lar co

mpone

cellu lar single metab oc -o an olic pr s proces biolrg ess og ism proc nt or sponical regu es ganire se za lations deve tion or to stim lopm biog ulus multic ental enes ellula is pr r orga oc loca es pr nismal prlizations mulre uctiv ti-orod es e proc ganism ocess immun proces s e syst s em prgrowth rhythm presyn ic ocess aptic locoprmocess proces otion biolog si s invl ical adgnalin oved g reprodhesion in ch emic uctio al sy ha n naptic cebe ll killlvior transm ing de is si to on cell ag xificat greg ion biolog ation ical ph ase cell pa rt orgaorganelle ne m lle pa marco embrembranrt molm an ecular e pa e extrac compl rt el lular supram ex gion extrac oleculcell jure nctio mem n compl bran ellular ar e-en region closed paex rt lumen mitoch ce ondr synanucleoidll ion-as pse sociat othe rt ed ad r orga synapa pse sm pa viralherensnico rt occlus mpl ion ex dy other symbo orgaplast nuclei virionnism c acid pa vi rt bindin bindrion g tran scrip catalying tion fa tic molec struct tran ctor ular fu ural m sporte ction olecul r sign regu e transc al tran lato riptio antiosducerr n fact or acmoleculeletron xidant tivity ar tran carrie , prot ein sducerr nu nding mettralient rebi loch servoi r translchemoraperon ation epelle e chem regu nt D-alaoattraclatator nyl ca nt morph rrier proteiogen n tag cellu

Percent of genes

E 0

100

Biological_process

100

Biological_process Cellular_component

A5_A1 Up

75

50 1490.5 1261.5

25

0

Cellular_component

A8_A1 Up

75

50 1711.5 1426

25

0

Cellular_component 0 0 cellula r proc single metab es -o an olic pr compo oc s biolrg og ism process nent es organiresponical regu zatio se to st lations deve n imulus lopmor biogen en es ta is multic l ellulareproduc locaprlizocess r orga tive ation multiproc ni ess al proc orgasm nism immun s proces ess e syst gr em pr owth rhythm ic ocess locoprmocess presyn biolog io aptic signot ical ad proces alinn g s invl reprodhesion oved uc detoxi in ch fic tion emic al syna be ation vior ptic tr cell ha ansmkillling cell ag is biolog gregatsion ical ph ion ase cell pa ga rt orgaor ne nelle memlle part marco m molecembranbrane ular co e part extrac ellula mplex mem br ece r regi ex an cl ll junc on acelen suprtram lu osed lumtion en oleclaulr region ar co part mitoch mplex ondrio synanucleoid n-asso pse pa ciated rt ad syna cell viralherens co pse occlus mpl ion boex dy other symplas orga t virionnism other pa organi viriort nuclei sm pa n c acid bind rt bindin catalying g tran tic scriptransporte molec struct n fact r ular fu uraltio ction moleculor signal re e transc trangulator riptio antiosducer n fact el xidant or acmolecularetron ca tivity, transd rrier protei ucer nutrie n bind m al nt rese ing lochap rvoir tranet slatio er chemn regulaone orepel tor chem lent D-alaoattracta nyl ca nt rrier morph proteiogen n tag cellula r

25

2981 Down 2523

0 0

Molecular_function

Molecular_function

Percent of genes 100

Molecular_function

D Biological_process

100

Biological_process

A3_A1 Up

75

50 1288 1118

25

0

Cellular_component

A7_A1

75

Up

50

25

1545 1291

0

Cellular_component

Number of genes

B

Number of genes

Biological_process 1117 Down 702

Percent of genes

558.5 351

cellu lar single metab oc -o an olic pr s proces biolrg ess og ism proc nt or sponical regu es ganire zatio se to lations deve n or stimul lopm biogen us multic ental ellula esis pr r orga oc nismlocalizatess re al pr ion multi-produc orga tive process nism oces immun process e syst grow s em rh ythm procesth presyn ic s aptic locoprmocess proces biolog ion signot s invl ical ad aling oved reprodhesion in ch emic uc tio al sy ha n naptic cebe ll killlvior transm ing toxificission cellde ation aggr biolog egatio ical ph n ase cell pa or rt orga ganelle ne m lle pa marco embrembranrt molm ecular ane pa e extrac rt mem ellulacomplex br ce r regi eex an cl ll junc on acelen suprtram lu osed lumtion oleclaulr region en pa ar mitoch compl rt ex ondr nucleo ion-as id sociat syna pse cell ed ad synapart viralherens co pse occlus mpl ion boex dy other sympl orga ast virionnism other part nuclei orga c acid nismvirion bindin bindpart g tran scrip catalying tion fa tic molec struct tran ctor ular fu ural m sporte ction olecul r sign regu e transc al tran la riptio sd tor molec n fact tio ucer ular an or ac tr xidant tivity eletroansducer , prot n ca ein bi rrier trient trannu nd sl io rese ing rvoir metat allonchregulato ap chem r oreperone chem ellent D-alaoattracta nyl ca nt morph rrier proteiogen n tag

50

Number of genes

75

mpone

Percent of genes

Up

Number of genes

cellu lar pr single metab oces -o an olic pr resprg ocess onseism proc s biolog lar co to im ess deve mpone ul lopmical rest ental gulatious orga mulnt proc n ticellu nizatio lar or n or localizatess re ganismbiogen ion multi-productival procesis or imm e process nism une ga system process pr ess rhythm grocess owth ic proc reprod es uctio s be ha n presyn gnalviinor detosi aptic xificat g proces ion ce s invl killlin locollm oved cell ag g in ch io biol gregot emic ation al sybiologicogical ph n naptic al ad ase he transm sion ission ce orgall part lle memmembrne bran ane e pa extracorgane ellula lle part marco r regi rt molec cell ju on extrac nctio mem ul ar el lu compl n br e- lar region supran ex amolenclosed part mitoch ecular lum ondr complen ion-as ex nucleo sociat othe id ed ad r orga cell sm pa viralherensnico rt occlus mpl synaion boex pse pady sy other mpl rt orga ast virionnism pa nuclei viriort c acid n syna bindin bindpse g tran scrip catalying tion fa tic molec sign tran ctor ular fu al transporte ction sduc r transc regu er riptio molec n fact tio lator ular an or ac tr xidant tivity eletroansducer protei n carr stru, ct n bind ier ur nutrieal molecing nt ul translchemorreservoie ation epelle r nt chem re metal oattrgulator ac nt D-alolachaperta nyl ca one morph rrier proteiogen n tag cellu

A1_A0

lar co

Percent of genes

C 100

cellu

mpone

A

Number of genes

lar co

cellu lar single metab oc -o an olic pr s proces biolrg ess og ism proc nt or sponical regu es ganire se za lations deve tion or to stim lopm biog ulus ental enes multic is pr re oc ellula prod loca s r orga uctive lizates multipr ion ni orgasmal process nism oces immun process e syst em prgrowths rhythm presyn ic ocess aptic locoprmocess proces biolog ion signot s invl ical ad aling oved reprodhesion in ch emic uctio al sy n havior naptic cebe ll ki transm llling detoxi is cell ag fic sion gregation biolog ation ical ph as cell pa e rt orgaorganelle nelle marco rt memmembrpa molec bran ane e pa extracular co ellula mplexrt mem br ce r regi eex an cl ll junc on acelen suprtram lu osed lumtion en oleclaulr region ar co part mitoch mplex ondr nucleo ion-as id sociat syna cell pse pa ed ad he sy viral rens conapsrt occlus mpl e ion boex dy sy other mpl orga ast virionnism other part nuclei orga c acid nismvirion bindin bindpart g tran scrip catalying tion fa tic molec struct tran ctor ular fu ural m sporte ction olecul r sign regu e transc al tran lato riptio antiosducerr n fact or acmoleculeletron xidant tivity ar tran carrie , prot ein sducerr nu nding mettralient rebi loch servoir translchemoraperon ation epelle e regu nt chem D-alaoattraclator ny ca tant rrier morl ph proteiogen n tag cellu

ZHU Ying-chun et al. Journal of Integrative Agriculture 2020, 19(2): 528–540 535

2576 Down 2236

0 0

Molecular_function

3090 Down 2582

0 0

Molecular_function

3423 Down 2852

0 0

Fig. 5 Gene ontology (GO) enrichment for differentially expressed genes (DEGs) upon heat stress in five comparison groups: A1_A0 (A), A3_A1 (B), A5_A1 (C), A7_A1 (D) and A8_A1 (E), respectively. The DEGs caused by heat stress matched various GO categories, including biological process, cellular components, and molecular function. A0, unfertilized ovaries heat shocked for 0 h before flowering; A1, A3, A5, A7, and A8, unfertilized ovaries heat shocked for 0, 4, 8, 12, and 24 h, respectively, at 37°C at the first day of the flowering stage.

function. In biological process genes in the A8_A1 group,

involved in the cellular process, cell part and binding were

DEGs were mainly mapped to cellular process (3 963),

metabolic process (3 439), single-organism process

the most significantly regulated by heat stress.

KEGG analysis revealed that 4 pathways in the A1_A0

(2 541), biological regulation (1 806) and response to stimulus

group, 5 pathways in the A3_A1 group, 2 pathways in

(1 501) functions. Of the cellular component DEGs, most

the A5_A1 group, 4 pathways in the A7_A1 group and 9

of those in the A8_A1 pair were mapped in cell part (4 889)

pathways in the A8_A1 group were significantly changed

and organelle (3 082) function. Most DEGs of the molecular

(Fig. 6). Based on significant change of DEGs, it should

function in the A8_A1 pair were mapped in binding (3 537)

be noted that cysteine and methionine metabolism (18,

and catalytic (3 124) functions (Fig. 5). Taken together, there

map00270), carbon fixation (15, map00710), ribosome

existed similarity among the comparison groups that DEGs

(41, map03010) and arginine and proline metabolism (16,

536

ZHU Ying-chun et al. Journal of Integrative Agriculture 2020, 19(2): 528–540

A

B

The enrichment of KEGG

●

Sulfur metabolism

Ribosome

Q_value 0.05 0.04 0.03 0.02 0.01 0.00

Carbon fixation in photosynthetic organisms

Size 1 2 3 4 5 6

Starch and sucrose metabolism

Description

Size 3 4 5 6

Cysteine and methionine metabolism

Description

The enrichment of KEGG

Plant hormone signal transduction

Q_value 0.05 0.04 0.03 0.02 0.01 0.00

Glycine, serine and threonine metabolism

Arginine and proline metabolism

Biosynthesis of amino acids 2.0

2.4

2.8

1.50

Rich_Ratio

C

1.75

2.00

Rich_Ratio

D

The enrichment of KEGG

The enrichment of KEGG

Glycine, serine and threonine metabolism

Description

Size 3 4 5 6 Q_value 0.05 0.04 0.03 0.02 0.01 0.00

Description

Ribosome

Glycine, serine and threonine metabolism

●

Size 2 3 4 5 6 Q_value 0.05 0.04 0.03 0.02 0.01 0.00

Biosynthesis of amino acids

Biosynthesis of amino acids

Alanine, aspartate and glutamate metabolism

1.44

1.48

1.52

1.56

1.4

Rich_Ratio

1.5

1.6

Rich_Ratio

The enrichment of KEGG

E Starch and sucrose metabolism Ribosome

Description

Plant hormone signal transduction

Size 2 3 4 5 6

Phenylpropanoid biosynthesis

Q_value

Glycolysis/ Gluconeogenesis

0.05 0.04 0.03 0.02 0.01 0.00

Glycine, serine and threonine metabolism Carbon fixation in photosynthetic organisms Biosynthesis of amino acids Amino sugar and nucleotide sugar metabolism 1.3

1.4

Rich_Ratio

1.5

Fig. 6 Specific significantly enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways for differentially expressed genes (DEGs) in five comparison groups: A1_A0 (A), A3_A1 (B), A5_A1 (C), A7_A1 (D) and A8_A1 (E), respectively. By correcting the P-value, the metabolic pathway with the smallest Q-value is selected to increase the reliability of the metabolic pathway enrichment. A0, unfertilized ovaries heat shocked for 0 h before flowering; A1, A3, A5, A7, and A8, unfertilized ovaries heat shocked for 0, 4, 8, 12, and 24 h, respectively, at 37°C at the first day of the flowering stage.

537

ZHU Ying-chun et al. Journal of Integrative Agriculture 2020, 19(2): 528–540

map00330) signaling pathways played key roles in response to the flowering phase in A1_A0 (Fig. 6-A). Under different periods of heat stress, in the comparison groups A3_A1, A5_A1, A7_A1, and A8_A1, the pathways of biosynthesis of amino acids signaling (map01230) and glycine, serine and threonine metabolism signaling (map00260) were found to be significantly enriched by P-value analysis (Fig. 6-B–E). Indeed, more regulated DEGs enriched in KEGG pathways occurred in longer periods of heat stress (group A8_A1). These were assigned to the top 3 functions as follows: plant hormone signal transduction (118, map04075), ribosome (112, map03010) and biosynthesis of amino acids (99, map01230). qRT-PCR

3 2 1

6

A3_A1 A5_A1 A7_A1 A8_A1 Cla007544

5 4 3 2 1 0

3.0

1.5 1.0 0.5

–2 –3 –4 –5

Log2FoldChange

4.5 4.0 3.5 3.0 2.5 2.0 1.5 1.0 0.5 0

3 2 1

10

4 2

8 6 4 2 0

A3_A1 A5_A1 A7_A1 A8_A1

A3_A1 A5_A1 A7_A1 A8_A1 Cla019599

12

6

Cla013275

4

12

8

A3_A1 A5_A1 A7_A1 A8_A1

5

14 10

Cla004923

6

0

A3_A1 A5_A1 A7_A1 A8_A1

0

–6

5.0 4.5 4.0 3.5 3.0 2.5 2.0 1.5 1.0 0.5 0

7

Cla018349

Log2FoldChange

Log2FoldChange

0

Cla011760

2.0

Cla017291 A3_A1 A5_A1 A7_A1 A8_A1

–1

A3_A1 A5_A1 A7_A1 A8_A1

2.5

0

A3_A1 A5_A1 A7_A1 A8_A1

Cla004223 Log2FoldChange

4

RNA-seq

Log2FoldChange

5

9 8 7 6 5 4 3 2 1 0

To test the validity of the data generated using the Illumina sequencing platform, 10 up-regulated and 1 down-regulated DEGs were analyzed by qRT-PCR (Fig. 7). Among the 10 selected up-regulated genes, Cla018349 (XP_004143369.1) exhibited the greatest enhancement when analyzed by both methods. Only Cla000855 (XP_004151060.1) and Cla013275 (XP_008463857.1) showed an obvious difference between the qRT-PCR and RNA-Seq methods. The generally similar results from the PCR analysis to the RNA-seq analysis confirmed the RNA-seq analysis results.

Log2FoldChange

Log2FoldChange

6

0

Log2FoldChange

Cla000855

Log2FoldChange

Log2FoldChange

7

3.4. qRT-PCR analysis

A3_A1 A5_A1 A7_A1 A8_A1

Cla021189

A3_A1 A5_A1 A7_A1 A8_A1

Fig. 7 Validation of differentially expressed genes (DEGs) by qRT-PCR in four comparison groups: A3_A1, A5_A1, A7_A1, and A8_A1. A1, A3, A5, A7, and A8, unfertilized ovaries heat shocked for 0, 4, 8, 12, and 24 h, respectively, at 37°C at the first day of the flowering stage.

538

ZHU Ying-chun et al. Journal of Integrative Agriculture 2020, 19(2): 528–540

4. Discussion Although heat stress has been identified as the major trigger in increasing the efficiency of gynogenesis for Cucurbitaceae (Kwack and Fujieda 1988; Gémesjuhász et al. 2002; Shalaby 2007), understanding of the related mechanisms in the physiological process is still limited. In our study, we utilized RNA-Seq to explore the transcriptomic changes of unfertilized ovaries with different periods of heat stress. The effect of different metabolic pathways on the production of ovules can be found during the optimal developmental stage of the female gametophyte. Some studies have shown that embryo sacs in or near maturity induced parthenogenesis more easily than in immature embryo sacs (Zhang et al. 2009; Zou et al. 2017). Further, the inductivity of unfertilized ovaries at 1 d pre-flowering was the highest in in vitro gynogenesis (Niu 2012). Similarly, we considered that the unfertilized ovary sample at 1 d post-flowering contributed to increasing ovule production and enlargement rate, which increased the choice of melon inbred lines to improve the breeding process. The occurrence of 2 363 DEGs in the group A1_A0 and 8 219 in the group A8_A1 suggested that a group of genes should participate in the development and maturation of ovules. In addition, DEGs in the group A1_A0 were mainly in the basic substance and energy metabolism and response actions functions related to a series of biological processes, such as the endogenous stimulus process, organic substance process and cellular metabolic process. Significantly changed pathways were present in the group A1_A0 and enriched into cysteine and methionine metabolism, carbon fixation, ribosome and arginine and proline metabolism pathways. These significantly changed pathways of group A1_A0 may be an important metabolic network during antheses that induces the difference in ovule enlargement rate between A1 and A0. The present study characterized an enhancement of serine acetyltransferase (SATase) (Cla003285 and Cla014346) in watermelon in the group A1_A0. As a conjunction step in cysteine biosynthesis and serine metabolism processes, SATase was found to have a catalytic activity on the formation of OAS (O-acetylserine) from acetyl-CoA and serine (Saito 1995). SATase is important for the biosynthesis of cysteine (Saito et al. 1997). Generally, OAS is the precursor of β-pyrazolealanine that is a secondary product and the precursor of cysteine in Curcurbitaceae plants (Brown 1993). SATase is also considered as a rate-limiting factor in the cysteine biosynthesis process (Saito et al. 1994). The gene expression of tyrosine aminotransferase (TAT) (Cla013546 and Cla023356) was significantly up-regulated in this study, which may be a response to conversion of tyrosine to 4-hydroxyphenylpyruvate (Buchanan et al. 2000), which participates in producing various important

ingredients, such as rosmarinic acid, plastoquinone and tocopherols (Riewe et al. 2012). Some previous studies have indicated that spermidine synthase mediated the synthesis of decarboxylated S-adenosyl-methionine (SAM) that is the precursor of aminocyclopropane carboxylic acid (Imai et al. 2004). It is well known that SAM is synthesized with methionine and ATP in plants and is a methyl group donor involved in DNA methylation for regulation of gene expression (Igarashi and Katoh 2012), playing a crucial role in plant development (Evans and Malmberg 1989). Here, the expression of spermidine synthase (Cla004242 and Cla019182) was significantly up-regulated into the cysteine and methionine metabolism pathway in the A1 group. Similarly, gene expression levels of DNA (cytosine5)-methyltransferase 1 (Cla001916), 1-aminocyclopropane1-carboxylate synthase (Cla014652), S-adenosylmethionine decarboxylase (Cla022875) and methionine-gamma-lyase (Cla008888) in the production of sulfur-containing amino acids was also enhanced in this study. Because of amino acids participating in many regulatory processes in flowering plants (Coruzzi and Last 2000), these results indicate that these enzymes which predominantly occurr at flowering (group A1) may function in amino acid synthesis and transformation. Transcriptome analysis of groups A3, A5, A7, and A8, cultured for different periods of time at 37°C, conveyed that with the extension of exposure time, more DEGs were significantly regulated. Among the different groups, DEGs were only enriched into biosynthesis of amino acids and glycine, serine and threonine metabolism signal pathways. Some of them were obviously up-regulated genes, such as hydroxypyruvate reductase (HPR) (Cla021283, Cla021283 and Cla021283) in cooperation with some related enzymes in chloroplasts and mitochondria to regulate the glycolate pathway of photorespiration (Sloan et al. 1993), aspartate kinase (Cla005922) and bifunctional aspartokinase/ homoserine dehydrogenase 1 (Cla011076) synthesizing some essential amino acids (Weisemann and Matthews 1993), and threonine dehydratase/deaminase (TD) (Cla018349, Cla018351 and Cla018352) as the key ratedetermining enzyme to Ile (Mourad and King 1995; Garcia and Mourad 2004). The expressions of phosphoserine aminotransferase (Cla008937), phosphoserine phosphatase (Cla019056), glycine hydroxymethyltransferase (Cla002746) and dihydrolipoamide dehydrogenase (Cla009256) genes increased in glycine and serine and threonine metabolism signal pathways with increasing time under heat stress, suggesting that heat stress played a significant influence on the these enzyme activates to promote some essential amino acid transformation and synthesis. Conversely, the expressions of probable phosphoglycerate mutase (Cla004748 and Cla019504), D-glycerate 3-kinase

ZHU Ying-chun et al. Journal of Integrative Agriculture 2020, 19(2): 528–540

(Cla015413), alanine-glyoxylate transaminase/serineglyoxylate-pyruvate transaminase (Cla006334 and Cla008570), sarcosine oxidase/L-pipecolate oxidase (Cla009971), and threonine aldolase (Cla000778 and Cla013620) were down-regulated in this signal pathway, which implied that the accumulation or transformation of these amino acids may be limited to ensure priority of the synthesis of other substances. In principle, plant morphogenesis and stress adaptation needs a variety of amino acids acting as a precursor of secondary metabolites and hormones (Jander et al. 2004). Therefore, the change of the above-mentioned genes was targeted for increasing amino acid contents for embryoric fission and embryo production. Among DEGs, primary-amine oxidase (Cla009405), amidase (Cla020306), 4-coumarate-CoA ligase (Cla002689) and caffeoyl-CoA O-methyltransferase (Cla016012) enriched into phenylalanine metabolism (map00360) were up-regulated existed in all comparable groups. Several studies have reported that many plants respond to pathogen infection or environmental stress conditions by the induction of their phenylpropanoid metabolism (Grimmig and Matern 1997), and plants also utilize it for inducing the reinforcement and lignification of the cell wall (Nicholson and Hammerschmidt 1992; Matern and Grimmig 1994). As a major component of secondary cell walls, lignin provides mechanical support and could also provide a defense for plants against pathogen attack through increased lignification (Senthil-Kumar et al. 2010; Bi et al. 2011). Caffeoyl-CoA O-methyltransferase is a key metabolite in the biosynthesis of guaiacyl lignin, specifically for the methylation of caffeoyl CoA (Wang et al. 2017). In our study, we propose that the gene changes were mainly a response to stress to protect plant cells from heat damage.

5. Conclusion We analyzed transcriptomic data from sequencing and found amounts of DEGs in response to treatment at 37°C in the unfertilized ovaries of watermelon. Gene Ontology and KEGG analyses showed similarities and differences in the regulation profiles of the molecular mechanism of heat stress. In addition, the expression patterns for a total of 11 genes were verified by qRT-PCR. The gene transcription data suggested some unrecognized changes, broadening our understanding of ovule production in response to temperature and exposure time, and provided further experimental direction for understanding the mechanism of in vitro gynogenesis.

Acknowledgements This research was supported by the earmarked fund

539

for China Agriculture Research System (CARS-25), the Fundamental Research Funds of the Chinese Academy of Agricultural Sciences (Y2018YJ15 and Y2019XK16-03), the Agricultural Science and Technology Innovation Program, Chinese Academy of Agricultural Sciences (CAAS-ASTIP2018-ZFRI), and the National Key R&D Program of China (2018YFD0201310). Appendices associated with this paper can be available on http://www.ChinaAgriSci.com/V2/En/appendix.htm

References Bi C, Chen F, Jackson L, Gill B S, Li W. 2011. Expression of lignin biosynthetic genes in wheat during development and upon infection by fungal pathogens. Plant Molecular Biology, 29, 149–161. Brown E G. 1993. Enzymes of N-heterocyclic ring synthesis. In: Lea P J, ed., Methods in Plant Biochemistry: Vol. 9. Enzymes of Secondary Metabolism. Academic Press, London. pp. 153–181. Buchanan B B, Gruissem W, Jones R L. 2000. Biochemistry & Molecular Biology of Plants. John Wiley & Sons, USA. Chambonner D, de Vaulx R D. 1985. Obtention of embryos and plants from in vitro culture of unfertilizered ovules of Cucurbita pepo. Cucurbit Genetics Coop, 8, 66. Coruzzi G M, Last R L. 2000. Amino acids. In: Buchanan R B, Gruissem W, Jones R, eds., Biochemistry and Molecular Biology of Plants. American Society of Plant Physiology Press, Rockville, MD. pp. 358–410. Du Q, Wang K, Xu C, Zou C, Xie C, Xu Y, Li W X. 2016. Strandspecific RNA-Seq transcriptome analysis of genotypes with and without low-phosphorus tolerance provides novel insights into phosphorus-use efficiency in maize. BMC Plant Biology, 16, 222. Evans P T, Malmberg R L. 1989. Do polyamines have roles in plant development? Annual Review of Plant Physiology and Plant Molecular Biology, 40, 235–269. Garcia E L, Mourad G S. 2004. A site-directed mutagenesis interrogation of the carboxy-terminal end of Arabidopsis thaliana threonine dehydratase/deaminase reveals a synergistic interaction between two effector-binding sites and contributes to the development of a novel selectable mark. Plant Molecular Biology, 55, 121–134. Gémesjuhász A, Balogh P, Ferenczy A, Kristóf Z. 2002. Effect of optimal stage of female gametophyte and heat treatment on in vitro gynogenesis induction in cucumber (Cucumis sativus L.). Plant Cell Reports, 21, 105–111. Grimmig B, Matern U. 1997. Structure of the parsley caffeoylCoA O-methyltransferase gene, harbouring a novel elicitor responsive cis-acting element. Plant Molecular Biology, 33, 323–341. Igarashi K, Katoh Y. 2012. Metabolic aspects of epigenome: Coupling of S-adenosylmethionine synthesis and gene regulation on chromatin by SAMIT module. Sub-Cellular Biochemistry, 61, 105–118. Imai A, Matsuyama T, Hanzawa Y, Akiyam T, Tamaoki M,

540

ZHU Ying-chun et al. Journal of Integrative Agriculture 2020, 19(2): 528–540

Saji H, Shirano Y, Kato T, Hayashi H, Shibata D. 2004. Spermidine synthase genes are essential for survival of Arabidopsis. Plant Physiology, 135, 1565–1573. Jander G, Norris S R, Joshi V, Fraga M, Rugg A, Yu S, Li L, Last R L. 2004. Application of a high-throughput HPLC-MS/ MS assay to Arabidopsis mutant screening; evidence that threonine aldolase plays a role in seed nutritional quality. The Plant Journal, 39, 465–475. Kwack S N, Fujieda K. 1988. Somatic embryogenesis in cultured unfertilized ovules of Cucurbita moschata. Engei Gakkai Zasshi, 57, 34–42. Livak K J, Schmittgen T D. 2001. Analysis of relative gene expression data using real-time quantitative PCR and the 2–ΔΔCt method. Methods, 25, 402–408. Lotfi M, Alan A R, Henning M J, Jahn M M, Earle E D. 2003. Production of haploid and doubled haploid plants of melon (Cucumis melo L.) for use in breeding for multiple virus resistance. Plant Cell Reports, 21, 1121–1128. Malik A A, Li C, Zhang S X, Chen J F. 2011. Efficiency of SSR markers for determining the origin of melon plantlets derived through unfertilized ovary culture. Horticultural Science, 38, 27–34. Matern U, Grimmig B. 1994. Natural phenols as stress metabolites. Acta Horticulturae, 381, 448–462. Metwally E I, Moustafa S A, El-Sawy B I, Haroun S A, Shalaby T A. 1998a. Production of haploid plants from in vitro culture of unpollinated ovules of Cucurbita pepo. Plant Cell Tissue and Organ Culture, 52, 117–121. Metwally E I, Moustafa S A, El-Sawy B I, Shalaby T A. 1998b. Haploid plantlets derived by anther culture of Cucurbita pepo. Plant Cell Tissue and Organ Culture, 52, 171–176. Mourad G, King J. 1995. L-O-methylthreonine-resistant mutant of Arabidopsis defective in isoleucine feedback regulation. Plant Physiology, 107, 43–52. Nicholson R L, Hammerschmidt R. 1992. Phenolic compounds and their role in disease resistance. Annual Review of Phytopathology, 30, 369–389. Niu M M. 2012. Study of gynogenesis in vitro to induce haploid plants in melon. MSc thesis, Northeast Agricultural University, China. Riewe D, Koohi M, Lisec J, Pfeiffer M, Lippmann R, Schmeichel J, Willmitzer L, Altmann T. 2012. A tyrosine aminotransferase involved in tocopherol synthesis in Arabidopsis. The Plant Journal, 71, 850–859. Saito K. 1995. Cysteine biosynthesis as a sulfur assimilation pathway in plants: Molecular and biochemical approach. In: Mathis P, ed., Photosynthesis: From Light to Biosphere. Kluwer Academic Publishers, Dordrecht, The Netherlands. pp. 347–352. Saito K, Inoue K, Fukushima R, Noji M. 1997. Genomic structure and expression analyses of serine acetyltransferase gene in Citrullus vulgaris (watermelon). Gene, 189, 57–63. Saito K, Kurosawa M, Tatsuguchi K, Takagi Y, Murakoshi I. 1994. Modulation of cysteine biosynthesis in chloroplasts of transgenic tobacco overexpressing cysteine synthase [O-acetylserine(thiol)-lyase]. Plant Physiology, 106,

887–895. San N L H. 1979. In vitro induction of gynogensis in higher plants. Proceedings of the Conference on Broadening Genetic Base of Crops, 51, 327–329. Senthil-Kumar M, Hema R, Suryachandra T R, Ramegowda H V, Gopalakrishna R, Rama N, Udayakumar M, Mysore K S. 2010. Functional characterization of three water deficit stress-induced genes in tobacco and Arabidopsis: An approach based on gene down regulation. Plant Physiology and Biochemistry, 48, 35–44. Shalaby T A. 2007. Factors affecting haploid induction through in vitro gynogenesis in summer squash (Cucurbita pepo L.). Scientia Horticulturae, 115, 1–6. Sloan J S, Schwartz B W, Becker W M. 1993. Promoter analysis of a light-regulated gene encoding hydroxypyruvate reductase, an enzyme of the photorespiratory glycolate pathway. The Plant Journal, 3, 867–874. Song J, Feng S J, Chen J, Zhao W T, Yang Z M. 2017. A cadmium stress-responsive gene AtFC1 confers plant tolerance to cadmium toxicity. BMC Plant Biology, 17, 187. Trapnell C, Williams B A, Pertea G, Mortazavi A, Kwan G, Van M B, Salzberg S L, Wold B J, Pachter L. 2010. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nature Biotechnology, 28, 511–515. Vinoth A, Ravindhran R. 2016. Efficient plant regeneration of watermelon (Citrullus lanatus Thunb.) via somatic embryogenesis and assessment of genetic fidelity using ISSR markers. In Vitro Cellular & Developmental BiologyPlant, 52, 107–115. Wang L, Feng Z, Wang X, Wang X, Zhang X. 2010. DEGseq: An R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics, 26, 136–138. Wang Z, Ge Q, Chen C, Jin X, Cao X, Wang Z. 2017. Function analysis of caffeoyl-CoA O-methyltransferase for biosynthesis of lignin and phenolic acid in Salvia miltiorrhiza. Applied Biochemistry and Biotechnology, 181, 562–572. Weisemann J M, Matthews B F. 1993. Identification and expression of a cDNA from Daucus carota encoding a bifunctional aspartokinase-homoserine dehydrogenase. Plant Molecular Biology, 22, 301–312. Yang N, Yang M Y, Wang W Y, Wu J X. 2019. Study on supply response and demand response of watermelon in China based on Nerlove model. China Cucurbits and Vegetables, 32, 50–53. (in Chinese) Zhang L L, Cui C S, Qu S P. 2009. Progress of plants in vitro gynogenesis. Journal of Northeast Agricultural University, 40, 133–136. (in Chinese) Zhu Q, Gao P, Liu S, Zhu Z, Amanullah S, Davis A R, Luan F. 2017. Comparative transcriptome analysis of two contrasting watermelon genotypes during fruit development and ripening. BMC Genomics, 18, 3. Zou T, Sun B, Wang Z, Gong S, Chen Z, Sun X. 2017. Review of in vitro gynogenesis of Cucurbitaceae plants. China Cucurbits & Vegetables, 30, 1–5. (in Chinese) Executive Editor-in-Chief Huang San-wen Managing editor WENG Ling-yun