Genome-Wide Analysis of the Peroxidase Gene Family and Verification of Lignin Synthesis-Related Genes in Watermelon

Watermelon (Citrullus lanatus) is an important horticultural crop worldwide, but peel cracking caused by peel hardness severely decreases its quality. Lignification is one of the important functions of class III peroxidase (PRX), and its accumulation in the plant cell wall leads to cell thickening and wood hardening. For in-depth physiological and genetical understanding, we studied the relationship between peel hardness and lignin accumulation and the role of PRXs affecting peel lignin biosynthesis using genome-wide bioinformatics analysis. The obtained results showed that lignin accumulation gradually increased to form the peel stone cell structure, and tissue lignification led to peel hardness. A total of 79 ClPRXs (class III) were identified using bioinformatics analysis, which were widely distributed on 11 chromosomes. The constructed phylogenetics indicated that ClPRXs were divided into seven groups and eleven subclasses, and gene members of each group had highly conserved intron structures. Repeated pattern analysis showed that deletion and replication events occurred during the process of ClPRX amplification. However, in the whole-protein sequence alignment analysis, high homology was not observed, although all contained four conserved functional sites. Repeated pattern analysis showed that deletion and replication events occurred during ClPRXs’ amplification process. The prediction of the promoter cis-acting element and qRT-PCR analysis in four tissues (leaf, petiole, stem, and peel) showed different expression patterns for tissue specificity, abiotic stress, and hormone response by providing a genetic basis of the ClPRX gene family involved in a variety of physiological processes in plants. To our knowledge, we for the first time report the key roles of two ClPRXs in watermelon peel lignin synthesis. In conclusion, the extensive data collected in this study can be used for additional functional analysis of ClPRXs in watermelon growth and development and hormone and abiotic stress response.


Introduction
The class III peroxidase gene family (PRXs) is a widely distributed isozyme family type that is well known to have various short names, e.g., PRX, POD, POX, and PER, and contributes to multiple significant physiological reactions in plants [1]. Its key functions are involved in causing reduction-oxidation reactions in electrons triggered by H 2 O 2 and other types of organic and inorganic compounds, e.g., eliminating the excess amount of H 2 O 2 produced in plant tissue, facilitating wound healing, cross-linking small molecules of oxidized poly-lignin within the cell wall, and providing protection against destructive insects and pathogens [2][3][4].
oxidase genes were finally screened by further deduplication, motif analysis, and functional domain analysis. The inferred genes obtained by HMMER (PF00141) were characterized by conserved and stable structures, and the searched genes contained two other kinds of results. The Cla97C01G007780 gene was obtained by keyword search and removed due to the absence of multiple motifs and the absence of a conserved peroxidase domain.
The 79 peroxidase genes were named ClPRX01 to ClPRX79 based on their genetic location across the watermelon chromosome. All these genes were distributed alone or in clusters on eleven chromosomes, among which chromosomes 4 and 10 had the least distribution with only two counts. Chromosome 2 showed the most abundant genes, with a maximum number of 22 counts. The length of the coding region was 966-14,306 bp, the length of the CDS region was 753-7195 bp; the number of exons was 1-14; the coding protein size was 250-632 aa; the molecular weight was 27,274.98-71,578.44 Da; the PI value was 4.42-9.57 (there were 36 acidic proteins of PI < 7 and 43 basic proteins of PI > 7). The secondary structure was mainly characterized as an irregular wavy pattern, accounting for 46.84-63.60%. The specific information for the corresponding numbering, chromosome location, gene starting and ending locations, coding region length, CDS length, number of exons, and encoded protein characteristics of each gene is shown below in Table 1.

Chromosomal Localization and Gene Replication Relationship of ClPRXs
A total of 79 ClPRXs were widely distributed over 11 chromosomes (Figure 1). Two chromosomes (4 and 10) exhibited the least distribution with the presence of two gene counts, and chromosomes 7 and 8 showed three gene counts. There was a total of five genes on chromosomes 3 and 5, and chromosome 6 showed the presence of seven genes. There was a maximum of eight genes on chromosome 11, and chromosome 9 had fourteen genes; however, chromosome 2 had the widest distribution, with 22 genes. Among them, genes on chromosomes 1, 2, 6, 9, and 11 were distributed in cluster form, with relatively dense genes.

Chromosomal Localization and Gene Replication Relationship of ClPRXs
A total of 79 ClPRXs were widely distributed over 11 chromosomes (Figure 1). Two chromosomes (4 and 10) exhibited the least distribution with the presence of two gene counts, and chromosomes 7 and 8 showed three gene counts. There was a total of five genes on chromosomes 3 and 5, and chromosome 6 showed the presence of seven genes. There was a maximum of eight genes on chromosome 11, and chromosome 9 had fourteen genes; however, chromosome 2 had the widest distribution, with 22 genes. Among them, genes on chromosomes 1, 2, 6, 9, and 11 were distributed in cluster form, with relatively dense genes.
The internal linearity of the97,103 watermelon genome species showed that ClPRXs displayed 12 pairs of fragment duplicated genes (Table S2) and 25 pairs of tandem duplicated genes (Table S3). After the calculation of the non-synonymous mutation rate (Ka) and synonymous mutation rate (Ks) values, two pairs of fragment duplicated genes (ClPRX03-ClPRX49 and ClPRX04-ClPRX67) were shown whose Ks value was zero All the gene pairs with synonymous mutations also had Ka/Ks ratios of less than one. The results indicated that the intraspecific evolution of ClPRXs existed by purified selection, which was relatively conservative in evolution, eliminating harmful mutations and keeping the translated proteins of ClPRXs unchanged.   The internal linearity of the97,103 watermelon genome species showed that ClPRXs displayed 12 pairs of fragment duplicated genes (Table S2) and 25 pairs of tandem duplicated genes (Table S3). After the calculation of the non-synonymous mutation rate (Ka) and synonymous mutation rate (Ks) values, two pairs of fragment duplicated genes (ClPRX03-ClPRX49 and ClPRX04-ClPRX67) were shown whose Ks value was zero All the gene pairs with synonymous mutations also had Ka/Ks ratios of less than one. The results indicated that the intraspecific evolution of ClPRXs existed by purified selection, which was relatively conservative in evolution, eliminating harmful mutations and keeping the translated proteins of ClPRXs unchanged.

Interspecific Collinearity of PRXs
The collinearity analysis results between watermelon and Arabidopsis revealed 13,489 collinearity gene pairs between the two species ( Figure 2), among which 35 PRX genes of Arabidopsis and 28 ClPRXs constituted 47 collinearity gene pairs (Table S4). The colinear gene pairs of ClPRXs reached a maximum (11 pairs) on Chr01. The number of gene pairs on Chr02, Chr09, Chr10, Chr05, Chr07, Chr03, Chr08, Chr06, Chr011, and Chr04 decreased sequentially, but Chr04 had at least one gene pair. These genes showed high similarity with the same function and might originate from a common ancestor. In addition, we found that Arabidopsis PRX genes can map 1-4 ClPRX homologous genes, indicating that ClPRXs may have fourfold replication in the evolutionary process [22]. After the Ka and Ks value calculations, thirteen pairs of collinear peroxidase genes of watermelon and Arabidopsis with a Ks value of zero were found. Synonymous mutations had occurred and could be synonymous mutations. The remaining genes for Ka/Ks had a value less than one, which showed more conservative results for the ClPRXs in the interspecific evolution method for the purification of choice. The gene pairs with a stable structure and collinearity between species were also more consistent in function.

Interspecific Collinearity of PRXs
The collinearity analysis results between watermelon and Arabidopsis reveale 13,489 collinearity gene pairs between the two species (Figure 2), among which 35 PR genes of Arabidopsis and 28 ClPRXs constituted 47 collinearity gene pairs (Table S4). Th colinear gene pairs of ClPRXs reached a maximum (11 pairs) on Chr01. The number o gene pairs on Chr02, Chr09, Chr10, Chr05, Chr07, Chr03, Chr08, Chr06, Chr011, and Chr0 decreased sequentially, but Chr04 had at least one gene pair. These genes showed hig similarity with the same function and might originate from a common ancestor. In add tion, we found that Arabidopsis PRX genes can map 1-4 ClPRX homologous genes, ind cating that ClPRXs may have fourfold replication in the evolutionary process [22]. Afte the Ka and Ks value calculations, thirteen pairs of collinear peroxidase genes of water melon and Arabidopsis with a Ks value of zero were found. Synonymous mutations ha occurred and could be synonymous mutations. The remaining genes for Ka/Ks had value less than one, which showed more conservative results for the ClPRXs in the inter specific evolution method for the purification of choice. The gene pairs with a stable struc ture and collinearity between species were also more consistent in function.

ClPRXs' Evolutionary Tree, Gene Structure, Protein Sequence Comparison, and Conserved Structure
The constructed phylogenetic tree of ClPRXs showed division into seven group ( Figure 3A). The number of genes in group I was the largest (30 counts), followed by grou III (14 counts), group IV (10 counts), groups II, V, and VII (7 counts), group VI (4 counts and finally, group IV (4 counts). The gene structure, amino acid ratio, conserved mot elements, and functional domain similarity of the genes were higher in the group. Th gene exon-intron structure is an important feature of gene evolution, but the exon struc ture is a significant part of the coding protein. The insertion of introns or structural loss i the process of exon evolution might lead to positive selection, functional loss, or mutatio in gene evolution. There were significant differences in the number of exons in ClPRXs as shown in Figure 3C. The number of exons in 79 ClPRXs ranged from one to fourteen among which sixty-nine genes had three and four exons and were the most commo genes. The structure of the four genes in group VI was relatively stable, with only one t  The constructed phylogenetic tree of ClPRXs showed division into seven groups ( Figure 3A). The number of genes in group I was the largest (30 counts), followed by group III (14 counts), group IV (10 counts), groups II, V, and VII (7 counts), group VI (4 counts), and finally, group IV (4 counts). The gene structure, amino acid ratio, conserved motif elements, and functional domain similarity of the genes were higher in the group. The gene exon-intron structure is an important feature of gene evolution, but the exon structure is a significant part of the coding protein. The insertion of introns or structural loss in the process of exon evolution might lead to positive selection, functional loss, or mutation in gene evolution. There were significant differences in the number of exons in ClPRXs, as shown in Figure 3C. The number of exons in 79 ClPRXs ranged from one to fourteen, among which sixty-nine genes had three and four exons and were the most common genes. The structure of the four genes in group VI was relatively stable, with only one to two exons, and three of them had no intrusions. The seven genes in group VII showed the most complex genetic structure, with a maximum of fourteen exon counts in ClPRX32 gene members and the fewest exons with nine counts in ClPRX21. To display the clear protein structural characteristics of ClPRXs, we combined amino acid comparative, protein structural characteristic, and protein functional domain binding site analyses. In accordance with previous studies, we found that ClPRXs had three highly conserved protein domains (I, II, and II, as shown in Figure 4A), which were the distalheme-binding domain, position domain, and proximal-heme-binding domain, and these domains also corresponded to the conserved protein motifs Motif 1, Motif 2, and Motif 3 in Figure 3D, respectively. From the comparison of the protein sequences of ClPRXs in Figure 4C, the region after the conserved domain of III demonstrated a large difference To display the clear protein structural characteristics of ClPRXs, we combined amino acid comparative, protein structural characteristic, and protein functional domain binding site analyses. In accordance with previous studies, we found that ClPRXs had three highly conserved protein domains (I, II, and II, as shown in Figure 4A), which were the distalheme-binding domain, position domain, and proximal-heme-binding domain, and these domains also corresponded to the conserved protein motifs Motif 1, Motif 2, and Motif 3 in Figure 3D, respectively. From the comparison of the protein sequences of ClPRXs in Figure 4C, the region after the conserved domain of III demonstrated a large difference due to the peroxidase characteristics in the multifunction catalysis. This region might be the catalytic domain specific to the peroxidase (shown as the gray area in Figure 4A). Peroxidase has four main conserved functions, and each functional site is marked in different colors, as shown in Figure 4B. The heme-binding site mainly consists of distal histidine (Hd) and proximal histidine (HP), and the proximal histidine (HP) can connect one Fe element. In addition, the active site, substrate-binding site, and Ca-binding site are responsible for connecting two calcium elements.

Comparative Analysis of PRXs between Watermelon and Arabidopsis
The phylogenetic tree of all genes was constructed and analyzed between the watermelon and Arabidopsis PRX families ( Figure 5), and the tree results were mainly divided into seven large groups and eleven subclasses. According to the results of the evolutionary tree analysis, the PRX gene families of watermelon and Arabidopsis could be classified into 11 subclasses (A, B, C, D, E, F, G, H, I, J, K) to better display the evolutionary and functional relationships between watermelon and Arabidopsis PRX genes. According to the internal evolutionary relationship of ClPRXs, we divided them into seven groups ( Figure 3A). When constructing the evolutionary tree with the Arabidopsis PRX gene family, we found that genes of group II and III ClPRX members had evolutionary relationship changes and showed a scattered distribution. Group II was divided into three subclasses (C, E, and I), and group III was divided into three subclasses (B, G, and J). All seven groups are indicated in Figure 5.
the internal evolutionary relationship of ClPRXs, we divided them into seven groups (Fig-ure 3A). When constructing the evolutionary tree with the Arabidopsis PRX gene family, we found that genes of group II and III ClPRX members had evolutionary relationship changes and showed a scattered distribution. Group II was divided into three subclasses (C, E, and I), and group III was divided into three subclasses (B, G, and J). All seven groups are indicated in Figure 5. The ClPRXs showed a unique group VII (also subclass A) in watermelon, and group VII had the most distant evolutionary relationship. In accordance with the results shown The ClPRXs showed a unique group VII (also subclass A) in watermelon, and group VII had the most distant evolutionary relationship. In accordance with the results shown in Figure 3A, compared with the other six groups, the gene structure of this group contained a peroxidase conserved domain; however, the number of introns in their gene structure was very large, and their conserved protein motifs were mostly missing. There were multiple pairs of homologous genes between PRXs of watermelon and Arabidopsis. The maximum PRXs of subclass K were 30 ClPRX genes and 27 AtPRX genes. Subclass E, subclass I, and subclass B were isolated from group II. Group III contained fewer genes, among which subclass E contained the minimum number of PRX genes, including one ClPRX gene and one AtPRX gene. However, subclass I and subclass B contained three PRX genes, and subclass I contained two ClPRX genes and one AtPRX gene. Subclass B included one ClPRX gene and two AtPRX genes. As shown in Figure 5, the phylogenetic tree was subdivided into 11 subclasses of watermelon and Arabidopsis PRX gene families. The same subclass showed the closest evolutionary relationship and had the most similar functions. The functions of unknown ClPRX members could be further predicted and verified based on the published functions of Arabidopsis PRX gene family members.

Prediction of Subcellular Localization of ClPRX Proteins
As shown in the attached Table S5, ClPRX proteins were mainly located in the chloroplast and cytoplasm. The first six groups of proteins were mainly localized in chloroplasts, which exceeded 82% of the full predicted value of ten. The subcellular localization results for the seventh group of proteins were different and were obviously divided into three types: ClPRX54, ClPRX79, and ClPRX53. These proteins were mainly located in chloroplasts with a predictive integral value of 5.87, followed by cytoplasm with a predictive integral value of 2.67. The subcellular localization of the ClPRX11, ClPRX07, and ClPRX21 proteins was mainly in the cytoplasm with predicted integral values of 4.31, 4.25, and 5.03, respectively, and the ClPRX11, ClPRX07, and ClPRX21 proteins were in vacuoles with predicted integral values of 2.81, 2.97, and 2.85, respectively. The subcellular localization of the ClPRX32 protein was mainly in the cytoplasm and nucleus with predicted integral values of 4.22 and 4.15, respectively.

Functional Gene Ontology and Kyoto Encyclopedia of Genes and Genomes Enrichment of ClPRXs
Gene Ontology (GO) enrichment analysis of all ClPRXs ( Figure 6) showed that the cellular component was mainly the extracellular region, accounting for 71 out of 79. Cl-PRX07, ClPRX10, ClPRX11, ClPRX21, ClPRX32, ClPRX53, ClPRX54, and ClPRX79 were not enriched. All the PRX genes of watermelon were enriched in the GO function of binding in the biological process. The main enriched functions were functional antioxidant activity, tetrapyrrole binding, heme binding, oxidoreductase activity, and oxidoreductase activity, with peroxide acceptor activity accounting for 78 out of 79 total functions, and neither ClPRX32 was enriched. The next most enriched functions were ion binding (total of 73 out of 79), cation binding, and metal ion binding (total of 72 out of 79). Among the molecular functions, all the PRX genes of watermelon were enriched in response to stimulus, cellular process, and biological process. Seven functions were prominently enriched: response to oxidative stress, response to chemicals, cellular detoxification, response to toxic substances, detoxification, cellular response to toxic substances. Cellular oxidant detoxification accounted for 78 out of 79 functions, while ClPRX32 was not detected at all. Cellular metabolic processes and metabolic processes accounted for 74 out of 79 functions. The specific PRX gene GO enrichment information is provided in Table S6. Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis showed that only seventy-five genes out of the seventy-nine ClPRXs were hits, and the main class was involved in six metabolic pathways (Figure 7): the biosynthesis of other secondary metabolites. The phenylpropanoid biosynthesis pathway accounted for 68 out of 75 metabolites: seventy-four metabolites participated in the metabolic pathway, and six out of Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis showed that only seventy-five genes out of the seventy-nine ClPRXs were hits, and the main class was involved in six metabolic pathways (Figure 7): the biosynthesis of other secondary metabolites. The phenylpropanoid biosynthesis pathway accounted for 68 out of 75 metabolites: seventy-four metabolites participated in the metabolic pathway, and six out of seventy-five were mainly involved in metabolism, ascorbate and aldarate metabolism, glutathione metabolism, and the metabolism of other amino acids. The subsequent main classes, A09130 Environmental Information Processing, A09150 Organismal Systems, A09180 Brite Hierarchies, A09130 Environmental Information Processing, and A09180 Brite Hierarchies were removed because of a p-value > 0.05. According to the KEGG enrichment analysis of seventy-nine ClPRXs, four (ClPRX16, ClPRX42, ClPRX56, and ClPRX61) were not involved in any pathway. The KEGG enrichment information for the specific PRX gene families is shown in the attached Table S7. Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis showed that only seventy-five genes out of the seventy-nine ClPRXs were hits, and the main class was involved in six metabolic pathways (Figure 7): the biosynthesis of other secondary metabolites. The phenylpropanoid biosynthesis pathway accounted for 68 out of 75 metabolites: seventy-four metabolites participated in the metabolic pathway, and six out of seventy-five were mainly involved in metabolism, ascorbate and aldarate metabolism, glutathione metabolism, and the metabolism of other amino acids. The subsequent main classes, A09130 Environmental Information Processing, A09150 Organismal Systems, A09180 Brite Hierarchies, A09130 Environmental Information Processing, and A09180 Brite Hierarchies were removed because of a p-value > 0.05. According to the KEGG enrichment analysis of seventy-nine ClPRXs, four (ClPRX16, ClPRX42, ClPRX56, and ClPRX61) were not involved in any pathway. The KEGG enrichment information for the specific PRX gene families is shown in the attached Table S7.

Promoter Elements of ClPRXs
According to the promoter analysis of ClPRX convenience elements (Figure 8), the family genes were mainly involved in environment response (photosensitive, low temperature, drought), hormone (auxin, gibberellin, abscisic acid, salicylic acid) regulation, growth and development (zein metabolic regulation, the biosynthesis of flavonoids), and plant defense and stress reactions, among other activities. It is involved in anti-stress, light
From Table S2, 28 genes were involved in the auxin core, TGA box, and TGA element, which are involved in hormone regulation. Fourteen and thirty-three genes were involved in gibberellin-responsive elements (TATC-box and P-box), respectively. There were 40 genes involved in salicylic acid response elements (TCA elements). There were 48 genes involved in ABRE and 23 genes involved in cis-acting regulatory elements participating in zein metabolism regulation (O2-site). Eight genes were involved in MYB binding sites participating in flavonoid biosynthetic gene regulation (MBSI). The cis-acting element was involved in defense and stress responsiveness (TC-rich repeats) and water responsiveness. Element WUN-motif was made up of twenty-nine and four genes, respectively. Information on the promoter homeopathic elements of specific genes is shown in Table S8.

Expression Analysis of ClPRXs in Watermelon Fruits
Transcriptome data of watermelon fruits at different periods were used to analyze the expression levels of ClPRXs in fruit peel and flesh, and the results showed that 51 genes were expressed in fruit tissues (Figure 9). Gene heatmaps were drawn according to the expression levels, and cluster analysis of PRX gene expression was performed by the expression patterns of ClPRX fruits, which were divided into three large groups (A, B, and C). Group A exhibited highly expressed genes, containing 22 PRX genes, and Group B contained genes with low expression, including 12 PRX genes. Group C contained 45 genes, including genes with no and weak expression. Peroxidase (EC1. 11. 1.7) is considered to be the last step of polymerization in the lignin macromolecule synthesis enzyme reaction because the lignification process can effectively enhance skin hardness, reducing dehiscent fruit, and our results showed that groups A and B contained 34 genes that might be involved in lignin accumulation and some key genes involved in the development of watermelon fruit, which mainly affected the formation of the watermelon peel stone cell structure [14,15].

Observation of the Microstructure, Hardness, and Lignin Content of Watermelon Peel at Different Periods
There were significant differences in peel structure development, peel hardness, and lignin content accumulation in crack resistance at different developmental stages of two watermelon materials. There was also a noteworthy correlation among peel hardness, structure, and endogenous lignin contents ( Figure 10). The peel microstructure analysis of crack-susceptible watermelon "812" displayed no presence of the stone cell structure, while the peel of material "1061" was crack-resistant because there was the presence of a cell structure, which began to develop at 14 DAP ( Figure 10A). We also noticed that peel hardness (g/cm 2 ) of crack-susceptible watermelon "812" did not change relatively (

Observation of the Microstructure, Hardness, and Lignin Content of Watermelon Peel at Different Periods
There were significant differences in peel structure development, peel hardness, and lignin content accumulation in crack resistance at different developmental stages of two watermelon materials. There was also a noteworthy correlation among peel hardness, structure, and endogenous lignin contents ( Figure 10). The peel microstructure analysis of crack-susceptible watermelon "812" displayed no presence of the stone cell structure, while the peel of material "1061" was crack-resistant because there was the presence of a cell structure, which began to develop at 14 DAP ( Figure 10A). We also noticed that peel hardness (g/cm 2 ) of crack-susceptible watermelon "812" did not change relatively ( Figure 10B), but hardness seemed to be stable until 35 DAP. Further, we measured the endogenous lignin contents (mg/g) of both watermelon materials ( Figure 10C) and observed that the lignin content of the fruit peel of crack-resistant material "1061" was significantly higher (0.84 ± 0.07, 0.56 ± 0.06, 1.63 ± 0.15, 1.60 ± 0.15, 1.73 ± 0.18) compared to the fruit peel of crack-susceptible material "812" (0.15 ± 0.01, 0.38 ± 0.04, 0.42 ± 0.04, 0.59 ± 0.05, 1.00 ± 0.11), respectively. However, the lignin content in watermelon material "1061" seemed lower at 14 DAP due to the cell structure differences, but steadily increased from 14-28 DAP; however, the lignin content became stable until 35 DAP. The lignin content in watermelon material "812" gradually increased at each period.

Expression Patterns of ClPRXs in Different Tissues
To further explore the peel lignin accumulation difference, we determined the expression patterns of ClPRX genes between two different watermelon materials. Higher expression was observed in group A and lower expression in group B (Figure 9). A total of 34 ClPRX expression patterns were analyzed in different tissues (leaf, petiole, stem, and peel) by qRT-PCR ( Figure 11). The obtained results indicated that 8 ClPRXs (ClPRX65, Cl-PRX14, ClPRX52, ClPRX79, ClPRX54, ClPRX39, ClPRX53, ClPRX11) were highly expressed in leaves, 12 ClPRXs were highly expressed in petioles (ClPRX67, ClPRX49, ClPRX75, ClPRX36, ClPRX70, ClPRX71, ClPRX34, ClPRX29, ClPRX06, ClPRX01, ClPRX05, ClPRX30), 6 ClPRXs were highly expressed in stems (ClPRX15, ClPRX24, ClPRX27, ClPRX28, ClPRX07, ClPRX51), and 8 ClPRXs were highly expressed in peel tissue (ClPRX66, ClPRX04, ClPRX21, ClPRX22, ClPRX32, ClPRX08, ClPRX35, ClPRX31). These ClPRX genes were strongly expressed in four different tissues, which might play a molecular role in the development and function of different watermelon tissues. Six ClPRXs (ClPRX54, ClPRX28, ClPRX35, ClPRX30, ClPRX06, and ClPRX51) had significantly higher gene expression in "1061" than in "812". Among them, ClPRX06 and ClPRX51 were located in the same phylogenetic tree branch as AT5G66390.1 (AtPRX72) in Arabidopsis (Figure 5), and colinear gene pairs existed (Table S4), further supporting their similar functions. Previous studies have shown that AT5G66390.1 (AtPRX72) plays an im- Six ClPRXs (ClPRX54, ClPRX28, ClPRX35, ClPRX30, ClPRX06, and ClPRX51) had significantly higher gene expression in "1061" than in "812". Among them, ClPRX06 and ClPRX51 were located in the same phylogenetic tree branch as AT5G66390.1 (AtPRX72) in Arabidopsis (Figure 5), and colinear gene pairs existed (Table S4), further supporting their similar functions. Previous studies have shown that AT5G66390.1 (AtPRX72) plays an important role in lignin biosynthesis [16,23]. ClPRX54 is an independent branch of group VII watermelon in the phylogenetic tree, suggesting that ClPRX54 might play a unique role in watermelon crops. These three genes are considered candidate genes for lignin synthesis in the peel of the materials, and further studies of the gene expression during different periods of peel development are needed.

Expression Patterns of Candidate ClPRXs at Different Developmental Stages of Peel
To further study the difference in ClPRX expression levels in different stages of watermelon peel development, we performed qRT-PCR analysis to examine the expression of three candidate genes ( Figure 12). Among the three candidate ClPRXs, both the Cl-PRX51 and ClPRX54 genes were significantly under expressed in the material with less lignin accumulation "812", while they were strongly expressed in material "1061" with more lignin accumulation that could form stone cell tissue, showing a declining trend in different periods. There was a significant difference in lignin accumulation, especially in the early stage (7 DAP, 14 DAP). However, ClPRX06 showed a trend of first increasing and then decreasing in the three periods of the two different materials. The expression level of material "1061" was higher than that of material "812" in the same period, but the difference between varieties in the same period did not reach significance. Therefore, the ClPRX51 and ClPRX54 genes are most likely to be involved in the synthesis of lignin in watermelon peel. in watermelon crops. These three genes are considered candidate genes for lignin synthesis in the peel of the materials, and further studies of the gene expression during different periods of peel development are needed.

Expression Patterns of Candidate ClPRXs at Different Developmental Stages of Peel
To further study the difference in ClPRX expression levels in different stages of watermelon peel development, we performed qRT-PCR analysis to examine the expression of three candidate genes ( Figure 12). Among the three candidate ClPRXs, both the ClPRX51 and ClPRX54 genes were significantly under expressed in the material with less lignin accumulation "812", while they were strongly expressed in material "1061" with more lignin accumulation that could form stone cell tissue, showing a declining trend in different periods. There was a significant difference in lignin accumulation, especially in the early stage (7 DAP, 14 DAP). However, ClPRX06 showed a trend of first increasing and then decreasing in the three periods of the two different materials. The expression level of material "1061" was higher than that of material "812" in the same period, but the difference between varieties in the same period did not reach significance. Therefore, the ClPRX51 and ClPRX54 genes are most likely to be involved in the synthesis of lignin in watermelon peel. The values represent the means ± SDs, and different bars indicate significant differences (p < 0.05) of genes' expression level. Statistical letters (a, b, c, and bc) indicate the significant differences of relative expression of genes in "812" and "1061" watermelon materials.

Discussion
Among the PRXs, class III peroxidases are well known on a large scale in plants; they are usually secreted into the cell wall or into the cytoplasmic base fluid and vacuoles and are widely involved in the physiological process of plant growth and development [23,24]. The PRX members of Arabidopsis, rice, maize, pear, cassava, potato, birch, and other plants have been identified and analyzed, but a detailed study of the watermelon PRX gene family has not been reported to date. Analyses of the genome-wide characteristics of ClPRX based on genome-wide watermelon database information could provide information for molecular improvement of watermelon quality in the future.
In this study, a total of 79 ClPRXs were identified based on the genomic information of Citrullus lanatus Subsp. vulgaris CV. 97103. The number of ClPRXs was mostly similar to that of Arabidopsis (73) [5], but lower than that of rice (138) [6] and maize (119) [7], a difference that can be explained by the status of Watermelon and Arabidopsis as dicotyledons and the close relationship of PRX family members to each other compared with rice and maize. The present results are also consistent with the finding that the PRX gene family extension differs between monocotyledons and eudicotyledons [7] A total of 47 collinear gene pairs were found between 35 PRX genes in Arabidopsis and 28 PRX genes in watermelon, among which the Arabidopsis PRX gene could map 1-4 watermelon PRX genes. Furthermore, a study on the gene replication relationship of ClPRXs reported 12 pairs of fragment replication genes and a total of 25 tandem duplicated genes. According  letters (a, b, c, and bc) indicate the significant differences of relative expression of genes in "812" and "1061" watermelon materials.

Discussion
Among the PRXs, class III peroxidases are well known on a large scale in plants; they are usually secreted into the cell wall or into the cytoplasmic base fluid and vacuoles and are widely involved in the physiological process of plant growth and development [23,24]. The PRX members of Arabidopsis, rice, maize, pear, cassava, potato, birch, and other plants have been identified and analyzed, but a detailed study of the watermelon PRX gene family has not been reported to date. Analyses of the genome-wide characteristics of ClPRX based on genome-wide watermelon database information could provide information for molecular improvement of watermelon quality in the future.
In this study, a total of 79 ClPRXs were identified based on the genomic information of Citrullus lanatus Subsp. vulgaris CV. 97103. The number of ClPRXs was mostly similar to that of Arabidopsis (73) [5], but lower than that of rice (138) [6] and maize (119) [7], a difference that can be explained by the status of Watermelon and Arabidopsis as dicotyledons and the close relationship of PRX family members to each other compared with rice and maize. The present results are also consistent with the finding that the PRX gene family extension differs between monocotyledons and eudicotyledons [7] A total of 47 collinear gene pairs were found between 35 PRX genes in Arabidopsis and 28 PRX genes in watermelon, among which the Arabidopsis PRX gene could map 1-4 watermelon PRX genes. Furthermore, a study on the gene replication relationship of ClPRXs reported 12 pairs of fragment replication genes and a total of 25 tandem duplicated genes. According to the above analysis, some ClPRXs were derived from the purification and selection of the PRX gene in the evolutionary process, and the gene sequence was highly conserved, maintaining the functional stability of the ClPRX protein. The other part exhibited fragment replication and tandem replication during the evolution of the watermelon genome. According to the Arabidopsis PRX gene family, there were at least four replicates, which may be the key reason for the amplification and functional differentiation of ClPRXs [6,25]. Our results provide significant new resources for understanding the evolution of the PRX gene family in different species.
It is well known that multigene family evolution leads to the diversification of gene structures, including protein sequences, exons and introns, promoters, and enhancers [26]. In this study, the structure of the ClPRX gene was studied, and it was found that these genes had a large number and different lengths of gene sequences. The comparison of primary structural sequences revealed a high variability in watermelon materials. Previous studies have also found a homology of total amino acid sequences of PRX gene family members in multiple specimens of less than 35% [5,27], but due to its coding region and motif research, we found that high similarity, indicating a relatively conserved status, was present in a number of structural domains (distal heme combined with domain structure, location structure, and proximal heme combining structural domains), including the heme binding site distal histidine (Hd) and proximal histidine (HP). Eight cysteines (C1-C8) are key amino acids in the disulfide bonds forming the secondary structure of peroxidase [28,29]. These PRX proteins have similar sequence structures only in specific domains, but not in other genetic regions, so we hypothesized that nonhomologous segments define their functions.
In addition, the ClPRX gene family was divided into seven groups based on the phylogenetic tree analysis. The gene structure (exon/intron) and motif composition of the same group were relatively conserved, indicating that genes in the same group might be more similar in function. Then, the ClPRXs and Arabidopsis PRX gene family members were analyzed using the constructed phylogenetic tree so that the ClPRX gene family could be divided into seven groups and eleven subclasses. Based on the existing functions of Arabidopsis on the same branch, we could better predict the function of ClPRXs. The ClPRXs contained a unique group VII (also subclass A) in watermelon, and group VII had the most distant evolutionary relationship. In accordance with the results shown in Figure 3A, compared with the other six groups, the gene structure of this group contained a peroxidase conserved domain, but a large number of introns were present in the gene structure. Moreover, the protein conserved motif was mostly missing. Studies have shown that introns were specifically inserted into the plant genome during plant evolution and were retained [30], which enabled these genes to undergo faster differentiation and pseudogenification [31].
Therefore, PRX genes of subclass A, as genes unique to watermelon, may have more important functions, which requires further study. In addition, members of the watermelon PRX gene promoter region family were found, as well as various gene family members mainly involved in temperature sensing (photosensitive, low temperature, drought), hormone regulation (auxin, gibberellin, abscisic acid, salicylic acid), growth and development (zein metabolic regulation, the biosynthesis of flavonoids), and plant defense and stress reactions, among other activities [32]. Studies have shown that protein sequences with similar structures and functions may exhibit different expression levels due to differences in regulatory sequences, which may be the result of long-term adaptation to environmental changes in plants [12]. Therefore, studying promoter homeopathic elements is one of the keys to comprehensively understanding gene-specific expression.
The accumulation of lignin in the plant cell wall results in thickening of the cell wall and wood (stem) hardening. This lignification is one of the important functions of the third type of peroxidase. In this study, we found that members of the ClPRX gene family were widely involved in the phenylpropanoid biosynthesis pathway in KEGG, accounting for 68/75 genes. Previous studies have shown that peroxidase (EC1.11.1.7) is the last reactive enzyme in the synthesis of large molecules of lignin via the polymerization of small molecules (coumaryl alcohol, coniferyl alcohol, and sinapyl alcohol) of lignin in the phenylpropanoid biosynthesis pathway, and the lignification process of plants occurs throughout the growth period of plants [14,15]. AtPrx2, AtPrx25, AtPrx71, and AtPrx72 play an important role in Arabidopsis lignification [16,17]. In our previous research, we found significant differences in the stone cell structure caused by lignification between the two materials with obvious differences in peel hardness, and no stone cell structure was detected in the materials with low peel hardness, which were susceptible to cracking [22]. In the latest study, herbaceous peony increased the layers of thickened secondary cells and lignin accumulation, resulting in enhanced stem strength and demonstrably straight stems [33]. There was a positive correlation between lignin, cellulose, and stone cell contents in 206 sand pear cultivars [34]. Therefore, we speculated that low lignin synthesis might lead to low peel hardness and easy cracking.
In this study, watermelon peel hardness was significantly correlated with stone cell development and lignin accumulation, and a significant correlation was identified between peel hardness and stone cell structure development at 14 DAP, 21 DAP, and 28 DAP. There was also a significant difference in endogenous lignin content between the two materials in each stage. Basically, stone cell development is the process of lignification caused by lignin accumulation in the cell wall, which can effectively enhance the hardness of tissue. We also obtained two watermelon peel materials with significant differences in endogenous lignin accumulation (28 DAP), which were subjected to qRT-PCR expression analysis. We screened six ClPRX genes (ClPRX54, ClPRX28, ClPRX35, ClPRX30, ClPRX06, ClPRX51). The expression of these genes was significantly higher in the "1061" peel material (with high lignin content) than the "812" peel material (with low lignin content). ClPRX51 and ClPRX54 showed significantly lower expression in the material with less lignin accumulation of "812" material, but strong expression in the material with more lignin accumulation capable of forming stone cell tissue of "1061" material, and they showed a declining trend in different periods. The significant differences were particularly observed in the early stage of lignin accumulation (7 DAP, 14 DAP), which may have represented the key genes affecting lignin synthesis in watermelon peel.
In particular, the ClPRX51 gene and Arabidopsis AT5G66390.1 gene were on the same phylogenetic tree branch (see Figure 5), and colinear gene pairs (Table S2) were identified in both of them, suggesting their similar functions. In previously published studies, knockout of the Arabidopsis gene (AT5G66390.1) exhibited significant reduction of lignin contents compared with the wild-type (WT) [23]. Our study also revealed that the ClPRX51 gene has similar functional effects regulating the lignin synthesis-related pathways in watermelon peel and needs to be further studied.

Plant Materials and Sampling
Two different types of watermelon parent materials, "1061" and "812", with distinct peel characteristics were planted in a plastic greenhouse at XiangYang Experimental Base of Northeast Agricultural University, Harbin, China. The differential peel characteristics of both materials have been described in our previously published study [21].
The plants of two watermelon materials were grown and checked daily for finding the new sprouting buds. When plants reached the full flowering stage, then selected flowers (to be opened the next morning) were marked and covered with a paper cap. The manual pollination was performed in the next morning (6.00 a.m. to 10.00 a.m.), and each pollinated flower was protected with the same paper cap and labeled with the date, to observe the developmental stages and fruit maturity at days after pollination (DAP). Peel materials were sampled at 7 DAP, 14 DAP, 21 DAP, 28 DAP, and 35 DAP, and subsequently, the peel hardness, peel cell microstructure, and peel endogenous lignin contents were observed, to study the lignin accumulation at different developmental stages. For the qRT-PCR analysis of ClPRXs, fruit peels were collected for each key time period of different developmental stages at 7 DAP, 14 DAP, and 28 DAP (when stone cell formation started). In addition, peel materials of 28 DAP and nearby stem, petiole, and leaf materials were obtained for RNA extraction, cDNA synthesis, and qRT-PCR analysis of ClPRXs.
In brief, all the experimental materials from the vegetative and reproductive growth stages of each plant were collected with 3 replications. A total of 1 g leaf, petiole, stem, and peel tissue material was sampled from the marked position of the second female flower. A total 3 undamaged and random samples from each fruit were also collected at 5 different stages with each time interval of 7 DAP until 35 DAP. For fruit peel hardness analysis, three random peel samples (3 cm × 3 cm × 3 cm) were also collected from each fruit, and texture analysis (TA-XT Plus) (Stable Micro System Division, U.K.) was performed. The peel was separated from the fruit flesh and cut into equal samples of 5 mm × 5 mm × 3 mm. Paraffin sectioning was performed by embedding the peel samples into FAA (1:1:18 mixture (v/v/v) of formalin:acetic acid:50% ethanol) fixative solution, and the endogenous lignin contents were visually quantified. In addition, total RNA extraction was also performed from subsequent frozen samples using the TRIzol method [35].

Preliminary Search and Identification of ClPRXs in the Watermelon Genome
Watermelon peroxidase gene family members (ClPRXs) were initially searched for identification, and three methods (searching keywords, HMM, Blast) were used to obtain the gene information and the preliminary screening. The ClPRX genes were termed according to the naming method of Arabidopsis as reported previously [5], e.g., the letter "Cl" for Citrullus lanatus (Thunb.) Matsum. et Nakai and "PRX" for peroxidase and by a number indicating the position of genes on the respective chromosomes.

Chromosome Localization, Gene Replication Relationship, and Interspecific and Intraspecific Collinearity Analysis
The One Step MCScanX module in the TBtools software [26] was used to obtain genome-wide replication events and identify the gene replication relationship of ClPRXs. Genome-wide collinearity analysis was performed among watermelon and watermelon intraspecies, watermelon, and Arabidopsis. The "Amazing Super Circos" module was used to visualize the gene localization and the gene linear relationship on the respective chromosomes. The "Ka/Ks Calculator" module in the TBtools software was used to further calculate the evolutionary relationship between peroxidase genes and fragment replication and tandem replication by incorporating the Ka and Ks of watermelon intraspecific peroxidase genes and comparative interspecific peroxidase genes of watermelon and Arabidopsis. The pressure selection analysis was carried out by calculating the Ka/Ks ratio.

ClPRXs' Structure and Protein Domain Analysis
Complete gene structure information for PRX, including gene length, CDS position, and gene functional domain prediction information, was obtained using the "Gene Structure View" module of the TBtools software. Then, the MEME online database (http://meme-suite.org/ (accessed on 2 April 2021)) was used to analyze the identified watermelon PRX protein conserved motifs. The motif length ranged from 6-200 amino acids; the number of motifs was set to 20; the analysis results were saved. Multisequence alignments of the PRX protein were filtered with MEGA 7.0 software in Muscle [39], and the comparison results were imported into the TBtools module "Quick Run TrimAL". Then, the composition of the parameter selection "NJ_STRICTPLUS" was used to display the results of the amino acid sequence alignment. The NCBI-CD (https://www.ncbi.nlm.nih.gov/Structure/cdd/WRPSB.Cgi (accessed on 2 April 2021)) functional domain and domain site analysis, functional domain binding site information, and NCBI functional module "Cn3D macromolecular structure viewer" were used to demonstrate the functional domain binding sites.

Construction of the Phylogenetic Tree for Watermelon and Arabidopsis
The obtained protein sequences of the watermelon and Arabidopsis PRX families along with homologous protein sequences were arranged in the MEGA 7.0.26 software, and an evolutionary tree was constructed. The neighbor-joining (NJ) adjacent method and 1000 bootstraps were used, and the other parameters were also set as default values for final construction of the evolutionary tree.

Functional Enrichment Analysis of ClPRXs
Gene functional enrichment analysis was performed using the online database (http: //cucurbitgenomics.org/goenrich (accessed on 2 April 2021)) for ClPRX GO annotations using the online TBtools software [26] "GO Enrichment" module downloaded GO annotationbased package hierarchy (GO-basic. Obo) and the GO Enrichment function. KEGG enrichment analysis was performed using the online KEGG database (http://cucurbitgenomics. org/pwyenrich (accessed on 2 April 2021)) for ClPRX annotation and the downloaded KEGG Enrichment using the "KEGG Enrichment" module. Finally, the "Enrichment Bar Plot" module was used to visually display the results.

Transcriptome Analysis of ClPRXs at Different Periods of Peel and Flesh
The online watermelon database (http://cucurbitgenomics.org/ (accessed on 2 April 2021)) and transcriptome registration number SRP012849 [40] were searched, and the downloaded ClPRXs of gene expression data were used to construct the heatmap using the TBtools software.

RNA Extraction, cDNA Synthesis, and Gene Expression Analysis
The samples were lyophilized at 72 h using a lyophilizer, and RNA was extracted using the common TRIzol method [35]. Total RNA was imaged on a 1% agarose gel, and the RNA concentration was detected using a NanoPhotometer ® P330 (IMPLEN, Munich, Germany). Then, the PrimeScript RT Master Mix Perfect Real-Time kit (TOYOBO, Osaka, Japan) was used to synthesize first-strand cDNA with 1 µg RNA. The final cDNA was placed in a -20 • C refrigerator for quantitative RT-PCR (qRT-PCR).
The Cla020175 gene was used as the reference gene, and all the gene primers used in the qRT-PCR analysis are shown in Table S2. The Primer Premier (V6.0) software was used to design the gene-specific primers [41]. Three biological replicates and three technical replicates were used for each cultivar tissue assayed. A 20 µL PCR mixture was prepared with 10 µL of SYBR Green Master mix (TOYOBO, Osaka, Japan), 1 µL of each primer pair, and 1 µL of cDNA templates. PCR amplification of target genes was performed in 96-well optical reaction plates on an iQ5 Gradient Real Time PCR system (Analytik, Jena, Germany). The PCR assay was set as follows: 95 • C for 60 s; 40 cycles of 95 • C for 15 s, 58 • C for 20 s, and 72 • C for 15 s; a final melt curve analysis in which the temperature was increased from 55 • C to 95 • C at a rate of 0.5 • C/5 s; a final hold at 4 • C. The specificity was verified by melt curve analysis and agarose gel electrophoresis of the qRT-PCR products. The relative expression levels of all different genes were determined using the 2 −∆∆CT method [42].

Determination of Watermelon Peel Hardness and Paraffin Sectioning
For the peel hardness, a 2 mm (P/2) diameter probe TA-XT Plus texture instrument was used to conduct puncture tests on the ripened fruit peel of different watermelon varieties as reported previously [21]. The paraffin-embedded tissue section sampling and preparation method were used as previously reported [31]. In brief, the material was embedded into FAA (1:1:18) fixed solution, which was changed for overnight incubation and dehydrated in an alcohol concentration gradient. The xylene was transparent after addition to the paraffin clastic and paraffin embedding at 65 • C using turning-wheel microtome sectioning at a slice thickness of 8 µm, drying under xylene concentration gradient dewaxing, an alcohol concentration gradient, red and green dye fixation, sealing, and imaging.

Endogenous Quantification of Lignin Contents in Watermelon Peel
A total of 1 g of leaf material was ground in liquid nitrogen; 0.1 g was weighed and put into a 2 mL precooled centrifugal tube, and then, the amount of lignin in the watermelon peel was determined by the thioglycolic acid method [43]. Simultaneously, commercial lignin (alkaline spruce lignin, Aldrich, Milwaukee, WI, USA) was used to configure standard samples with different gradient concentrations. The final concentration of all standard samples was determined using the same method, which was used to calibrate the lignin content curve and to calculate the endogenous lignin contents in both watermelon peel materials.

Conclusions
In this study, a total of 79 ClPRX members were identified by genome-wide bioinformatics analysis of watermelon. These identified genes were divided into seven groups and eleven subclasses based on the constructed phylogenetic tree. To better understand these genes, we also analyzed the chromosome distribution, gene replication, gene structure, conserved motifs, GO annotation, KEGG pathway, and promoter element prediction and combined them with the expression analysis of the ClPRX gene family in fruits using online transcriptome data. In addition, we focused on the key role of the peel lignin synthesisrelated ClPRX gene family, and qRT-PCR was used to analyze the expression patterns in different tissues during different development periods. The obtained results showed that deletion and replication events occurred in the process of gene amplification of the ClPRX family, and the overall protein sequences had low homology, but consistently, four conserved functional loci were present. ClPRXs showed different tissue-specific abiotic stress and hormone response expression patterns, providing strong evidence that the ClPRX gene family participates in a variety of different physiological processes occurring in plants.
To our knowledge, this is the first study to predict the key roles of two ClPRXs in peel lignin synthesis. In conclusion, the extensive data collected in this study can be used for additional functional analyses of ClPRXs in watermelon growth and development and hormone and abiotic stress responses.