Identification of the P/A genes

The whole-genome sequences and annotations of F. vesca Hawaii 4 (v4.0.a2) (Li et al. 2019b) and F. × ananassa Camarosa (v1.0.a1) (Edger et al. 2019) were downloaded from the Genome Database for Rosaceae (GDR) ( https://www.rosaceae.org). Subsequently, the whole-genome nucleotide coding sequences (CDSs, 34 006 CDSs with 2953 bp average length) of F. vesca (Li et al. 2019b) were used as query sequences for BLASTN against the whole-genome sequences of F. × ananassa with a default e-value by using local BLAST+. The query CDSs without BLAST hits were determined to be candidate P/A genes; that is, present in F. vesca but absent in F. × ananassa. In contrast, the CDSs with BLAST hits were considered as candidate non-P/A genes. Furthermore, the distributions of these non-P/A genes in F. nipponica, F. nilgerrensis, F. iinumae, F. viridis and F. nubicola of Fragaria, and Rosa chinensis, Rubus occidentalis, Prunus persica, P. mira, Malus × domestica, Pyrus pyrifolia and P. betulifolia were analyzed by BLASTN searches with a default e-value using non-P/A gene CDSs as the query sequences and every genome sequence as the database (Hirakawa et al. 2014; Edger et al. 2020; Zhang et al. 2020; Feng et al. 2021; Hardigan et al. 2021).

PCR examinations of the P/A genes

To verify the accuracy of identification of the candidate P/A genes, 32 of them were randomly chosen for three-primer PCR amplifications. Two of the three primers for each selected P/A gene were designed to be located in the two flanking sequences of the P/A regions, and the third primer was located in the middle of the P/A regions ( Supplementary Fig. S1 (cjps-2020-0316suppla.pdf) ¹1). Subsequently, the genomic DNA of F. × ananassa Camarosa was extracted and used as the PCR template, and that of F. vesca Hawaii 4 was used as the positive control template. According to the detection on the gel electrophoresis, the true P/A genes could be confirmed based on the production of two PCR products in F. vesca but only one product with a different nucleotide length in F. × ananassa. In contrast, non-P/A genes exhibited the same PCR results in both species.

Classification of multi-genes and single genes

An all-vs.-all BLASTN search was performed on the whole-genome CDSs of F. vesca with a default e-value. The genes were divided into multi-genes families based on the BLAST results with the criteria of coverage ≥60% and identity ≥60%. The remaining genes were defined as single genes. Subsequently, the P/A genes were classified as multi-genes and single genes based on the classification of the genes within the whole genome.

Duplication types of the P/A genes

The different genome-wide duplication types of the genes in F. vesca were defined by using the comparative genomic tool DupGen_finder (Qiao et al. 2019). Afterwards, the duplication types of the P/A genes were detected based on their genome-wide duplication types.

Functional analysis of the P/A genes

The proteins encoded by the identified P/A genes were found by searching the Pfam database ( http://pfam.xfam.org/) with an e-value cutoff equal to 1.0. GO analysis was performed according to the sequencing annotations of F. vesca. KEGG analysis was performed by using the KEGG Automatic Annotation Server (KAAS,  https://www.genome.jp/tools/kaas/) and KEGG Mapper ( https://www.kegg.jp/kegg/mapper.html).

Chromosomal locations of the P/A genes

The chromosomal positions of the P/A genes were determined from the sequencing annotations of F. vesca. Each of the seven chromosomes was divided into different windows in units of 1 Mb, and the gene numbers in each window of all seven chromosomes were counted. The chromosome locations of the P/A genes were performed by MapChart v2.32 software.

Identification and domain preference of the P/A genes

A total of 334 genes were discovered to be present in the F. vesca genome but missing in the F. × ananassa genome. This analysis also demonstrated that 0.98% of the F. vesca genes were absent from the F. × ananassa Camarosa genome. The genomes of the octoploid F. × ananassa were donated by the diploid ancestor species A, B, C, and D (Hardigan et al. 2021). However, except for F. vesca (A) and F. iinumae (B) have been verified as the diploid ancestors of F. × ananassa, the other two subgenomes still remain unknown (Edger et al. 2019; Feng et al. 2021; Hardigan et al. 2021). In this study, the distribution of non-P/A genes was investigated by BLAST+analysis for some diploid Fragaria species used as reference genomes. The non-P/A genes exhibited similar distribution trends in these diploids Fragaria species, including 98.31% (33 104/33 673) of the non-P/A genes in F. nipponica, 96.80% (32 597/33 673) in F. nilgerrensis, 96.52% (32 501/33 673) in F. iinumae, 95.85% (32 275/33 673) in F. viridis and 94.39% (31 784/33 673) in F. nubicola. These results indicated that the non-P/A genes might be commonly prevalent among Fragaria species. In addition, to further excavate whether the non-P/A genes were widely distributed in Rosaceae family, some Rosaceae species also served as the reference genomes to perform BLAST + searches. There were 26 600 (79.00%) non-P/A genes in R. chinensis, 23 382 (69.44%) in R. occidentalis, 16 293 (48.39%) in P. persica, 16 096 (47.80%) in P. mira, 16 020 (47.58%) in M. domestica, 15 769 (46.83%) in P. pyrifolia and 15 659 (46.50%) in P. betulifolia. Therefore, about 46.50%–79.00% of the non-P/A genes widely distributed in the Rosaceae family, and approximately 15.39% to 19.31% of the non-P/A genes were exclusive to Fragaria species.

To verify the true situation of these candidate P/A genes, three-primer PCR amplifications were performed for the 32 randomly selected genes ( Supplementary Table S1 (cjps-2020-0316supplc.xlsx);  Fig. S2 (cjps-2020-0316supplb.pdf)¹). Among the 32 genes, 31 were found to be true P/A genes, and only one was a false P/A gene. For this false P/A gene, we checked the whole-genome sequences of F. × ananassa and found no corresponding nucleotide sequence of this gene. This indicated that the misjudgment of this gene was due to sequencing or assembly errors rather than to our methods. Therefore, the 333 genes were regarded as the candidate P/A genes for further analysis ( Supplementary Table S2 (cjps-2020-0316supplc.xlsx)¹).

The P/A genes were involved in a total of 208 protein domains, in which some specific domains showed a relative high occurrence frequency (Fig. 1). Among the top 20 domains in occurrence frequency, the PPR, PKs, NB-ARC, F-box and EF-hand domains are all involved in the response to biotic or abiotic stresses and adaptation to the environment. For example, an Arabidopsis PPR protein is closely related to ecological adaption (Zsigmond et al. 2008); PKs and NBS-LRR genes are key plant disease-resistance genes (Li et al. 2016; Zhong et al. 2018; Yang et al. 2019); and EF-hand proteins are important in response to ecological stress in the soybean genome (Zeng et al. 2017).

Fig. 1.

The top 20 protein domains of the P/A genes.

Multi-gene families and duplication types of the P/A genes

Based on the criteria of coverage and identity both being larger than 60%, a total of 2667 multi-gene families were defined among the whole genome CDSs of F. vesca. These families contain 7578 multi-genes, representing 22.28% of the whole-genome genes sorted in multi-gene families. Only 27 P/A genes belonged to multi-gene families, indicating that 8.11% (27/333) of the identified P/A genes were multi-genes. The vast majority of them (306/333) were single genes without homologs in F. × ananassa. However, 7551 non-P/A genes were classified as multi-genes, and 26 122 non-P/A genes were single genes. The non-P/A genes classified as multi-genes accounted for 22.42% of all non-P/A genes, which was larger than the proportion of multi-genes among the P/A genes.

The duplication types of F. vesca genes were uncovered at the genome-wide level and included WGDs, tandem duplications and transposed duplications. In total, the smallest number of genes (3038) were WGDs, and 3061 genes were considered transposed duplications. Finally, the maximum number of genes (3385) in the F. vesca genome were found to be tandem duplications (Table 1). Among the identified P/A genes, there were 5, 34 and 15 genes from WGDs, tandem duplications and transposed duplications, respectively. These results demonstrated that 1.50% (5/333), 10.21% (34/333) and 4.50% (15/333) of all the P/A genes were generated by WGDs, tandem duplications and transposed duplications, respectively. In addition, 9.00% (3033/33 673), 9.95% (3351/33 673) and 9.05% (3046/33 673) of all the non-P/A genes were involved in WGDs, tandem duplications and transposed duplications. The majority of P/A genes were tandem duplications compared with the other two duplication types. This demonstrated that tandem duplications might play more important roles in P/A gene expansions than WGDs and transposed duplications.

Table 1.

Number of P/A genes and non-P/A genes in different duplication types.

Functional analysis of the P/A genes

Three categories of GO terms, biological process, cellular component and molecular function, were analyzed. Most P/A genes were categorized in terms of molecular function, followed by the biological process and cellular component categories (Fig. 2). In the category of molecular function, many of the P/A genes exhibited binding functions with proteins, nucleotides, metal ions and other molecules, such as the dominant subcategory of protein binding (23, 6.91%), as well as ATP binding (15, 4.50%), DNA binding (13, 3.90%), nucleic acid binding (7, 2.10%) and zinc ion binding (5, 1.50%). A variety of activities were also in this category, including protein kinase activity (7, 2.10%), oxidoreductase activity (3, 0.90%) and catalytic activity (2, 0.60%), etc. In the category of biological process, the P/A genes seemed to be enriched in protein phosphorylation (7, 2.10%) and oxidation–reduction process (5, 1.50%). Within the cellular component category, more P/A genes were associated with the subcategories of chromosome (8, 2.40%), nuclear pore (2, 0.60%) and ribosome (2, 0.60%).

Fig. 2.

Gene Ontology (GO) analysis of the P/A genes. The abscissa represents the P/A gene numbers in each category.

KEGG pathway analysis showed that the P/A genes were preferentially involved in signal transduction pathways such as the MAPK signaling pathway, phosphatidylinositol signaling system, and plant hormone signal transduction (Table 2). The P/A genes also participated in metabolism-related pathways, including amino sugar and nucleotide sugar metabolism, biosynthesis of secondary metabolites and metabolic pathways. In addition, some P/A genes might take part in the response to pathogens for their clustering in the plant–pathogen interaction pathway.

Table 2.

KEGG pathways of the P/A genes.

Chromosomal locations of the P/A genes

Except for the two genes (FvH4_c10g00010.1 and FvH4_c1g00300.1) that have not been assembled on the chromosome, the remaining 331 P/A genes were unevenly distributed among the seven chromosomes. The most P/A genes were located on the chromosome 3 (59), 4 (59) and 6 (62), respectively, and the fewest of 24 P/A genes were on the chromosome 7. Some P/A genes were discovered clustering near the telomeric regions (Fig. 3). In particular, fewer P/A genes were in the bottom telomeres than in the upper telomeres on the chromosomes 2, 3, 4 and 6. The P/A genes displayed no conspicuous gene cluster on chromosome 7. However, the P/A genes located in different sizes of gene clusters on the other six chromosomes. For instance, 13 P/A gene constituting a cluster were found on chromosomes 3 around 31–32Mb; a cluster containing 9 P/A genes was near 4–5Mb on chromosome 2.

Fig. 3.

The locations of the P/A genes on chromosomes. Gene IDs of the P/A genes were shown on the right of each chromosome, and the physical distance is listed on the left of each chromosome (Mb).

Number variations of the P/A genes in different species

F. virginiana and F. chiloensis were evolved from the four diploids of ancestral species (subgenomes A, B, C, and D) since about 1 million years (Njuguna et al. 2013; Hardigan et al. 2021). The allo-octoploid F. × ananassa species were originated from spontaneous hybrids between F. virginiana and F. chiloensis and underwent about 300 years of domestication and breeding selection (Edger et al. 2019). F. virginiana and F. chiloensis might be the intermediate genomes between the diploid and octoploid species and most progenitor genes from these two octoploid ancestral species were remained during the domestication process of Fragaria species (Edger et al. 2019; Feng et al. 2021). However, except F. vesca (A) and F. iinumae (B) were verified as the diploid ancestors of F. × ananassa, the other two subgenomes still remain unknown (Edger et al. 2020; Hardigan et al. 2021). Our study identified 333 candidate P/A genes in F. vesca Hawaii 4 and F. × ananassa Camarosa genomes, and the remaining 33 673 non-P/A genes were commonly distributed in Fragaria genus even across Rosaceae family. This finding was supported by massive highly conservative genes existed across the diploid Fragaria species (Hardigan et al. 2019; Feng et al. 2021). Although F. vesca genome contributed a large number of genes for the F. × ananassa genome, the F. × ananassa chromosomes also maintained more genes from the other ancestral diploids (Edger et al. 2019). The subgenome F. vesca accounts for dominant roles in the speciation of the allopolyploid cultivated strawberries, and follows its own independent evolutionary mechanism and stricter selection restrictions (Comai 2005; Edger et al. 2019). However, the diploid species have undergone widespread recombination in intermediate polyploids and a variety of changes have taken place during the evolution of the ancestral species (Edger et al. 2018; Hardigan et al. 2021). As for the 333 candidate P/A genes present in F. vesca and absent in F. × ananassa, the functions could be acquired by the other subgenomes.

The octoploid F. × ananassa contains four subgenomes, A, B, C, and D (Hardigan et al. 2021). F. vesca is considered to be the diploid ancestor of subgenome A of F. × ananassa and has made great contributions to the adaptation of cultivated strawberry to various environments (Edger et al. 2019; Edger et al. 2020; Hardigan et al. 2021). PAPs commonly appear in nature, but different species have different P/A gene numbers. Using the melon cultivar DHL92 as the reference genome, 1.5% of the genes in the analyzed varieties of C. melo were PAPs (Gonzalez et al. 2013). Out of the 10 255 single-copy genes surveyed in the stony coral (A. digitifera), 5% of them were determined to be P/A genes (Takahashi-Kariyazono et al. 2020). There was 2% of the genes exhibiting PAPs in 38 strains of Microbotryum lychnidis-dioicae, while only 0.6% of the autosomal genes were identified as P/A genes in 19 M. silenes-dioicae strains (Hartmann et al. 2018). Furthermore, among the 94 013 genes in B. napus genome, 38% of them exhibited PAPs (Hurgobin et al. 2018).

In our study, 333 genes were identified as PAPs between the two Fragaria species by using a strict criterion. The results demonstrated that approximately 0.98% of the F. vesca genes were absent from the F. × ananassa Camarosa genome. Only one gene was misidentified among the 32 randomly selected P/A genes checked for correctness by PCR. This false identification was due to a sequencing or assembly error in the whole-genome sequences of F. × ananassa. In previous studies, 75% of the P/A genes identified by PCR approaches in Arabidopsis were categorized correctly (Tan et al. 2012), and the accuracy of P/A gene identification by PCR amplification in B. napus was 86% (Gabur et al. 2020). The strict standards used here led to a higher identification accuracy in the studied strawberry species.

Gene duplications are significant for evolution and adaptation to environmental changes in organisms (Magadum et al. 2013). The emergence of P/A genes was closely related to gene duplications in many species. For example, a number of P/A multi-genes were caused by duplications in fungi (Hartmann et al. 2018), C. melon (Gonzalez et al. 2013) and Arabidopsis (Tan et al. 2012). In our results, 8.11% of the P/A genes were classified into multi-gene families, showing a relatively lower percentage of the P/A genes derived from duplications. This might be due to the small number of genome-wide duplication events in the woodland strawberry genome (Shulaev et al. 2011). Moreover, tandem duplications contributed to the formation of P/A gene clusters on chromosomes in the Arabidopsis genome (Kaul et al. 2000; Tan et al. 2012). In this study, 10.21% of the P/A genes were produced by tandem duplication, higher than the proportion of genes from WGDs and transposed duplications.

Distributions of the P/A genes on chromosomes

The genes with PAPs were commonly found to be non-randomly located on each chromosome, with a particular tendency to cluster at centromeres and/or telomeres. A large number of P/A genes are situated near centromeres in Arabidopsis (Tan et al. 2012) and at subtelomeric regions and centromeres in the fungal genome (Hartmann et al. 2018). P/A genes also occurred at the telomeric regions on the chromosomes (Winzeler et al. 2003). The locations of P/A genes between maize lines were not randomly distributed on chromosomes, and more P/A genes were located in the telomeric regions than in the centromeric regions (Darracq et al. 2018). The P/A genes tended to cluster in telomeres on the 3L chromosomal arms of D. melanogaster and D. simulans (Kern and Begun 2008). Similarly, the studied P/A genes in F. vesca were nonrandomly distributed along the seven chromosomes. Some of them tended to cluster near the telomeres, especially the chromosomes 2, 3, 4 and 6, with more P/A genes in the upper telomere regions than in the bottom telomeres.

The P/A genes in response to biotic and abiotic stresses

The GO analysis showed that the investigated P/A genes were involved in a variety of biological processes, molecular functions and cellular components. The P/A genes were relatively well-annotated for molecular functions and were especially involved in the binding of proteins, ATP and DNA. Moreover, among the biological process categories, most P/A genes participated in protein phosphorylation and oxidation–reduction processes (redox). Protein phosphorylation is involved in many cells physiological activities and actively participates in cell information transmission (Watanabe and Osada 2016). For example, protein phosphorylation regulates antioxidant enzyme activities to meet the challenges of waterlogging in maize (Xu et al. 2009). Phosphorylation of photorespiration enzymes is closely associated with a series of metabolic reactions in Arabidopsis (Hodges et al. 2013). Stimulating redox reactions enhances the response to biotic and abiotic stresses in plants (Cornic and Fresneau 2002; Gonzalez-Bosch 2018). Therefore, it could be inferred that some of the P/A genes might be related to adaptability to environmental stresses in strawberries.

The studied P/A genes encoding NB-ARC, PPR, PKs, F-box and EF-hand domains were functionally enriched for increased stress-related resistance (Kepinski and Leyser 2005; Wang et al. 2016; Zhu 2016; Zeng et al. 2017; Zhang et al. 2019). Similarly, some of the P/A genes were previously shown to belong to the plant resistance gene families (Bush et al. 2014; Rosa et al. 2015; Weisweiler et al. 2019; Gabur et al. 2020), such as NBS-LRR genes in Arabidopsis and C. melon (Shen et al. 2006; Tan et al. 2012; Gonzalez et al. 2013). In addition, PPR genes improved the resistance to cold stress in rice organs, and PKs enhanced drought tolerance in maize (Chang et al. 2017; Li et al. 2018a; Li et al. 2019a). F-box proteins can increase the abiotic stress resistance of peppers (Yu et al. 2007; Chen et al. 2014), and EF-hand proteins regulate calcium ions (Ca²⁺) to improve the plant defense and adaptability to environmental stresses in soybeans (Zeng et al. 2017). Therefore, it could be hypothesized that the P/A genes in strawberries might participate in the response to biotic and abiotic stresses.