text:
15 years of facts and figures on veterinary disciplinary measures in the Netherlands

In most countries, a veterinary disciplinary system is in force to ensure the quality of the veterinary profession and to offer an objective platform for complaints. We present an analysis of 15 years of veterinary disciplinary verdicts (2001–2016) to compare facts and figures and identify which factors are of major influence on the outcome of the verdicts. Rulings were collected from both paper files and the digital database of the veterinary disciplinary council (VDC), categorized, and used to create a database that enabled a statistical analysis. The results showed that complaints pertaining to companion animals are filed predominantly by owners, whereas complaints about livestock are mostly filed by the governmental civil servant (CS). CS complaints mostly address compliance issues. For the complaints made by owners (client complaints, CCs), reporting, communication, and veterinary mistakes appeared to be of statistical significance. Further studies are needed to investigate the impact of the complaints on veterinarians in general and how we can improve the veterinary disciplinary system.

Introduction
In most countries, like the Netherlands, a veterinary disciplinary system is in force to ensure the quality of the veterinary profession and to offer an objective platform for complaints. When first established in the Netherlands, its primary goal was animal welfare and food safety (1) but from the start, it served as a platform for client complaints (CCs) (2). Owners may perceive a veterinary disciplinary system as a body that they can use to file their complaints but it is actually a system to ensure quality and improve veterinary care (3,4). Owners rarely realize the effect of their complaints. CCs have a huge impact on veterinarians and their co-workers (4)(5)(6). Recent studies investigating the effect of complaints in both human and veterinarian medicine have shown that complaints contribute to psychological distress such as depression, burnout, anxiety, and even suicidal ideation (7, 8), and health problems such as cardiovascular disease, .
/fvets. . insomnia, headaches, anger, and relationship problems (5,9). Although the possibility of complaining may stimulate health workers to improve their level of care, the detrimental effect on job satisfaction and mental health (4-6) may warrant a critical analysis of these disciplinary systems as well. The Dutch veterinary disciplinary system was first established in 1990. In contrast to the disciplinary system of, for example, the UK (10) or the USA (4), where complaints are first evaluated by a committee and only taken into consideration after careful scrutinizing, nearly all complaints sent to the Dutch veterinary disciplinary council (VDC) are taken into consideration. The original intent of the Dutch government, when establishing the VDC, had been the importance of veterinary care for livestock as it is one of the most important agricultural economic activities in Dutch (11). For this reason, the legal system was not designed to separate complaints in different ways. Regardless of the nature of the complaint, veterinarians feel forced to respond. This may have, like in other countries, a huge impact on their well-being (4-6).   Although VDC and VAC are in function for over 29 years, a scientific analysis of the nature of the complaints and the rulings of the VDC and VAC has not yet been performed. This study aimed to describe the sources of complaints, the reasons for the complaints, and the kind of rulings made.

Materials and methods
For the analysis in this study, we chose the VDC rulings from the years 2001-2016. Rulings prior to this date were either incomplete or not available for this analysis and hence had to be excluded. At the time of analysis of this study, most of the rulings beyond this time frame were, for various reasons, not yet available to us and will therefore be analyzed in a different study.
This study was neither subjected to ethical approval as it does not involve any use or experimental use of animals nor will any data be disclosed that are subject to any privacy laws.
The rulings were systematically arranged in an inventory data set. The data were collected both from paper files and the digital database of the VDC. Only the data of complaints handled by the VDC and VAC are included in this Excel database. Complaints that have not been handled, or resolved for various reasons, such as administrative errors, have not been registered and are therefore unknown to us.
All rulings contained data that were used to establish different value labels which could be categorized into excel rows. Different value labels were created such as animal species, size, gender, age, type of complaint, outcome, etc. A complete list of these different value labels and their descriptions are found in Table 1. The value labels were coded into subcategories to enable further analysis. A complete list of these value labels and subcategories is listed in Appendices I and II.
The complete data set consisted of rulings of the VDC and VAC. The data sets were first checked for errors and missing values, and when possible, assigned a numeric property and imported into SPSS. All rulings, whether from the VDC or VAC, used the same value labels. The VDC rulings were grouped into six categories: (1) the complaint is either dismissed or inadmissible for instance due to administrative reasons, (2) unfound, (3) found with either no measure or a warning or reprimand, (4) a fine conditional or unconditional, (5) Microsoft Excel .

SPSS version for statistical analyses.
Frontiers in Veterinary Science frontiersin.org . /fvets. . a fine and/or (un)conditional suspension, and (6) partial or complete removal of the veterinary profession. The VAC data are measured in equal, heavier, or lower rulings, and an additional parameter was the person who appealed to the ruling of the VDC, such as either the CS, veterinary, or the owner. The veterinary disciplinary council data set was first analyzed using all rulings. First, the effect of time on the outcome was studied: found or unfound, and second, whether there was a time effect for the ruling itself. We used Kendall's tau association test because the rulings were non-parametric ordinal data.
As the complaints filed by the CS address mostly legal issues and hence are of a completely different nature compared to those filed by owners, the CS and CC rulings were analyzed separately. Each subcategory (Appendices I and II) was analyzed for the two types of rulings: CS and CC rulings. To see whether there was a statistically significant effect for each subcategory for either "found" or "unfound, " a Fisher Exact test with a pvalue for the significance of <0.05 was used. The CS rulings were analyzed for (1) animal species, (2) status (deceased or not), (3) influence of reporting, and (4) the summary of the complaint. If the information was missing (for instance whether the animal had died or not), that ruling was excluded from the analysis.
The CC rulings were analyzed for (1) animal species, (2) sex, (3) age, (4) status (deceased or not), (5) influence of reporting, (6) dog breed, and (7) summary of the complaint. If the information was missing (for instance age or breed), that ruling was excluded from the analysis.
The veterinary disciplinary appeal council rulings were analyzed for (1) the person who appealed, (2) the type of appeal, and (3) the influence of reporting with as a dependent parameter the new ruling, using a Fisher Exact test with a pvalue for significance <0.05. If applicable, a two-tailed Pearson correlation coefficient was calculated.

Descriptive analysis rulings VDC
A total of 1,197 verdicts were available for analysis. The number of complaints involving companion animals was 890 (74%) out of 1,197 verdicts; 157 (13%) complaints concerned equines, and 150 (13%) were related to livestock. In total, 1,023 complaints were filed by owners (CCs) (85%) representing a vast majority in contrast with the 174 cases filed by the CS (15%). The total number of cases did not vary greatly during this period. The mean number was 79 cases per year with a range of 68 to 92 Frontiers in Veterinary Science frontiersin.org . /fvets. . ( Figure 1). During the analyzed period, the number of CS cases gradually increased and the number of CC cases decreased. In both groups, the number of inadmissible cases declined during the analyzed period (Figures 2, 3). Cases filed by the CS addressed more than 96% of compliance cases (Table 2). During the analyzed period, the number of CS cases gradually increased. During the period from 2008 to 2011, more rulings were made with (un)conditional fines and suspension (Figure 2). After this period, most rulings concerned fines (Figure 2), and this difference was statistically significant (p < 0.001). A similar statistically significant effect was observed for the year, irrespective of whether the ruling was found or unfound (p = 0.011). Cases filed in the period after 2008 were more likely to be found. The majority of CS cases concerned livestock (79%), 8% concerned companion animals, and 13% concerned equines (Table 2). There was no statistically significant difference for animal species, status, or influence of reporting, irrespective of whether the ruling was unfound or found (Table 2). However, there was a clear statistically significant difference if the summary of the case addressed compliance issues (p = 0.01).
The number of CCs gradually decreased as mentioned above ( Figure 3 and Table 4). This was partly due to the decline in inadmissible cases. There was a clear effect for year and ruling (p < 0.001). In the first analyzed period, the number of cases that were considered to be unfound decreased, whereas the number of cases that were considered to be found remained constant ( Table 3). The CC cases were analyzed for several categories, irrespective of whether they were judged to be found or not. Animal species, sex, age, status (deceased or not), or dog breed had no statistically significant effect on the outcome of the ruling (Table 4). However, there was a clear statistically significant effect for reporting and the summary of the complaint (mistakes and communication) (Table 4). Poor reporting resulted more likely in a "found" ruling (p < 0.001). The statistically significant effect was even greater for communication and veterinary mistakes (p < 0.001).

Descriptive analysis rulings veterinary appeal council
A total of 177 appeals were submitted during the analyzed period (Table 5). In the majority of these cases, the VAC came to a similar verdict as the VDC (70%). In 18% of the cases, the verdict was lower compared to the VDC and higher in the remaining 12%. Owners, veterinarians, and the CS could appeal, where owners appealed in 53% of the cases and .
/fvets. .  veterinarians appealed in 42%. The variables appellant, type of appeal complaint, and influence of reporting proved to be of statistically significant influence (Table 5).

Discussion
This study is an analysis of 15 years of rulings made by the VDC and VAC. The goal of these rulings is 2-folds: to improve veterinary care and to serve as an objective platform for complaints (2, 3). However, the majority of cases concerned CCs from companion animal owners. Only around 12% of the cases concerned livestock animals, and the majority of the cases were filed by the CS. These CS cases mostly addressed compliance issues and if the veterinarians did not fulfill their duty of care, they are at great risk of being fined, or temporally suspended as a veterinarian. Based on the outcome of the CS cases, we concluded that if it concerns compliance issues, the system works. However, as the number of cases only increased during this period, either the CS became more active in controlling or a learning effect is absent. Although it was not the intent of the Dutch government to create a platform, especially for CCs (2), the vast majority of cases concerned CCs. CCs have a huge effect on the well-being of veterinarians and co-workers (4-6), .
/fvets. .  Although it should be possible to complain as it may also improve our level of care (6), the detrimental effect on job satisfaction and mental health (4, 5) justified a critical analysis of the rulings.
Looking at the performed analysis, there are, in light of the effect of a CC, some recommendations. In principle, every CC is taken into consideration by the VDC. The legal system does not provide an option to treat complaints of a less heavy character, and/or rooted in the owner's frustration, through a different complaint route. The "unfound" ruling can only be made after a full disciplinary procedure, which plays a major role in the veterinarian's well-being and job satisfaction (4-6). It seems safe to conclude that the current method of working could be adapted.
This study has limitations. The actual impact of a CC on the veterinarian could not be studied. It is not possible to establish the exact amount of learning effect. It would have been very interesting to compare the outcome of this study with different veterinary disciplinary systems in other countries. There is, one recent Californian study providing insights from the veterinary disciplinary system, but this study addresses different issues (4). But all of these questions will be subject to additional studies.

Conclusion
It is concluded that compliance, performing its duty of care, proper record keeping, good communication, and avoiding veterinary mistakes can reduce the number of CS complaints and CCs. Additional studies are needed to investigate the impact and learning effect of CCs on veterinarians.

Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.
. /fvets. . of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.


provenance: AgriculturalAndFoodSciences,2022-2022-Agricultural And Food Sciences-2022.jsonl.gz.gz

text:
HuNAC20 and HuNAC25, Two Novel NAC Genes from Pitaya, Confer Cold Tolerance in Transgenic Arabidopsis

NAC transcription factors are one of the largest families of transcriptional regulators in plants, and members of the gene family play vital roles in regulating plant growth and development processes including biotic/abiotic stress responses. However, little information is available about the NAC family in pitaya. In this study, we conducted a genome-wide analysis and a total of 64 NACs (named HuNAC1-HuNAC64) were identified in pitaya (Hylocereus). These genes were grouped into fifteen subgroups with diversities in gene proportions, exon–intron structures, and conserved motifs. Genome mapping analysis revealed that HuNAC genes were unevenly scattered on all eleven chromosomes. Synteny analysis indicated that the segmental duplication events played key roles in the expansion of the pitaya NAC gene family. Expression levels of these HuNAC genes were analyzed under cold treatments using qRT-PCR. Four HuNAC genes, i.e., HuNAC7, HuNAC20, HuNAC25, and HuNAC30, were highly induced by cold stress. HuNAC7, HuNAC20, HuNAC25, and HuNAC30 were localized exclusively in the nucleus. HuNAC20, HuNAC25, and HuNAC30 were transcriptional activators while HuNAC7 was a transcriptional repressor. Overexpression of HuNAC20 and HuNAC25 in Arabidopsis thaliana significantly enhanced tolerance to cold stress through decreasing ion leakage, malondialdehyde (MDA), and H2O2 and O2− accumulation, accompanied by upregulating the expression of cold-responsive genes (AtRD29A, AtCOR15A, AtCOR47, and AtKIN1). This study presents comprehensive information on the understanding of the NAC gene family and provides candidate genes to breed new pitaya cultivars with tolerance to cold conditions through genetic transformation.

Introduction
Transcription factors (TFs) and cis-elements function in the promoter region of different stress-related genes, and thereby alter their expression in response to the stress tolerance [1,2]. Many plant TFs, including MYB, bHLH, AP2, MYC, WRKY, and NAC, have been identified and play diverse functions in various biological processes [3]. Among those TFs, NACs are one of the largest families and play important roles in diverse developmental processes in plants [4,5]. The name of the NAC gene family was derived from the three earliest characterized proteins with a particular domain (NAC domain) from petunia NAM (no apical meristem), Arabidopsis ATAF1/2, and CUC2 (cup-shaped cotyledon) [6,7].
Protein sequences of this family reveal that a typical NAC TF has a highly conserved NAC domain with about 160 amino acid residues at the N-terminal region while the C-terminal region is highly diversified in length and sequence, which is considered the transcriptional activation domain [8]. The NAC domain is further divided into five subdomains (A-E) that represent motifs for both DNA-binding and protein-protein interactions [9].

Genome-Wide Identification of NAC Family Genes
To identify pitaya NAC TF encoding genes, all proteins were annotated from the H. undatus genome [36]. In total, 64 NAC genes were identified and named HuNAC1 to HuNAC64 according to their chromosomal position. Each of these HuNAC proteins contained an NAM domain (PF02365. 15), a specific conserved domain of NAC TF protein family (Supplementary File S1). These genes encoded predicted proteins ranging from 154 to 684 AA (amino acids) with isoelectric point (PI) values ranging from 5.33 to 9.19 and molecular weights from 17.95 to 75.87 KDa (Table S1). Subcellular location of these genes was predicted using an online tool from Molecular Bioinformatics Center (http://cello.life.nctu.edu.tw/) (Accessed on 13 January 2022). The subcellular localization and the protein sequences of the 64 HuNACs were listed in Supplementary File S2. Among the 64 NAC proteins, five (HuNAC16, HuNAC17, HuNAC21, HuNAC52, and HuNAC55) were predicted to be located in the cytoplasmic; three (HuNAC3, HuNAC34, and HuNAC48) were chloroplast; two (HuNAC4 and HuNAC10) were mitochondrial; and the rest were localized in the nucleus.

Phylogenetic Analyses of NAC Family Members
To explore the evolutionary relationships among HuNAC genes, a Maximum Likelihood (ML) phylogenetic tree was constructed according to NAC protein sequences from H. undatus and A. thaliana (Figure 1). Based on the ANAC classification and NAC domain alignments of HuNACs, all members of the HuNACs and ANACs were divided into two major groups: Groups A and B. The tree was divided into fifteen subgroups in Group A (A1-15) according to similarities in NAC domain structures. The number of HuNAC genes in each subgroup varied greatly. Subgroups A14 and A13 had only one and two genes, respectively. Subgroup A1 constituted the largest clades with 12 pitaya NAC members, followed by Subgroup A3 with 10 sequences. However, no members were detected in Group B in pitaya, suggesting that they might have been lost in these subfamilies. Phylogenetic tree analyses of NAC transcription factor family in pitaya (64 genes) and A. thaliana (105 genes). The full-length sequences of the NAC proteins were aligned using ClustalW, and the phylogenetic tree was constructed using the maximum-likelihood method in the MEGAX software. The Bootstrap value was 1000 replicates. HuNACs were indicated by blue stars. All members from both species were designated as Group A (A1-A15) and Group B. Each of the subfamily is indicated in a specific color.

Gene Structure and Conserved Motif Analyses of HuNACs
To better understand the similarity and diversity of the HuNAC genes, the exon/intron distributions of HuNAC and conserved motifs were analyzed. As shown in Figure 2, the 64 NAC genes were divided into 15 subcases. Genes within the same subcase had a similar exon/intron structure in terms of intron number and exon length. For example, among the 12 members in subgroup A3, 11 had two introns, while six genes in subgroup A4 had 2-5 introns. The conserved motifs of NAC family proteins in pitaya were investigated using MEME online software. Based on this program, 20 distinct motifs were identified ( Figure S1). According to frequencies of occurrence, motifs 1, 4, and 3 were the three most frequently presented motifs which were observed in all HuNAC proteins, and most of the conserved motifs were located at the N-terminus. Among the 20 motifs, motif 2, motif 4, motif 3, and motif 5 had the subdomains A, B, D, and E, respectively, and motifs 1 and 6 contained subdomain C ( Figure 2D). These results are consistent with the finding that a connection existed between subfamilies and motifs [37].

Chromosomal Localization and Synteny Analyses of HuNACs
Genome chromosomal location analyses revealed that the 64 HuNAC genes were unevenly scattered on all 11 chromosomes (Figure 3). Chromosomes 1, 3, 4, and 8 contained 12, 8, 8, and 8 HuNAC genes, respectively. Chromosomes 6, 10, and 11 had only three HuNAC genes. The average numbers of HuNAC genes on the chromosome were approximately 6.0. The phenomenon of gene duplication has been recognized to occur throughout plant evolution, and plays an important role in the expansion of the large gene families in plant [38]. Gene duplication events were investigated to clarify the expansion patterns of the NAC family in pitaya. Twenty segmental duplicated events with 33 HuNAC genes were identified in the pitaya genome ( Figure 4; Supplementary File S3). The high-identity segmental duplication events in pitaya suggested that this duplication type likely plays a crucial role in the expansion of the pitaya NAC gene family.

Expression Analyses of HuNAC Genes under Cold Treatment
In plants, many NAC genes are involved in the response to abiotic stresses. Therefore, the expression patterns and putative functions of all 64 HuNAC genes were analyzed under cold conditions, and 49 HuNAC genes were induced under cold treatment ( Figure S2). As shown in Figure 5, compared with normal temperature, higher expression levels of HuNAC7, HuNAC20, HuNAC25, and HuNAC30 were detected at a low temperature. The expression of HuNAC7, HuNAC20, HuNAC25, and HuNAC30 showed a trend of gradual increase. HuNAC7, HuNAC25, and HuNAC30 were strongly induced at 24 h after cold treatment, while expression levels of HuNAC20 significantly increased at 36 h after cold treatment. Expressions of HuNAC7, HuNAC20, HuNAC25, and HuNAC30 reached their maximum levels at 60 h after cold treatment.

Subcellular Localization and Transcriptional Activation Analyses of HuNACs
To analyze the subcellular localization of HuNAC7, HuNAC20, HuNAC25, and HuNAC30, their full-length coding sequences were fused with the GFP to construct 35S-HuNACs-GFP vectors, and transiently expressed in leaves of N. benthamiana. As shown in Figure 6, the fluorescence of HuNACs-GFP were predominately observed in the nucleus of epidermal cells, while the GFP signal of positive control was detected around the cytoplasm and the nucleus. Those results suggested that the HuNAC7, HuNAC20, HuNAC25, and HuNAC30, like the other reported NAC proteins [39], encoded nuclear proteins. Full-length coding regions of HuNAC7, HuNAC20, HuNAC25, and HuNAC30 were fused to the GAL4BD to generate pGBKT7-HuNACs fusion plasmids to study the transcriptional activation abilities of HuNACs. Like the positive control (pGBKT7-53+pGADT7-T), yeast cells expressing the HuNAC20, HuNAC25, and HuNAC30 grew well on SD/-Trp-His-Ade and showed α-galactosidase activity, indicating that HuNAC20, HuNAC25, and HuNAC30 had transactivation ability in yeast cells; however, yeast cells expressing the HuNAC7 could not grow on SD/-Trp-His-Ade (pGBKT7) ( Figure 7A). Trans-activation of HuNACs were further verified in leaves of N. benthamiana using the dual-luciferase reporter system. Compared with the negative control, the expression of the positive control (BD-62SK-VP16), HuNAC20, HuNAC25, and HuNAC30 resulted in a higher value of the LUC/REN ratio. However, HuNAC7 significantly repressed the expression of the LUC reporter in comparison to the negative control (pBD) and the ratio of LUC/REN of HuNAC7 was 0.3-fold compared with the pBD ( Figure 7C). These results demonstrated that HuNAC20, HuNAC25, and HuNAC30 are transcriptional activators while HuNAC7 is a transcriptional repressor.

Phylogenetic and Sequence Analyses of HuNAC20 and HuNAC25
HuNAC20 and HuNAC25 were cloned based on the sequences of the pitaya genome [36]. The full-length coding DNA sequences (CDSs) of HuNAC20 and HuNAC25 were both 891 bp, encoding 296 amino acids with molecular weights of 33.74 and 33.65 kDa and pIs of 7.69 and 6.26, respectively. Homologous analyses of HuNAC20, HuNAC25, and 18 NAC protein sequences from different plant species were aligned with DNAMAN software. The similarity between HuNAC20 and HuNAC25 was 17.7% ( Figure S3a). Both HuNAC20 and HuNAC25 had a conserved subdomains (A-E) in the N-terminal region and a diverse activation domain in C-terminal, respectively. HuNAC20 and HuNAC25 both contained a conserved nuclear localization signal (NLS) (Figure S3b,c). HuNAC20 shared 78.2% similarity with CqNAC2 (XP_021772603.1) from Chenopodium quinoa and 77.6% similarity with BvATAF2 (QGZ00533.1) from Beta vulgaris. HuNAC25 had 67.2% similarity with BvNAC (XP_010679326.1) from Beta vulgaris and 64.0% similarity with MeNAC (XP_021630142.1) from Manihot esculenta. The trans-activation ability of HuNACs is indicated by the ratio of LUC to REN. The LUC/REN ratio of the empty pBD vector (negative control) was used as a calibrator (set as 1). pBD-VP16 was used as a positive control. The asterisk indicates a significant difference at the 1% level compared to the pBD.

Overexpression of HuNAC20 and HuNAC25 in Arabidopsis Enhanced Tolerance to Freezing Stress
To further explore the function of HuNAC20 and HuNAC25, we generated transgenic A. thaliana plants constitutively expressing HuNAC20 and HuNAC25 driven by the 35S promoter. Homozygous T 3 lines were obtained on the basis of 3:1 segregation for kanamycin resistance phenotype. Expression levels of the target genes in the homozygous transgenic lines were analyzed by RT-qPCR. Two independent T 3 transgenic lines with a relatively high expression of HuNAC20 (NAC20-L1 and NAC20-L3) and HuNAC25 (NAC25-L2 and NAC25-L3) were selected for cold tolerance experiments.
At normal conditions, no significant difference in morphology was observed between transgenic lines and WT plants. After cold acclimation (48 h at 4 • C) and recovery for 6 d at 22 • C, transgenic plants displayed better performance than wild-type under freezing treatment (−6 • C for 6 h) ( Figures 8A and 9A). Most WT plants died with a survival rate at around 8.0%. The survival rates of the transgenic lines (90.0% for NAC20-L1, 88.0% for NAC20-L3, 82.0% for NAC25-L2 and 78.0% for NAC25-L3) were significantly higher than that of WT plants (8.0%) (Figures 8C and 9C). These results indicated that overexpression of HuNAC20 and HuNAC25 in A. thaliana enhances tolerance to cold stress. − contents in the transgenic lines were significantly lower than those in WT plants ( Figure 10). These results suggested that overexpression of HuNAC20 and HuNAC25 enhances A. thaliana tolerance to cold stress by altering the ion leakage, MDA, H 2 O 2 , and O 2 − accumulation. . Asterisks indicate significant differences (** p < 0.01) between the transgenic lines and WT plants.

Overexpression of HuNAC20 and HuNAC25 Activated the Expression of Cold-Responsive Genes under Cold Stress
The transcript levels of cold-responsive genes including AtRD29A, AtCOR15A, AtCOR47, and AtKIN1 were analyzed in HuNAC20 and HuNAC25 transgenic lines and WT plants under cold treatment. The expressions of these genes in transgenic lines and WT plants were relatively low under normal conditions (22 • C). Compared with WT plants, higher transcripts of AtRD29A, AtCOR15A, AtCOR47, and AtKIN1 were observed in HuNAC20 and HuNAC25 transgenic plants after 24 h under cold treatment ( Figure 11). These results indicated that HuNAC20 and HuNAC25 could activate expression levels of AtRD29A, AtCOR15A, AtCOR47, and AtKIN1 responsible for stronger tolerance to cold stress in transgenic A. thaliana lines.

Discussion
NAC is one of the largest families of TFs unique to plants and plays important roles in defense against harsh environmental conditions [40,41]. NAC TFs also play crucial roles in plant physiological processes [42]. With the development of whole genome sequencing technology, great progress has been made in identification of NAC genes from Arabidopsis [8], rice [8], cucumber [24], pepper [32], and apple [34]. However, little information is available about the NAC family in pitaya. In the present study, a total of 64 NAC genes distributed in 11 chromosomes were identified from pitaya ( Figure 3). Great differences in the size, sequence, physical, and chemical properties of the proteins encoded by HuNACs were detected, which was consistent with NACs from the other plant species [22,24,33,38]. The phylogenetic analyses showed that these HuNACs were classified into 15 subgroups with ANACs in the Group A, and the subgroup A1 had the most HuNACs (19%), followed by the A3 (16%), whereas the A14 subgroup had the fewest genes (2%) (Figure 1). However, Group B had no member of HuNACs similar to the results from quinoa [43]. Most HuNACs had three exons and two introns (Figure 2), indicating that the genetic structural diversity of NAC in pitaya, which was similar to that of the other species such as kiwifruit [44], Musa Acuminata [45], cucumber [24], and white pear [22]. High segmental duplication events resulted in the expansion of the NAC gene family in H. undatus (Figure 4), which was consistent with the evolutionary pattern in pineapple [38]. Therefore, we speculated that these pitaya-specific NAC genes may have special roles in pitaya.
The growth and development of pitaya are frequently affected by abiotic stresses such as drought, high salinity, or extreme temperatures [46,47]. Cold stress is one of the most important influencing factors restricting pitaya production in South China. Previous studies have suggested that NAC TFs can be induced by cold stress. Overexpression of SlNAC35 enhanced chilling tolerance of transgenic tomato [48]. Compared with non-cold control fruits, the expression level of MaNAC1 in the fruits directly stored at 7 • C was significantly increased [49]. VvNAC17, a novel NAC TF, was expressed in various tissues under cold treatment and enhanced freezing tolerance in transgenic Arabidopsis [50]. In this study, expressions of HuNAC7, HuNAC20, HuNAC25, and HuNAC30 were strongly induced under cold-stress ( Figure 5), indicating that they may be involved in plant responses to cold stress. Subcellular localization is important to elucidate protein function. In our study, HuNAC7, HuNAC20, HuNAC25, and HuNAC30 were nuclear proteins (Figure 6), which is consistent with the other reported NAC proteins [16,39,45].
NAC proteins are involved in plant responses to cold stress. Overexpressing LlNAC2 in A. thaliana enhanced tolerance to cold stress. Compared with WT plants, the transgenic plants had lower electrolyte leakage and higher levels of soluble sugars [28]. Overexpression of GmNAC20 confers tolerance of freezing stress in transgenic plants by regulating downstream stress-responsive genes such as cor6.6, cor78, and cor15A [14]. Overexpressing MbNAC25 in Arabidopsis increased cold tolerance via enhanced scavenging capability of reactive oxygen species (ROS) under low-temperature stress (4 • C) [51]. In the present study, two novel NAC TFs, i.e., HuNAC20 and HuNAC25, were cloned from pitaya. HuNAC20 and HuNAC25 were nuclear-localized NAC-type DNA-binding proteins. Overexpression of HuNAC20 and HuNAC25 in Arabidopsis enhanced cold tolerance. Ion leakage, MDA contents, and H 2 O 2 and O 2 − accumulation are commonly used to assess the severity of membrane lipid damage under stress conditions [52]. In our study, overexpression of HuNAC20 and HuNAC25 in Arabidopsis lines resulted in lower contents of ion leakage, MDA, and H 2 O 2 and O 2 − under cold stress. Compared with WT, higher expression levels of the cold-responsive genes (AtRD29A, AtCOR15A, AtCOR47, and AtKIN1) were detected in the transgenic plants under cold stress ( Figure 11). These results indicated that overexpression of HuNAC20 and HuNAC25 in A. thaliana could activate expressions of AtRD29A, AtCOR15A, AtCOR47, and AtKIN1 under cold stress ( Figure 12). The results are consistent with those NAC TFs from the other plants species such as Lilium lancifolium [28], Solanum lycopersicum [48,53], and Malus baccata [51]. Putative binding cis-elements of NAC proteins have been identified to localize in the promoter sequences of these genes [54]. However, the mechanism of interaction between HuNAC20 and HuNAC25 and the cold-responsive genes is still not clear yet, and further study is necessary to elucidate it.

Conclusions
In summary, this study provides the first report on identification and characterization of the NAC gene family based on the genome-wide analyses of the H. undatus genome. A total of 64 HuNAC were identified and divided into fifteen subfamilies. The 64 HuNAC were unevenly scattered on all 11 chromosomes. HuNAC7, HuNAC20, HuNAC25, and HuNAC30 may be involved in cold stress tolerance according to their expression patterns. HuNAC7, HuNAC20, HuNAC25, and HuNAC30 were nucleus proteins. HuNAC20, HuNAC25, and HuNAC30 had transcriptional activity while HuNAC7 acted as a transcriptional repressor. Overexpressing HuNAC20 and HuNAC25 in Arabidopsis enhanced tolerance of cold stress by altering the expression of cold-responsive genes in transgenic plants. The results of the present study provide valuable information for a better understanding of NAC TFs involved in cold-stress response in pitaya.

Identification of the Pitaya NAC Gene
The Hidden Markov Model (HMM) profiles of the NAM domain PF02365.15 were downloaded from the Pfam database (https://pfam.xfam.org/) (Accessed on 13 January 2022). HMM was used to search NAM (PF02365.15) domains from the pitaya genome with values (e-value) cut-off at 1.0 [36]. The integrity of the NAM domain was determined using the online program SMART with an e-value < 0.1. The length, molecular weight, and isoelectric point of each NAC protein were predicted by the online ExPasy program (https://web.expasy.org/protparam/) (accessed on 13 January 2022).

Phylogenetic Analyses of the Pitaya NAC Gene
To investigate the phylogenetic relationship of the NAC gene families in Arabidopsis and H. undatus, the Arabidopsis NAC protein sequences were downloaded from the Arabidopsis Information Resource (https://www.arabidopsis.org/) (accessed on 13 January 2022). All NAC TFs were aligned using the MUSCLE set at default parameters [55]. The Maximum Likelihood (ML) phylogenetic tree was constructed using the MEGA7 program [56]. Bootstrapping was performed with 1000 replications. The tree was exhibited by the EVOLVIEW online tool (https://www.evolgenius.info/evolview/) (accessed on 13 January 2022).

Analyses of Gene Structure and Conserved Motifs in the HuNAC Family
The gene structure of each HuNAC gene was drawn using the TBtools software [57]. The MEME (http://meme-suite.org/tools/meme) (accessed on 13 January 2022) was used to identify the unknown conserved motifs with the following parameters: site distribution: zero or one occurrence (of a contributing motif site) per sequence, maximum number of motifs: 25, and optimum motif width ≥6 and ≤200.

Chromosomal Locations and Synteny Analyses of HuNAC Genes
The chromosome location of each HuNAC gene was obtained from the pitaya genome [36]. These data were then integrated and plotted using Mapchart software [58]. The Multiple Collinearity Scan toolkit (MCscanX) was applied to analyze the duplication pattern for each HuNAC followed the operation manual [59].

Plant Material
Stems (6-8 cm in height) in vitro from pitaya cultivar 'Hongguan No. 1' (Hylocereus monacanthus) on rooting medium (MS containing 0.5 mM IBA) for 15 d were used as materials [60]. The pitaya plants were cultured in a climate chamber at the South China Agricultural University (Guangzhou, China) under a temperature range of 23-25 • C and a 16 h light (50 µmol m −2 s −1 ) and 8 h dark. Plantlets were placed in a climate cabinet at 5 • C to analyze the expression patterns of the NAC family. The stems from treatments and control were collected respectively at 12, 24, 36, 48, and 60 h after treatment, immediately frozen in liquid nitrogen and stored at −80 • C until future analysis.
Arabidopsis thaliana Columbia-0 (Col-0) was used for genetic transformation of HuNAC20 and HuNAC25. Transgenic lines and WT plants were grown in 9 cm × 9 cm plastic pots containing a 1:1 mixture of sterile peat soil and vermiculite with normal man-agement in a growth chamber [ Three-week-old seedlings of transgenic Arabidopsis lines and WT were treated at −6 • C for 6 h after cold acclimation at 4 • C for 48 h, followed by 4 • C in the dark for 12 h. The plants were transferred to normal conditions (22 • C) for recovery for 6 d, and then the survival rates were counted. Photos were taken before freezing and after recovery. The rosette leaves were collected for analyses of ion leakage [61], malondialdehyde (MDA) content, and H 2 O 2 and O 2 − accumulation [62,63]. For chilling treatment, the seedlings were put in a 4 • C incubator, and leaves were collected respectively at 0, 6, 12, and 24 h for expression analyses of cold-responsive genes (AtRD29A, AtCOR15A, AtCOR47, and AtKIN1) [64], and AtACTIN2 (AT1G13320) were used as internal control [61]. Specific primers are listed in Supplementary File S4.

Gene Cloning and Expression Analyses
Total RNA was isolated using the EASYspin Plus Complex Plant RNA Kit (RN53) (Aidlab Biotechnology, Beijing) according to the manufacturer's protocol. Single-stranded cDNA was synthesized using the PrimeScript™ RT Reagent Kit with gDNA Eraser (TaKaRa, Shiga, Japan). The RT-qPCR primers were designed by BatchPrimer3 (https://probes.pw. usda.gov/cgi-bin/batchprimer3/batchprimer3.cgi) (accessed on 13 January 2022), and the Actin(1) reference gene was used as the internal control [65]. RT-qPCR was performed in an CFX384-Real-Time system (C1000 Touch Thermal Cycler, Bio-Rad, CA, USA) using the Re-alUniversal Color PreMix (SYBR Green) (TIANGEN, Beijing, China). Specific primers are in Supplementary File S4. Each experiment was repeated in triplicate using independent RNA samples. Relative gene expression levels were calculated using the 2 −∆∆C T method [66].
The full-length coding sequences of HuNACs were cloned using I-5TM2× High-Fidelity Master Mix (MCLAB, San Francisco, CA, USA) with specific primers (Supplementary File S5).

Subcellular Localization Analyses
The full lengths of HuNAC7/20/25/30 without a stop codon were subcloned into the pGreen-35S-GFP vector to fuse with the gene sequence of green fluorescent protein (GFP) (primers are listed in Supplementary File S5). Then, the pGreen-HuNAC7/20/25/30-35S-GFP and the control pGreen-35S-GFP vector were transferred into the Agrobacterium tumefaciens strain GV3101 (pSoup-p19), and injected into the abaxial side of 4-to 6-week-old N. benthamiana leaves. After 48 h of infiltration, infected leaf tissues were collected for analyses. The GFP signal was captured under a fluorescence microscope (ZEISS LCM-800, Oberkochen, Germany). All assays were repeated three times.

Transcriptional Activation Analyses in Yeast Cells
The pGBKT7-HuNACs, pGBKT7-p53 and pGBKT7 empty plasmids were transferred into the Y2HGold yeast strain independently using the lithium acetate method (PT1172-1, Clontech) (primers are listed in Supplementary File S5). The transformed yeast cells were cultured on SD/-Trp and SD/-Trp-His-Ade medium. The growth status of yeast cells and the activity of X-α-galactosidase (X-α-Gal) were observed after incubation with 20 mg/mL X-α-Gal for 10-30 min.

Dual-Luciferase Reporter Assays in N. benthamiana Leaves
For transcriptional activity analyses of HuNACs in N. benthamiana leaves, coding sequences of HuNAC7/20/25/30 were cloned into the 35S promoter-driven pBD vectors to fuse with the yeast GAL4 DNA-binding domain (GAL4BD) as an effector (pBD-HuNACs). The double-reporter vector contained a firefly luciferase (LUC) driven by five copies of the GAL4-binding elements (5×GAL4) and minimal TATA region of CaMV 35S. Renilla luciferase (REN) in the same vector driven by the CaMV 35S promoter was used for normalization. The primers are listed in Supplementary File S5.

Arabidopsis thaliana Transformation and Phenotypic Analyses
The open read frame (ORF) of HuNAC20 and HuNAC25 were subcloned into the pPZP6K90 vector under the control of the 35S promoter and introduced into Agrobacterium tumefaciens strain GV3101. The recombinant vectors were transformed into Arabidopsis Col-0 using the floral dip method. Positive plants were screened on MS medium containing 100 mg L −1 kanamycin and identified by PCR detection. The expression levels of HuNAC20 and HuNAC25 in the homozygous T 3 transgenic lines were analyzed by RT-qPCR, and the primers used were listed in Supplementary File S5. Two homozygous T 3 transgenic lines were used for cold tolerance experiments and photographed with a digital camera (G16, Canon, City, Japan).


s2fieldsofstudy:
Agricultural And Food Sciences
Biology

text:
Dissection of the Genetic Basis of Yield Traits in Line per se and Testcross Populations and Identification of Candidate Genes for Hybrid Performance in Maize

Dissecting the genetic basis of yield traits in hybrid populations and identifying the candidate genes are important for molecular crop breeding. In this study, a BC1F3:4 population, the line per se (LPS) population, was constructed by using elite inbred lines Zheng58 and PH4CV as the parental lines. The population was genotyped with 55,000 SNPs and testcrossed to Chang7-2 and PH6WC (two testers) to construct two testcross (TC) populations. The three populations were evaluated for hundred kernel weight (HKW) and yield per plant (YPP) in multiple environments. Marker–trait association analysis (MTA) identified 24 to 151 significant SNPs in the three populations. Comparison of the significant SNPs identified common and specific quantitative trait locus/loci (QTL) in the LPS and TC populations. Genetic feature analysis of these significant SNPs proved that these SNPs were associated with the tested traits and could be used to predict trait performance of both LPS and TC populations. RNA-seq analysis was performed using maize hybrid varieties and their parental lines, and differentially expressed genes (DEGs) between hybrid varieties and parental lines were identified. Comparison of the chromosome positions of DEGs with those of significant SNPs detected in the TC population identified potential candidate genes that might be related to hybrid performance. Combining RNA-seq analysis and MTA results identified candidate genes for hybrid performance, providing information that could be useful for maize hybrid breeding.

Introduction
With the increase in human population, it is expected that the production of staple food crops should be doubled to feed the growing population by 2050 [1]. Hybrid crop breeding can increase crop yield and meet human demand for food. As a successful example of hybrid breeding, hybrid maize has played an important role in maize yield increase in the last century [2]. To further increase the yield of hybrid maize by using molecular design breeding, it is necessary to dissect the genetic and molecular mechanism of hybrid performance in hybrid populations.
In maize, hybrid breeding generally requires selection of breeding materials according to line per se (LPS) and testcross (TC) performance [3]. However, LPS performance can only indirectly assess TC performance [4], and the accuracy of selection based on LPS performance depends on the relationship between the LPS and TC population. The genomic differences between an LPS population and its TC counterpart are the replacement of one half

Phenotypic Data of the Three Populations
The PH6WC TC population had the highest HKW and YPP across all the environments ( Figure 1, Table 1), followed by the Chang7-2 TC population, suggesting that PH6WC generally had better combining ability with the LPS population than Chang7-2. Both traits of the LPS population had a higher coefficient of variance and broad-sense heritability, indicating that the genetic variation was greater and more stable than that of the two TC populations. Furthermore, HKW had higher heritability than YPP, indicating that HKW was more stable than YPP across environments (Table 1). The correlation of either HKW or YPP across environments was generally significant for each population ( Figure S1), suggesting that the genetic basis plays a major role in determining the two traits across environments. SD-standard deviation, N-population size, CV-coefficient of variance, H 2 -broad-sense heritability.

Genotypic Data Analysis and Genetic Dissection of Yield Traits of the Three Populations
In total, 15,386 SNPs were obtained after genotype processing, and these SNPs were distributed evenly across the maize physical map. The number of SNPs ranged from 1056 on chromosome 10 to 2477 on chromosome 1, and the SNP density ranged from 6.62 SNP/Mb on chromosome 2 to 8.80 on chromosome 9, with a mean density of 7.5 SNP/Mb (Table S1). The SNP density was sufficiently high for MTA analysis [17]. MTA found that the numbers of significant SNPs for HKW were 24, 30 and 121 in the LPS population, Chang7-2 and PH6WC TC populations, respectively (Table S2)

Genotypic Data Analysis and Genetic Dissection of Yield Traits of the Three Populations
In total, 15,386 SNPs were obtained after genotype processing, and these SNPs were distributed evenly across the maize physical map. The number of SNPs ranged from 1056 on chromosome 10 to 2477 on chromosome 1, and the SNP density ranged from 6.62 SNP/Mb on chromosome 2 to 8.80 on chromosome 9, with a mean density of 7.5 SNP/Mb (Table S1). The SNP density was sufficiently high for MTA analysis [17]. MTA found that the numbers of significant SNPs for HKW were 24, 30 and 121 in the LPS population, Chang7-2 and PH6WC TC populations, respectively (Table S2)  the density of significant SNPs within a 1 Mb interval; red asterisk indicates common QTL between LPS and Chang7-2 TC populations; blue asterisk indicates common QTL between LPS and PH6WC TC populations; green asterisk indicates common QTL between Chang7-2 TC and PH6WC TC populations; dark asterisk indicates common QTL among LPS, Chang7-2 TC and PH6WC TC populations. (b-d) PVE of the large-effect SNPs controlling HKW of the LPS, Chang7-2 and PH6WC TC populations, respectively. (e-g) PVE of the large-effect SNPs controlling YPP of the LPS, Chang7-2 and PH6WC TC populations, respectively. In (b-g), only SNPs with PVE larger than 1% are shown; PVE_LPS, PVE_Chang7-2 and PVE_PH6WC indicate phenotypic variance explained by the significant SNPs in the LPS, Chang7-2 and PH6WC TC populations, respectively.

Genetic Features of the Significant SNPs
We calculated the cumulative effects of favorable genotypes of the significant SNPs. The correlation coefficients between the number of favorable genotypes and HKW were 0.28, 0.

Genetic Features of the Significant SNPs
We calculated the cumulative effects of favorable genotypes of the significant SNPs. The correlation coefficients between the number of favorable genotypes and HKW were 0.28, 0.49 and 0.36 for the LPS, Chang7-2 and PH6WC TC populations, respectively ( Figure 3a). The correlation coefficients between the number of favorable genotypes and YPP were 0.30, 0.22 and 0.38 for the LPS, Chang7-2 and PH6WC TC populations, respectively ( Figure 3b). The strong positive correlation indicates that the LPS performance and TC performance increased with the accumulation of favorable genotypes, which further proved the reliability of the MTA results.
To test the effect of these significant SNPs in predicting LPS and TC performance, we performed genomic prediction (GP) and marker-assisted selection (MAS) analysis. The analysis revealed that the prediction accuracies (PAs) of GP models were larger than those of the MAS.Sig model (Figure 4a-c), suggesting that some genetic factors were not identified due to the problem of false negatives for both traits in each population. We also found that the PAs of the MAS.Sig model were larger than those of MAS.Random model (Figure 4d-f) for each trait in each population, further proving that the significant SNPs were in linkage disequilibrium with the genes controlling the tested traits. To test the effect of these significant SNPs in predicting LPS and TC performance, we performed genomic prediction (GP) and marker-assisted selection (MAS) analysis. The analysis revealed that the prediction accuracies (PAs) of GP models were larger than those of the MAS.Sig model (Figure 4a-c), suggesting that some genetic factors were not identified due to the problem of false negatives for both traits in each population. We also found that the PAs of the MAS.Sig model were larger than those of MAS.Random model ( Figure  4d-f) for each trait in each population, further proving that the significant SNPs were in linkage disequilibrium with the genes controlling the tested traits.

Identification of Common QTLs between LPS and TC Populations
Given the significant correlation between the tested traits of the LPS population and the TC population (Figure 5a), we considered that there should be common QTLs controlling LPS and TC performance. To prove this hypothesis, we examined whether the signif-

Identification of Common QTLs between LPS and TC Populations
Given the significant correlation between the tested traits of the LPS population and the TC population (Figure 5a), we considered that there should be common QTLs controlling LPS and TC performance. To prove this hypothesis, we examined whether the significant SNPs of the LPS population take effect in the TC populations. The results showed that the significant SNPs for HKW of the LPS population explained 25.82% and 12.24% of phenotypic variance in the Chang7-2 and PH6WC TC populations, respectively ( Figure 5b). Furthermore, the significant SNPs controlling YPP of the LPS population explained 16.58% and 16.51% of phenotypic variance in the Chang7-2 and PH6WC TC populations, respectively (Figure 5b). The analysis proved that the significant SNPs controlling LPS performance also controlled TC performance.
the LPS (d), Chang7-2 TC (e) and PH6WC TC (f) populations. GP indicates GP models usi SNPs; MAS.Sig indicates MAS models using the significant SNPs; MAS.Random indicates models using the same number of randomly selected SNPs. P values indicate significant lev Wilcox.test.

Identification of Common QTLs between LPS and TC Populations
Given the significant correlation between the tested traits of the LPS population the TC population (Figure 5a), we considered that there should be common QTLs con ling LPS and TC performance. To prove this hypothesis, we examined whether the si icant SNPs of the LPS population take effect in the TC populations. The results sho that the significant SNPs for HKW of the LPS population explained 25.82% and 12.24 phenotypic variance in the Chang7-2 and PH6WC TC populations, respectively (F 5b). Furthermore, the significant SNPs controlling YPP of the LPS population expla 16.58% and 16.51% of phenotypic variance in the Chang7-2 and PH6WC TC populat respectively (Figure 5b). The analysis proved that the significant SNPs controlling performance also controlled TC performance. We further compared the locations of significant SNPs detected for each trait. Two common QTLs on chromosomes 4 and 6 were associated with HKW of the LPS, Chang7-2 and PH6WC TC populations. In addition, one QTL at the end of chromosome 6 was commonly detected in LPS and Chang7-2 TC populations; three common QTLs on chromosomes 2, 3 and 9 were identified in LPS and PH6WC TC populations; and two common QTLs on chromosomes 1 and 7 were identified between Chang7-2 TC and PH6WC TC populations ( Figure 2a, Table S2). For YPP, one common QTL on chromosome 1 was detected in LPS and Chang7-2 TC populations, one common QTL on chromosome 2 was detected in LPS and PH6WC TC populations, and two common QTLs on chromosomes 1 and 2 were detected in Chang7-2 and PH6WC TC populations (Figure 2a). The analysis revealed that there were common QTLs between each pair of the three populations, reflecting their strong phenotypic correlations (Figure 5a).

RNA-seq Analysis Identified the Candidate Genes in the Surrounding Region of the Significant SNPs
To find the candidate genes associated with hybrid performance of Chang7-2 TC lines, we found common DEGs between ZD958 and each of its parents (Zheng58 and Chang7-2) and compared the locations of these DEGs with those of significant SNPs. RNA-seq analysis identified 4593 common DEGs (Figure 6a). According to the candidate genes found in our previous article, the orthologs of 57 and 102 DEGs (Table S3) were related to the control of kernel weight and yield, respectively [18]. Among the 57 DEGs, the locations of nine DEGs were close to four HKW QTLs detected in the Chang7-2 TC population (Figure 6a, Table S3), including GRMZM2G159456, GRMZM2G399072, GRMZM2G445634, GRMZM2G420357, GRMZM2G034876, GRMZM2G092749, GRMZM2G059939, GRMZM2G328988 and GR-MZM2G034647. Meanwhile, five of the 102 DEGs were found in the surrounding regions of one YPP QTL detected in the Chang7-2 TC population, including GRMZM2G095968, GR-MZM2G159456, GRMZM2G399072, GRMZM2G445634 and GRMZM2G420357 (Figure 6a, Table S3). the locations of nine DEGs were close to four HKW QTLs detected in the Chang7-2 TC population (Figure 6a, Table S3), including GRMZM2G159456, GRMZM2G399072, GRMZM2G445634, GRMZM2G420357, GRMZM2G034876, GRMZM2G092749, GRMZM2G059939, GRMZM2G328988 and GRMZM2G034647. Meanwhile, five of the 102 DEGs were found in the surrounding regions of one YPP QTL detected in the Chang7-2 TC population, including GRMZM2G095968, GRMZM2G159456, GRMZM2G399072, GRMZM2G445634 and GRMZM2G420357 (Figure 6a, Table S3). Figure 6. Co-localization of significant SNPs and DEGs identified candidate genes for hybrid performance. (a) Co-localization of significant SNPs detected in the Chang7-2 TC population and the DEGs between ZD958 and its parental lines (Zheng58 and Chang7-2). The five circles, from inner to outer, indicate the location of common DEGs identified between ZD958 and each of its parental lines, the MTA results for YPP in the Chang7-2 TC population (red points indicate significant SNPs), the MTA results for HKW in the Chang7-2 TC population (blue points indicate significant SNPs), SNP density heatmap and the candidate genes in the surrounding regions of significant SNPs, and maize chromosomes. (b) Co-localization of significant SNPs detected in the PH6WC TC population and the DEGs identified between XY335 and its parental lines (PH6WC and PH4CV). The five circles, from inner to outer, indicate the location of common DEGs between XY958 and each of its parental lines, the MTA results for YPP in the PH6WC TC population (red points indicate significant SNPs), the MTA results for HK in the PH6WC TC population (blue points indicate significant SNPs), SNP density heatmap and the candidate genes in the surrounding regions of significant SNPs, and maize chromosomes.
To find the candidate genes in the surrounding regions of the significant SNPs detected in the PH6WC TC population, we first found 1801 DEGs that were commonly detected between XY335 and each of its parents. Among the 1801 DEGs, 12 and 20 DEGs (Table S4) were related to the control of yield traits [18]. Among the 12 DEGs, two candidate genes (GRMZM2G007288 and GRMZM5G875502) were found in the surrounding regions of two HKW QTLs detected in the PH6WC TC population. Meanwhile, three DEGs (GR-MZM2G050305, GRMZM2G034876 and GRMZM2G463904) were found in the surrounding regions of two YPP QTLs detected in the PH6WC TC population (Figure 6b, Table S4). The genes mentioned in this section could be considered candidate genes for the QTLs controlling HKW and YPP.

Discussion
Maize has rich genetic diversity and rapid linkage disequilibrium, and MTA analysis of various traits has been performed in maize [19][20][21][22]. Many methods have been developed to increase the calculation speed and statistical power of MTA, such as the general linear model (GLM), mixed linear model (MLM), etc., but most of them only work for inbred line populations [23]. At present, there are only a few published MTA methods suitable for performing MTA analysis in hybrid populations, including EMMAX software, which uses different encoding schemes to discriminate additive, dominant and recessive effects [10]; the linear mixed model, which fits additive and dominant effects [24]; and the PEPIS pipeline, which contains all genetic effects [12]. However, various genetic effects (additive, dominant and recessive) are dissected in the former two methods, which complicates the results of MTA analysis. Moreover, the epistasis effect is not considered in the former two methods. The PEPIS pipeline comprehensively dissects the main effect and calculates the LRT values of each SNP in a user-friendly manner. Therefore, PEPIS was used for MTA analysis in this study.
Identification of QTLs and candidate genes controlling agronomic traits is the basis for developing functional markers and molecular design breeding in maize [25,26]. Although many QTLs have been identified using family-based QTL mapping or association mapping [18,27], the QTLs or genes identified using family-based or homozygous lines are different from those detected using hybrid populations [26,28], and most of these QTLs were not proven functional in hybrid lines. In this study, we not only detected common QTLs in the LPS and TC populations but also detected specific QTLs in TC populations. The results showed that the genetic basis of LPS and TC populations is not completely different. The two common QTLs for HKW on chromosomes 4 and 6 with effects in all the three populations required further investigation. Although no common QTLs were detected in the three populations for YPP, there are still common QTLs shared between at least two populations, indicating that YPP might have a complex genetic basis [28]. These common QTLs indicated that manipulating QTLs in the LPS population could increase the yield traits of the hybrid population. Additionally, because the four parental lines of the tested populations were parents of the two most popular hybrid varieties in China, the detected QTLs for yield-related traits could explain why the two hybrid varieties are high-yield and popular in China.
Compared with the results of previous studies, it was found that the significant SNPs Chr3_104753320 (on chromosome 3) and Chr4_9699802 (on chromosome 4) associated with HKW were detected in the LPS and TC populations, respectively, which coincides with the results of a previous study [29]. A significant locus, Chr2_130338518, associated with YPP was detected in the LPS population, and this locus was also detected in an RIL population [29]. Because research on genetic mapping is scarce in maize hybrid populations, most QTL detected in this study were specific. This study has some limitations; although some significant epistatic QTLs (additive × additive, additive × dominant and dominant × dominant) were associated with hybrid performance [29][30][31], we only considered the main effect and did not identify their modes of inheritance [12]. However, this drawback did not influence the identification of common QTLs in LPS and TC populations or the process of finding candidate genes by colocalization.
GP has been proven as a reliable method for predicting both LPS and TC performance [32,33]. Because GP relies on the genetic basis of the population [17,34], the distance between the molecular markers and the QTLs of the target traits could influence the PAs of GP models. In this study, we used GP models to prove that the significant SNPs are reliable because the PAs of GP models fitting significant SNPs were larger than those fitting random SNPs [32]. However, the PAs of GP models fitting significant SNPs were lower than those of GP models fitting genome-wide SNPs, indicating that some QTLs were not detected in each population, which might be related to the high threshold level used in MTA analysis. Furthermore, the realness of the detected QTLs was also supported by the results that both HKW and YPP increased with increased favorable genotypes.
RNA sequencing analysis has been used to dissect the genetic basis of crop traits in combination with genetic methods such as association mapping and linkage-based QTL mapping [15,16]. In this study, in order to find candidate genes underlying hybrid performance, we identified the DEGs between F1 and its parental lines and compared the locations of these DEGs with those of the significant SNPs. The DEGs between F1 and the parental lines might be candidate genes, especially the 14 candidate genes for yield traits (Figure 6b, Tables S3 and S4). Among the 14 candidate genes for kernel weight, GRMZM2G159456, GRMZM2G399072, GRMZM2G445634, GRMZM2G034876, GRMZM2G007288 and GRMZM5G875502 are orthologous genes of rice BU1 [35], SNB [36], TIFY [37], GL2 [38], GW2 [39] and OsiEZ1 [40] genes, respectively, which have been reported to regulate rice seed weight; GRMZM2G420357, GRMZM2G092749, GRMZM2G059939 and GRMZM2G034647 are orthologs of Arabidopsis IKU1 [41], FERONIA [42], DPA4 [43] and CYCB1:4 [44] genes, respectively, and these genes are associated with Arabidopsis seed weight. Furthermore, GRMZM2G328988, GRMZM2G463904, GRMZM2G095968 and GRMZM2G050305 are orthologous to UPL3 in oilseed rape [45], to RLK7 in maize [46], to IbEXP1 in sweet potato [47] and to GmMYB14 in soybean [48], respectively. These genes were also reported to be associated with seed weight or overall yield. Therefore, the 14 DEGs for kernel weight might be candidate genes because they were close to the positions of the significant SNPs identified in the TC populations.

Population Construction, Phenotype Evaluation and Phenotypic Data Analysis
Four elite inbred lines were used in this study, including Zheng58, Chang7-2, PH6WC and PH4CV. Zheng58 and Chang7-2 were the female and male parents of ZD958, respectively. PH6WC and PH4CV were the female and male parents of XY335, respectively. ZD958 and XY335 are popular hybrid varieties in China [49,50]. The 481 BC1F3:4 families were introduced in detail in [32]. Briefly, PH4CV was used as the recurrent parent, and Zheng58 was used as the donor parent to develop a BC1F3 population, which was self-pollinated to develop the BC1F3:4 families. The BC1F3:4 families were defined as an LPS population. The 481 BC1F3 plants testcrossed to Chang7-2 and PH6WC in the winter of 2015 in Sanya (Hainan province), producing Chang7-2 and PH6WC TC populations, respectively. These materials were frequently used in maize genetic improvement research (https://maizedata.cn/, accessed on 3 October 2021) [33,51,52]. The materials and populations were provided by the molecular genetic improvement group of the Institute of Crop Sciences, Chinese Academy of Agricultural Sciences.
The LPS and TC populations were sown in Shunyi (Beijing municipality) and Changji (Xinjiang Uygur Autonomous region) in the summer of 2016 and 2017. The two TC populations were also sown in Xinxiang (Henan province) in the summer of 2017. The five environments were identified as 16BJ, 17BJ, 16XJ, 17XJ and 17HN, where BJ, XJ and HN indicate Beijing, Xinjiang and Henan locations, respectively, and 16 and 17 indicate the years 2016 and 2017, respectively. The field experimental design was an incomplete block design, as explained in detail in in our previous publication [33]. The row length and row space were 5 m and 60 cm, respectively, and the planting density was 4444 plants per mu, where mu is a traditional Chinese unit for measuring field size. At the harvest stage, the yield of each plot was measured and adjusted to 14% water content. YPP was calculated by dividing the plot yield into the number of plants in the plot. The hundred kernel weight (HKW) of each plot was measured manually. For HKW, the Chang7-2 and PH6WC TC populations were evaluated in five environments (16BJ, 17BJ, 16XJ, 17XJ and 17HN), whereas the LPS population was evaluated in four environments (16BJ, 17BJ, 16XJ and 17XJ). For YPP, the Chang7-2 and PH6WC TC populations were evaluated in four environments (16BJ, 17BJ, 17XJ and 17HN), whereas the LPS population was evaluated in three environments (16BJ, 17BJ and 17XJ). All experimental research on plants, including collection of plant materials, complied with institutional, national or international guidelines. Field studies were conducted in accordance with local legislation.
The model for calculating BLUEs is as follows [53,54]: where y ikmb is the phenotypic data of the i th genotype in the b th block nested in the m th replication that is nested in the k th environment, µ is the overall mean, g i is the genotype effect, τ k is the environmental effect, gτ ik is the G × E effect, δ (k)m is the replication effect nested in each environment, β (m)b is the block effect nested in the replication effect, and ε ikmb is the residual error. When calculating the BLUEs, the other variables (except the genotype) are treated as random effects and assume to follow normal distributions. All factors are random effects when calculating broad-sense heritability; the formula is as follows [55]: where σ 2 g is the genotype variance; σ 2 gτ is the variance of G × E; σ 2 ε is the error variance; and n and r are the number of environments and replicates, respectively. The models for calculating BLUEs and broad-sense heritability are fitted by using the lme4 package [56]. The phenotypic data of the LPS and TC populations are available in Table S5.

Genotype Processing
Leaf samples of the 481 BC1F3 plants and two testers (Chang7-2 and PH6WC) were sampled, and DNA of these samples was extracted following the cetyltrimethyl ammonium bromide method [57]. DNA was sent to Beijing Capital Bio for genotyping using DNA chips containing 55,000 SNPs [58]. The genotypic data were filtered by following the steps: (1) SNPs with a calling rate lower than 97% were removed; (2) SNPs with no physical position information were removed; (3) SNPs with a missing rate greater than 5% were removed; (4) SNPs with minor allele frequencies lower than 0.05 were removed; and (5) the missing genotypes were input using the codeGen function of the R package "synbreed", the method "beagle" was used and the other settings were default [59,60]. The minor and major alleles were coded as 2 and 0, respectively, and the heterozygous genotypes were coded as 1 [60]. The genotypes of the testcross population were deduced from the testers (Chang7-2 and PH6WC) and 481 BC1F3 parents. Because some loci of the BC1F3 population were heterozygous ( Figure S2), we deduced the genotype of the testcross population as described by Cui et al. [61]. Assuming the genotype code of a loci is defined as aa = 0, Aa = 1, AA = 2, if the genotype of the tester is aa, the genotypes of the testcross progeny would be 0, 0.5 and 1; if the genotype of the tester is AA, the genotypes of the testcross progenies would be 1, 1.5 and 2. The genotypic data of the LPS, Chang7-2 TC and PH6WC TC populations are available in Table S6.

Marker-Trait Association Analysis and Calculation of PVE
In this study, we used PEPIS software to perform MTA analysis. PEPIS software is one of the few public user-friendly tools for performing genetic mapping of hybrid populations [12]. The PEPIS software package is based on a linear mixed model [62], and the statistical method of PEPIS is as follows: First, the genotype of individual j in marker k is encoded into two numerical variables: where Z jk and W jk are indicators of additive and dominant effects, respectively. A (the first homozygous genotype), H (heterozygous genotype) and B (the second homozygous genotype) indicate genotypes of each marker. Then, the following statistical model is used: where y is the n × 1 vector of the phenotypic data (BLUE); Xβ is a non-genetic effect; and a k and d k are the additive and dominance effects, respectively. For markers k and k , (aa) kk , (ad) kk , (da) kk , (dd) kk are additive × additive, additive × dominant, dominant × additive and dominant × dominant epistatic effects, respectively [63]. For each population, we first constructed the additive genotype matrix and the dominance genotype matrix. Then, we input the two matrices and BLUE data into PEPIS software to run marker-trait association analysis (http://bioinfo.noble.org/PolyGenic_QTL, accessed on 7 June 2019). The LRT threshold for declaring significance was −log10(0.05/the number of markers) according to the PEPIS pipeline [12]. The LRT threshold should be 5.49, given that the number of markers is 15,386. The SNPs with LRT values over 5.49 were identified as significant SNPs.
To calculate the phenotypic variance explained by SNPs (PVE), the significant SNPs were fitted in a multiple linear model [64], from which SSreg and SStol for each SNP were computed. SSreg is the sum of square of each SNP, whereas SStol is the sum of square of the linear model. The PVE of each SNP was calculated by dividing SSreg into SStol.

GP and MAS Analysis
We used the ridge regression best linear unbiased prediction (rrBLUP) model to run GP analysis. The rrBLUP model is [65]: where y is the BLUEs, β is a vector of the fixed effects including only the overall mean, u is the vector of random effects including only additive effect, ε is the residual error, X and Z are the design matrices. GP was implemented by running five-fold cross validation for 200 repeats. The effects of genome-wide markers were estimated, and the predicted phenotypic values were calculated by inputting the effects of genome-wide markers into the rrBLUP model. PA was calculated as the Pearson correlation coefficient between the observed and predicted phenotype. The R package rrBLUP was used to implement the GP model [66], and the code is available in File S1.
To calculate the PA of MAS model fitting the significant SNPs (defined as MAS.Sig model), a multiple regression model was fitted using the lm function in R. The phenotype was predicted using the predict function [53]. The PA was calculated by running five-fold cross validation for 200 repeats. In order to prove the effect of the MAS model fitting significant SNPs, we also calculated the PA of the MAS model fitting the same number of randomly selected SNPs. The model was defined as the MAS.Random model. We used Wilcox.test to compare the differences among the PAs of the three prediction models (GP, MAS.Sig and MAS.Random models).

Identification of Common QTLs among LPS and two TC Populations
Previously, all QTLs within a 20 cM interval were considered a single QTL [67]. According to a previous report, the average recombination rate was 1 cM/Mb [68], which is approximately 1 Mb. Therefore, we defined SNPs within 20 Mb as in linkage with one QTL.

RNA-seq Analysis and Identification of Differentially Expressed Genes around Significant SNPs
The six materials, including two hybrids (ZD958 and XY335), and their parents (Zheng58, Chang7-2, PH6WC and PH4CV) were sown in July 2018 in Haidian, Beijing. Decapitated shoot tips at the V7 stage were used for RNA extraction. For each material, three biological replicates were used, with each replicate containing three samples. We extracted RNA using an RNeasy plant mini kit (Qiagen, Germany) and checked RNA purity using a kaiaoK5500 spectrophotometer (Kaiao, Beijing, China). Then, we assessed RNA integrity and concentration using an RNA Nano 6000 assay kit of a Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA). After purifying mRNA from total RNA using poly-T oligo-attached magnetic beads, we generated the sequencing libraries using NEBNext ® Ultra™ Directional RNA Library Prep Kit for Illumina (NEB, Ispawich, MA, USA). The libraries were sequenced using the Illumina Novaseq system with a read length of 150 bp (pair end) at Annoroad Gene Technology (Beijing, China).
The RNA sequencing data were analyzed according to the procedure used in our laboratory [69]. Briefly, the raw data were filtered to remove low-quality reads, adaptor-polluted reads and reads with more than 5% N bases. The filtered clean reads were mapped to the B73 RefGen_V3 genome (www.maizegdb.org, accessed on 5 August 2019) using Hisat2 with default settings. The expression level of each sample was estimated using FPKM (fragments per kilobase of transcripts per million fragments mapped), which is calculated by normalizing raw reads. The threshold for identifying the DEGs was false discovery rate (FDR, p value < 0.05), which was computed by using Cufflinks. The code is available in File S1.
Some studies have shown that differentially expressed genes between hybrids and their parental lines contribute to hybrid performance [70][71][72][73]. Therefore, we first identified the DEGs between ZD958 and its parental lines (Zheng58 and Chang7-2) and between XY335 and its parental lines (PH6WC and PH4CV). For the Chang7-2 TC population, the physical positions of common DEGs between ZD958 and each of its parents were compared to those of the significant SNPs identified in the Chang7-2 TC population. Those DEGs located within 20 Mb of significant SNPs were identified as potential candidate genes. In the same way as stated above, common DEGs between XY335 and each of its parental lines were identified, and the physical positions of these DEGs were compared to those of the significant SNPs identified in the PH6WC TC population. In the same way, those DEGs located within 20 Mb of significant SNPs were identified as potential candidate genes. To prioritize the candidate genes for HKW, the DEGs in the surrounding regions of significant SNPs associated with HKW were compared to the candidate genes for kernel weight [18]. In the same way, in order to prioritize the candidate genes for YPP, the DEGs in the surrounding regions of significant SNPs associated with YPP were compared to the candidate genes for yield and yield-related traits [18].  Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: The phenotype and genotype data used in this study are available as supplementary files. The RNA-seq data were deposited in NCBI, (BioProject ID PRJNA766146). The code is available in File S1.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.