Identification of POD genes
To identify members of POD family in B. pendula, we used the 73 POD genes of Arabidopsis to obtain the best hits in the B. pendula genome by BLASTP. A total of 90 putative PODs were identified in the B. pendula genome. We further examined the conserve domains of proteins encoded by these genes using Pfam  and SMART  database. The results revealed that all the genes have classical POD domain structures, which demonstrates the reliability of the results. The B. pendula genome contains more PODs than Arabidopsis (73) , but fewer than Populus euphratica (93) , Pyrus bretschneideri (94) , and rice (138) . We defined the BpPODs as BpPOD1 to BpPOD90. The isoelectric point (PI) varied from 4.28 to 9.6 with a mean of 7.25 and >7.0 of 52.2% POD proteins. In addition, subcellular locations of these BpPODs are mainly in the cytoplasm, cell membrane, vacuole, chloroplast and nucleus. Their detailed information, including chromosome location, gene name, subcellular location and molecular weight (MW) gene size of each BpPOD gene/protein, was listed in Table 1.
Phylogenetic analyses of POD proteins in B. pendula
To investigate the evolutionary history and phylogenetic relationships among the members of POD family in B. pendula, a phylogenetic tree was constructed with the Neighbor-Joining method based on multiple sequence alignment of the 90 BpPODs, with 1000 bootstrap replicates (Figure 1). The BpPOD proteins were divided into twelve major subgroups with high bootstrap probabilities, designated group I to group XII. The POD genes of each subgroup is unevenly distributed, with the number of members varies from 4 to 15. Subgroup VIII contains the most members (15), subgroup X, XI, XII contains the least number of members, with only 4 members.
To further gain insights into the structural diversity of the POD genes, we subsequently performed exon-intron analysis in BpPODs. The result reveals several variations (Figure 2), including five BpPODs (BpPOD9, BpPOD11, BpPOD16, BpPOD57 and BpPOD61) lacking an intron in their gene structures. In the remaining BpPODs, the number of introns varies from one to six, while the major members have one to three introns. In addition, BpPOD76 and BpPOD87 have the most introns (6), followed by BpPOD24 and BpPOD51 (5). Furthermore, genes in the same subfamily display similar exon/intron structures. For example, BpPOD20, BpPOD22 and BpPOD82 have three exons and two intron, both of which belong to Group V; BpPOD73 and BpPOD74 have two exons and one intron, both of which belong to Group XI. However, in some cases, the number of exons/introns varies among genes clustered together in the phylogenetic tree. For example, BpPOD52 has one more intron than BpPOD42, BpPOD56 has two more introns than BpPOD43. These differences may be derived from a single intron loss or gain events during the long evolutionary period .
Analysis of conserved amino acid motifs
To understand the functional regions of BpPODs, conserved amino acid motifs analyses of BpPOD proteins were performed. A total of eight conserved amino acid motifs were identified in the BpPOD proteins (Figure 3). All BpPOD proteins contain at least one conserved amino acid motif. For example, BpPOD55 only contains motif 8, BpPOD83 contains motif 1 and 7, while BpPOD10 proteins contain all the eight conserved amino acid motifs.
Most of the closely related members have the same motif compositions, suggesting that there are functional similarities between POD proteins within the same subgroup . We found that motifs 1, 2, 3, 4, 5 and 7 appeared in nearly all members of BpPOD proteins, these motifs might be important for the functions of BpPOD proteins.
Chromosomal location and evolution analysis of BpPODs
To investigate the genome organization and distribution of BpPODs on different chromosomes of B. pendula, a chromosome map was constructed. The results show that the 90 BpPODs were distributed among 14 chromosomes, as shown in Figure 4, the physical locations of these BpPODs on chromosomes were scattered and uneven. Chromosome 1 and 8 contains the most BpPODs (14), followed by chromosome 13 (10). Eight BpPODs were simultaneously distributed on chromosomes 5 and 7, whereas chromosomes 14 had only one and chromosomes 11 does not include the POD gene. In addition, some chromosomes exhibit a relatively high density of BpPODs, such as the bottoms of chromosomes 13 and the top of chromosome 8.
Gene duplication, including segmental and tandem duplication, is considered to be one of the primary driving forces in the evolution of genomes [34, 35]. In this study, among the 90 BpPODs identified, a large number of BpPODs have the same duplicated regions (Figure 5). Generally, a gene cluster is the result of gene tandem duplication . In this study, we found that some BpPODs were adjacent to each other (Figure 4). For instance, BpPOD17-20, BpPOD22-29 and BpPOD11-15 were located sequentially in tandem on chromosomes 5, 8, and 13, respectively, implying that these genes might arise from recent tandem duplication events . The result indicated that tandem duplications play main contributors in the expansion of the BpPOD gene family. The result was consistent with the previously reported for Populus euphratica POD genes, tandem duplications also contributed significantly to the expansion of POD gene family in Populus euphratica . However, in previous studies, many species also have produced some different results. For example, segmental duplication events were the major contributors to the expansion of the pear POD family  and segmental duplication and tandem duplication were identified in maize POD family . These results indicate that there are significant differences in the POD genes expansion pattern in B. pendula, maize and Chinese pear, which strongly implied that POD family members have different expansion patterns among different species. It may be the reason why the POD family members (90) in B. pendula were less than those in the maize (119) .
To explore the selection pressures among BpPOD duplicated genes, we calculated the Ka, Ks and Ka/Ks values for the 23 gene pairs (Table 2). In the general, Ka/Ks > 1 indicates positive selection, Ka/Ks = 1 indicates neutral selection and Ka/Ks < 1 indicates negative or purifying selection . The Ka/Ks ratios of most BpPOD gene pairs were < 1, suggesting that these gene pairs evolved under negative selection in B. pendula. The results of this Ka/Ks analysis suggest that negative selection was vital to the functional divergence of BpPODs.
To further elucidate the evolution mechanisms of BpPODs family, we constructed the comparative syntenic maps of birch associated with representative species containing two dicots (Arabidopsis thaliana and Populus trichocarpa) and one monocots (Vitis vinifera) (Figure 6; Table3). A total of 17 (Arabidopsis thaliana), 49 (Populus trichocarpa) and 43 (Vitis vinifera) orthologous pairs were found between birch and the other species. Interestingly, among these gene pairs, some BpPOD genes (BpPOD3, BpPOD7, BpPOD16, BpPOD21, BpPOD40, BpPOD4, BpPOD48, BpPOD52, BpPOD57 and BpPOD84) were shown to have collinear relationships with all of the above three species, showing that these orthologous pairs might already have been present before the divergence of dicotyledonous and monocotyledonous plants .
Tissue-specific expression of BpPODs
To better understand the functions of POD genes in the growth and development of Betula Platyphylla × Betula Pendula, their expression profiles in different tissues (including root, xylem, young leaf and flower) were analyzed with publicly available transcriptome data. Of the 90 BpPODs, 69 genes were expressed in one or more birch tissues, while 21 BpPOD genes exhibited no expression in various individual tissues. The heat map (Figure 7) demonstrated that most BpPODs had tissue-specific or preferential expression patterns. As shown in Figure 7, BpPOD6, BpPOD21 and BpPOD37 were highly expressed in xylem. Several BpPODs were expressed in root during development, such as BpPOD62, BpPOD63 and BpPOD65. BpPOD78 and BpPOD19 showed higher expression levels in young leaf and flower, respectively. The expression level of BpPOD6 was high in xylem and low in root, leaf and flower. In contrast, BpPOD67, BpPOD68, BpPOD80 and BpPOD81 had no expression in any of the investigated tissues. BpPOD21, BpPOD59 and BpPOD62 were highly expressed in developing xylem, root, leaf and flower. In conclusion, the variations in the expression of BpPODs in different tissues revealed that POD genes may be involved in several processes during B. pendula growth and development.
Responses of BpPODs expression to cold treatment
Several roles have been attributed to plant peroxidases in response to biotic and abiotic stresses . In recent years, the number of studies on POD genes response to abiotic stress have been reported . For example, Arabidopsis overexpressing AtPOD3 showed an increase in dehydration and salt tolerance, whereas the antisense suppression of AtPOD3 exhibited dehydration and salt sensitive phenotypes . The expression of POD genes is induced by various environmental stresses, such as metal, pathogens, humidity, temperature, anoxia and potassium deficiency , suggesting that POD genes are involved in plant defense. In this study, we examined the expression levels of the BpPODs in response to low temperature stress. As shown in Figure 8, the result indicated that the expression of most BpPODs was altered under cold treatment. After cold treatment, the expression levels of BpPOD4, BpPOD13, BpPOD15, BpPOD17 and BpPOD21 were significantly induced at a relatively early stage (0.5 h after treatment), and with the increase of cold treatment time, the relative expression level of these genes was also at a high level. As shown in Figure 8, the expression levels of BpPOD19, BpPOD21, BpPOD39 and BpPOD47 were increased after 1.5 h treatment of low temperature. BpPOD50 and BpPOD58 did not response to cold treatment at the beginning (0.5 h), and were slightly increased after 2 h exposure to low temperature. In general, the low temperature responsive BpPODs may play important roles in birch under cold stress.
Validation of DEGs identified by RNA-seq using qRT-PCR
We used qRT-PCR (Quantitative Real-time PCR) to validate the expression profiles of BpPODs after low temperature treatment. A total of 6 cold-responsive BpPODs were selected for qRT-PCR analysis. As shown in Figure 9, the qRT-PCR results indicated that all the 6 BpPODs were up-regulated by low-temperature stress, which was consistent with the results derived from the RNA-seq data. BpPODs including BpPOD15, BpPOD47 and BpPOD49 showed the highest transcript level when exposed to a low temperature for 1.5h. BpPOD4, BpPOD17 and BpPOD26 showed the highest transcript level at 3h. In general, this qRT-PCR result supports the reliability of the RNA-seq analysis.