Comparative Transcriptome Combined with Proteome Analyses Revealed Key Factors Involved in Alfalfa (Medicago sativa) Response to Waterlogging Stress

Alfalfa (Medicago sativa) is the most widely grown and most important forage crop in the world. However, alfalfa is susceptible to waterlogging stress, which is the major constraint for its cultivation area and crop production. So far, the molecular mechanism of alfalfa response to the waterlogging is largely unknown. Here, comparative transcriptome combined with proteomic analyses of two cultivars (M12, tolerant; M25, sensitive) of alfalfa showing contrasting tolerance to waterlogging were performed to understand the mechanism of alfalfa in response to waterlogging stress. Totally, 748 (581 up- and 167 down-regulated) genes were differentially expressed in leaves of waterlogging-stressed alfalfa compared with the control (M12_W vs. M12_CK), whereas 1193 (740 up- and 453 down-regulated) differentially abundant transcripts (DATs) were detected in the leaves of waterlogging-stressed plants in comparison with the control plants (M25_W vs. M25_CK). Furthermore, a total of 187 (122 up- and 65 down-regulated) and 190 (105 up- and 85 down-regulated) differentially abundant proteins (DAPs) were identified via isobaric tags for relative and absolute quantification (iTRAQ) method in M12_W vs. M12_CK and M25_W vs. M25_CK comparison, respectively. Compared dataset analysis of proteomics and transcriptomics revealed that 27 and eight genes displayed jointly up-regulated or down-regulated expression profiles at both mRNA and protein levels in M12_W vs. M12_CK comparison, whereas 30 and 27 genes were found to be co-up-regulated or co-down-regulated in M25_W vs. M25_CK comparison, respectively. The strongly enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways for co-up-regulated genes at mRNA and protein levels in M12_W vs. M12_CK comparison were ‘Amino sugar and nucleotide sugar metabolism’, ‘Arginine and proline metabolism’ and ‘Starch and sucrose metabolism’, whereas co-up-regulated protein-related pathways including ‘Arginine and proline metabolism’ and ‘Valine, leucine and isoleucine degradation’ were largely enriched in M25_W vs. M25_CK comparison. Importantly, the identified genes related to beta-amylase, Ethylene response Factor (ERF), Calcineurin B-like (CBL) interacting protein kinases (CIPKs), Glutathione peroxidase (GPX), and Glutathione-S-transferase (GST) may play key roles in conferring alfalfa tolerance to waterlogging stress. The present study may contribute to our understanding the molecular mechanism underlying the responses of alfalfa to waterlogging stress, and also provide important clues for further study and in-depth characterization of waterlogging-resistance breeding candidate genes in alfalfa.


Introduction
Alfalfa (Medicago sativa) is the most widely grown and most important forage crop in the world [1,2]. It is of great economic, environmental, and social values, such as providing food for livestock, enhancing soil fertility, N 2 fixation, eliminating water contamination, and biofuels production [3,4]. However, the alfalfa is susceptible to waterlogging or water stress caused by excess rainfall and excessive irrigation, and this susceptibility is the major constraint for its cultivation area and crop production [2,5]. Therefore, to overcome the effects of waterlogging for alfalfa, it is essential to explore how alfalfa responds to the waterlogging.
As a serious abiotic stress, waterlogging threatens the growth and survival of plants because of oxygen deprivation [6]. Plants try to make metabolic or morphological changes, including glycolysis and formation of aerenchyma, to cope with the low oxygen conditions [7,8], and the process is always accompanied by specific alteration of transcription and translation. Genome-scale transcription studies have been performed in Arabidopsis, rice, maize, cotton, and other species in waterlogging stress [9][10][11]. For example, many transcription factors, belonging to AP2/ERF, WRKY, TGA, MYB, and bZIP families, were differentially expressed under waterlogging stress in kiwifruit [12]. Particularly, the ERF always accounted for the highest number of differentially expressed transcription factors under waterlogging conditions [13]. As reported, genes in response to waterlogging are mainly involved in ethylene synthesis, carbohydrate catabolism, lipid metabolism, glycolysis, ethanol fermentation, auxin-mediated processes, and reactive oxygen species (ROS) generation [14,15].
In addition, proteomic analysis has been widely used to detect differentially regulated proteins under waterlogging in maize, tomato, soybean, wheat [16][17][18][19][20]. Many anaerobically induced polypeptides (ANPs) which are essential for hypoxic tolerance have been identified by proteomics research in diverse plants [21]. Most of these ANPs are involved in the fermentation and glycolysis process, including sucrose synthase, glucose-6-phosphate isomerase, aldolase, and enolase [22,23]. Alam et al. discovered several novel waterlogging-responsive proteins that were not known previously as being waterlogging responsive except for the above well-known classical anaerobically induced proteins. The novel proteins were involved in several processes, i.e., signal transduction, programmed cell death, RNA processing, redox homeostasis, and metabolisms of energy [24].
Despite the molecular mechanism of the response to waterlogging have been reported in many studies, large differences in the transcriptional and translational regulation of above pathways are observed between species [25], and there have been relatively few studies of transcriptomic and proteomic changes in alfalfa until now. In this study, to understand the mechanism of the waterlogging response in alfalfa, the high throughput RNA-sequencing and iTRAQ analyses of two cultivars (M12 and M25) of alfalfa showing different tolerances to waterlogging stress were performed. Lists of candidate genes which may be related to waterlogging response in alfalfa were obtained by correlation analysis of the differentially abundant transcripts and proteins. This work will enhance understanding the mechanisms of waterlogging response in alfalfa, and may aid in stress resistance breeding of alfalfa.

Physiological Response to Waterlogging
Under well-watered control conditions, there was no significant difference between two cultivars ( Figure 1A,B). When subjected to 12 days of waterlogging stress, the M25 showed a remarkable symptom of leaf senescence and chlorosis, as indicated by the decreased chlorophyll content ( Figure 1B), whereas a slight reduction in chlorophyll content was observed in M12 after 12 days of waterlogging. Leaf Fv/Fm significantly declined for M25 after waterlogging, which reduced to 50% of the control level. For M12, however, Fv/Fm decreased only 10% when compared to the control plants ( Figure 1C). Leaf P n remarkably declined for M12 and M25 after waterlogging, which decreased by 31% and 66% when compared to the respective control plants ( Figure 1D). Int

Transcriptome Sequencing and Assembly
To investigate the genes associated with waterlogging stress response, four cDNA libraries were constructed from total RNA extracted from leaves of alfalfa (M12 and M25) with or without waterlogging treatment. The libraries were then sequenced using the Illumina HiSeq™ 2000 platform. An overview of sequence assembly after illumina sequencing was shown in Table 1. A total of 180,507 transcripts were obtained from the clean reads with a mean length of 726 bp and length ranging from 201 to 15,720 bp (Table 2). Furthermore, 112,464 unigenes were obtained with an average length of 995 bp. The N50 and N90 for unigenes was 1448 and 456 bp, respectively.

Transcriptome Sequencing and Assembly
To investigate the genes associated with waterlogging stress response, four cDNA libraries were constructed from total RNA extracted from leaves of alfalfa (M12 and M25) with or without waterlogging treatment. The libraries were then sequenced using the Illumina HiSeq™ 2000 platform. An overview of sequence assembly after illumina sequencing was shown in Table 1. A total of 180,507 transcripts were obtained from the clean reads with a mean length of 726 bp and length ranging from 201 to 15,720 bp (Table 2). Furthermore, 112,464 unigenes were obtained with an average length of 995 bp. The N50 and N90 for unigenes was 1448 and 456 bp, respectively.
There were 13,371 genes assigned to euKaryotic Ortholog Groups (KOG) classification and divided into 25 specific categories ( Figure 3). The largest number of genes belonged to 'general functional prediction only' (1796) category, while only a few genes were divided into 'Cell motility' category.  There were 13,371 genes assigned to euKaryotic Ortholog Groups (KOG) classification and divided into 25 specific categories ( Figure 3). The largest number of genes belonged to 'general functional prediction only' (1796) category, while only a few genes were divided into 'Cell motility' category.

Differentially Abundant Transcripts (DATs) Analysis
Totally, 748 (581 up-and 167 down-regulated) genes were differentially expressed in leaves of waterlogging-stressed alfalfa compared with control (M12_W vs. M12_CK) ( Figure 4A Figure 4C). Further hierarchical clustering method was adopted to observe the overall expression pattern of DATs ( Figure 4D). Most of the DATs showed remarkable differences in expression levels in waterlogging-treated conditions as compared to controls ( Figure  4D).   Figure 4C). Further hierarchical clustering method was adopted to observe the overall expression pattern of DATs ( Figure 4D). Most of the DATs showed remarkable differences in expression levels in waterlogging-treated conditions as compared to controls ( Figure 4D).

Function Annotation of DATs Using the KEGG Database
The significantly enriched pathway terms were determined based on the p-Value and corrected p-Value according to the number of DATs in each term of the KEGG pathway in comparison with the background number. The top-four KEGG pathways associated with 'Photosynthesis-antenna proteins', 'Arginine and proline metabolism', 'α-Linolenic acid metabolism', and 'Nitrogen metabolism' were significantly enriched in the waterlogging-treated alfalfa versus the controls (M12_W vs. M12_CK) (q < 0.05). By contrast, the significantly enriched pathway terms in M25_W vs. M25_CK comparison were as follows: 'Photosynthesis', 'Photosynthesis -antenna proteins',

Function Annotation of DATs Using the KEGG Database
The significantly enriched pathway terms were determined based on the p-Value and corrected p-Value according to the number of DATs in each term of the KEGG pathway in comparison with the background number. The top-four KEGG pathways associated with 'Photosynthesis-antenna proteins', 'Arginine and proline metabolism', 'α-Linolenic acid metabolism', and 'Nitrogen metabolism' were significantly enriched in the waterlogging-treated alfalfa versus the controls (M12_W vs. M12_CK) (q < 0.05). By contrast, the significantly enriched pathway terms in M25_W vs. M25_CK comparison were as follows: 'Photosynthesis', 'Photosynthesis -antenna proteins', 'Carbon fixation in photosynthetic organisms', 'Nitrogen metabolism', 'Glyoxylate and dicarboxylate metabolism', 'Valine, leucine and isoleucine degradation', 'Arginine and proline metabolism', 'Porphyrin and chlorophyll metabolism', 'Fatty acid degradation','Ubiquinone and other terpenoid-quinone biosynthesis', and 'Cysteine and methionine metabolism' (Table 3).

Proteome-Wide Analysis of Differentially Abundant Proteins by Waterlogging Treatment
To further explore the differentially abundant proteins (DAPs) regulated by waterlogging, the iTRAQ method was applied. The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium via the Proteomics Identifications (PRIDE) partner repository with the dataset identifier PXD013025. Totally, 3977 proteins were identified from alfalfa, among which 3436 proteins were quantified (

Gene Ontology (GO) Analysis of DAPs
During the M12_W vs. M12_CK comparison, 'metabolic process', 'cellular process', 'single-organism process', and 'response to stimulus' were the predominant components for 'Biological Process' category; while 'cell', 'macromolecular complex' and 'organelle' were highly represented in the 'Cellular Component' category. In the 'Molecular Function' category, proteins involved in 'catalytic activity', 'binding', 'structural molecule activity', and 'antioxidant activity' were the core components (Supplementary File 4). The similar results were also revealed in M25_W vs. M25_CK comparison (Supplementary File 5).

Gene Ontology (GO) Analysis of DAPs
During the M12_W vs.
M12_CK comparison, 'metabolic process', 'cellular process', 'single-organism process', and 'response to stimulus' were the predominant components for 'Biological Process' category; while 'cell', 'macromolecular complex' and 'organelle' were highly represented in the 'Cellular Component' category. In the 'Molecular Function' category, proteins involved in 'catalytic activity', 'binding', 'structural molecule activity', and 'antioxidant activity' were the core components (Supplementary File 4). The similar results were also revealed in M25_W vs. M25_CK comparison (Supplementary File 5).

Transcriptomes and Proteomics Crosstalk Analysis
Compared dataset of proteomics and transcriptomics in M12_W vs. M12_CK comparison, 3853 proteins or transcripts were identified both in proteomics and transcriptomics research. Among which transcript and protein expression level of 3681 genes remained unchanged; 27 and eight genes displayed jointly up-regulated or down-regulated expression profiles at both mRNA and protein levels, respectively (Table 4). On the other hand, there were 3851 proteins or transcripts identified both in proteomics and transcriptomics research in M25_W vs. M25_CK comparison. Transcript or protein expression levels of 3475, 30, and 27 genes were found to be unchanged, co-up-regulated or co-down-regulated, respectively (Table 5). To confirm the transcript expression both at proteomic and transcriptomic levels, nine co-up-regulated or co-down-regulated target genes were selected to validate the RNA-seq and iTRAQ results by qPCR, and the results showed the same expression tendency as the RNA-seq and iTRAQ results (Supplementary File 6).

Transcriptomes and Proteomics Crosstalk Analysis
Compared dataset of proteomics and transcriptomics in M12_W vs. M12_CK comparison, 3853 proteins or transcripts were identified both in proteomics and transcriptomics research. Among which transcript and protein expression level of 3681 genes remained unchanged; 27 and eight genes displayed jointly up-regulated or down-regulated expression profiles at both mRNA and protein levels, respectively (Table 4). On the other hand, there were 3851 proteins or transcripts identified both in proteomics and transcriptomics research in M25_W vs. M25_CK comparison. Transcript or protein expression levels of 3475, 30, and 27 genes were found to be unchanged, co-up-regulated or co-down-regulated, respectively (Table 5). To confirm the transcript expression both at proteomic and transcriptomic levels, nine co-up-regulated or co-down-regulated target genes were selected to validate the RNA-seq and iTRAQ results by qPCR, and the results showed the same expression tendency as the RNA-seq and iTRAQ results (Supplementary File 6).

Discussion
Waterlogging restricts O 2 diffusion and thereby inhibits aerobic respiration to plants. One of the adaptive strategies that plants survive periods of hypoxia (limited O 2 ) or anoxia (no O 2 ) is to generate energy through fermentative metabolism [26]. The fermentation, which is the main carbohydrate metabolism pathway during anaerobic condition, generates ATP and essential metabolites such as Ala, lactate, ethanol, and acetate in plants. During waterlogging stress, the energy metabolism of the plant is in a state of demand greater than supply. Moreover, the content of soluble sugar in plants is limited. Thus, the content of soluble sugar is vital for plants against waterlogging stress-induced injury. Studies have revealed that amylases play major roles in catalyzing carbohydrates conversion into soluble sugar [27]. In the present study, the expression of one gene (Cluster-1252.44338) encoding beta-amylase was induced 2.8-fold in waterlogging-resistant cultivar M12 after waterlogging treatment, while its expression was not significantly changed in waterlogging-sensitive cultivar M25. The enriched KEGG pathways consisting of co-up-regulated genes at both mRNA and protein levels were 'Amino sugar and nucleotide sugar metabolism', 'Arginine and proline metabolism', and 'Starch and sucrose metabolism' in M12_W vs. M12_CK comparison, and 'Arginine and proline metabolism' and 'Valine, leucine and isoleucine degradation' pathways in M25_W vs. M25_CK comparison. These results suggested that pathways related to carbohydrate and amino metabolism were activated when Alfalfa was exposed to waterlogging stress.
Plants adaptation to submergence stress includes two different opposing strategies, the quiescence strategy and the escape strategy. Plants temporarily stop shoot elongation to conserve energy as a quiescence strategy and resume growth when the water level is reduced. Conversely, as for escape strategy, plants keep their leaves above the water surface by stem elongation. The SUB1 locus on chromosome 9 containing a cluster of three group VII ethylene response factor (ERF) genes (SUB1A, SUB1B, and SUB1C), was identified to be involved in regulating the quiescence strategy via quantitative trait locus (QTL) mapping. The SUB1A-1 suppresses ethylene production, resulting in reduction of GA synthesis; the SUB1A-1 allele presence restricts undersurface shoot growth. Interestingly, two ERF family genes in deepwater rice named SNORKEL1 and SNORKEL2 (SK1/2), were revealed as positive regulators of internode elongation. In the present study, Ethylene-responsive transcription factor ERF110 (Cluster-1252.13837) was significantly induced by waterlogging treatment in both M12 and M25 cultivars, but the induction fold was higher in M25 (Log 2 -fold change 3.47) than that of M12 (Log 2 -fold change 2.29). Besides, ERF061 (Cluster-1252.45693) gene expression was up-regulated in M25 cultivar, but not in M12 cultivar. These results revealed that ERF genes play essential roles in alfalfa response to waterlogging. The differential expression levels of ERF110 gene in waterlogging resistant cultivar M12 and waterlogging sensitive cultivar M25, and specifically expressed ERF061 gene in M25 under waterlogging treatment condition may result in the different flood resistance ability. The roles of these ERF genes in waterlogging stress in alfalfa need to be further investigated.
Calcium signals are involved in plant responses to various stimuli, including abiotic-and biotic stresses, and regulate a wide range of physiological processes [28]. Ca 2+ has also been proposed to regulate low O 2 signaling in plants as the studies revealed that cytosolic Ca 2+ concentration was transiently increased after flooding of maize roots or anoxic or hypoxic treatment of Arabidopsis [29,30]. Calcineurin B-like (CBL) proteins and their interacting protein kinases (CIPKs) transduce plant Ca 2+ signaling through a complex network [31]. Lee et al. (2009) reported that the CIPK15-SnRK1A-MYBS1-mediated sugar-sensing pathway contributes to O 2 deficiency tolerance during rice seed germination [32]. A very recent study has revealed that natural variation in OsCBL10 promoter is involved in flooding stress response during seed germination among rice subspecies [33]. Here, the expressions of two CIPK genes (Cluster-1252.49669, Cluster-1252.47248) were obviously induced by waterlogging stress in waterlogging-resistant cultivar M12, whereas only one CIPK gene transcripts was accumulated by waterlogging treatment in waterlogging-sensitive cultivar M25 (Cluster-1252.43363). These results suggested that CIPK15-mediated calcium signals also play crucial roles in waterlogging stress in alfalfa.
Different abiotic stresses including waterlogging are known to cause the accumulation of reactive oxygen species (ROS) such as singlet oxygen, superoxide anion, hydrogen peroxide, and hydroxyl radicals which lead to membrane lipid peroxidation [34]. To minimize ROS induced injury for plant cells, plants adopt an antioxidant defense system, including non-enzymatic and enzymatic components to scavenge the excess ROS [35,36]. Among enzymatic antioxidants, Glutathione peroxidase (GPX) and Glutathione-S-transferase (GST) play important roles in protecting organisms from oxidative stress. In the present research, transcripts of five genes encoding GST and one gene encoding GPX were all obviously induced by waterlogging treatment in both M12 and M25 cultivars; however, the induction folds of these six genes were higher in M12 than those of M25 (Supplementary File 7). These results revealed that GST and GPX genes may play an important function in conferring M12 cultivar with enhanced tolerance to waterlogging.

Plant Materials and Waterlogging Treatment
Two alfalfa cultivars with contrasting waterlogging tolerance obtained from Beijing RYTWAY, M12 (tolerant) and M25 (sensitive), previously determined in our preliminary test, were used in this study. Equal amount of seeds for two cultivars were planted in pots with nursery substrate: river sand mix (2:1 v/v) in a greenhouse at Hunan Agricultural University under natural lighting in March 17th, 2017, with 18 pots per species. Plants were thinned to 15 seedlings in each pot two weeks after sowing. Four weeks' old plants were moved to a growth chamber (Percival, Boone, Iowa, USA) with the following growth conditions: 25 • C/23 • C temperature (day/night), 14/10 h photoperiod, 500 µmol·m −2 ·s −1 light intensity, and 60-85% relative humidity. After one week of acclimation, plants were randomly assigned to "control" and "waterlogging" groups. For treatment, plants were subjected to waterlogging by immersing the plastic pots into water-filled plastic tubs by maintaining 1 cm water layer above soil surface, whereas the control pots were watered regularly to keep as 100% field soil water capacity. After 12 days' treatment, samples were frozen by liquid nitrogen and were stored in −80 • C until use.

Leaf Chlorophyll Content, Maximum Quantum Yield of Photosystem II Efficiency (Fv/Fm), and Net Photosynthetic Rate (Pn) Determination
Leaf chlorophyll content was measured on the forth leaves from top by using a hand-held chlorophyll meter (SPAD-502, Spectrum Technologies Inc., Plainfield, IL, USA). Leaf maximum quantum yield of photosystem II efficiency (Fv/Fm) was evaluated by using a chlorophyll fluorometer (OS1-FL, Opti-Sciences, Hudson, NH, USA). Plants were adapted in darkness for 30 min and the then the measurements were made on intact leaves with the fluorometer.

RNA Preparation, Sequencing, and Data Analysis
Leaves total RNA was extracted using Trizol reagent (Invitrogen, Carlsbad, CA, USA) and purified using the RNeasy Plant Mini kit (Qiagen, Hilden, Germany) according to the Handbook. The quality and integrity of RNA was checked by Agilent Bioanalyzer 2100 system (Agilent Technologies, Palo Alto, CA, USA) and agarose gel electrophoresis.
A total amount of 1.5 µg RNA per sample was used as input material for the RNA sample preparations. Sequencing libraries were generated using NEBNext ® Ultra™ RNA Library Prep Kit for Illumina ® (NEB, Ipswich, MA, UK) following manufacturer's recommendations and index codes were added to attribute sequences to each sample. Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. Fragmentation was carried out using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer (5×). First strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase (RNase H). Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylation of 3 ends of DNA fragments, NEBNext Adaptor with hairpin loop structure were ligated to prepare for hybridization. In order to select cDNA fragments of preferentially 150-200 bp in length, the library fragments were purified with AMPure XP system (Beckman Coulter, Beverly, MA, USA). Then 3 mm 3 USER Enzyme (NEB) was used with size-selected, adaptor-ligated cDNA at 37 • C for 15 min followed by 5 min at 95 • C before PCR. Then PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers and Index (X) Primer. At last, PCR products were purified (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumia, Santiago, CA, USA) according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq platform and paired-end reads were generated. Raw data (raw reads) of fastq format were firstly processed through in-house perl scripts. In this step, clean data (clean reads) were obtained by removing reads containing adapter, reads containing ploy-N and low quality reads from raw data. At the same time, Q20, Q30, GC-content and sequence duplication level of the clean data were calculated. All the downstream analyses were based on clean data with high quality. Transcriptome assembly was accomplished using Trinity [37] with min_kmer_cov set to 2 by default and all other parameters set default.
Gene expression levels were estimated by RSEM [38] for each sample. Prior to differential gene expression analysis, for each sequenced library, the read counts were adjusted by edge R program package through one scaling normalized factor. Differential expression analysis of two samples was performed using the DEGseq [39] R package. p-value was adjusted using q value [40]. The q-value was set as the threshold for significantly differential expression.
Gene Ontology (GO) enrichment analysis of the differentially abundant transcripts (DATs) was implemented by the GOseq R packages based Wallenius non-central hyper-geometric distribution [41], which can adjust for gene length bias in DATs. KEGG [42] is a database resource for understanding high-level functions and utilities of the biological system, such as the cell, the organism and the ecosystem, from molecular-level information, especially large-scale molecular datasets generated by genome sequencing and other high-throughput experimental technologies (http://www.genome.jp/ kegg/). We used KOBAS [43] software to test the statistical enrichment of differential expression genes in KEGG pathways.

Protein Extraction and Trypsin Digestion
The sample was ground using a mortar and a pestle in liquid nitrogen into cell powder and then transferred to a 5 cm 3 centrifuge tube. After that, four volumes of lysis buffer (8 M urea, 1% Triton-100, 10 mM dithiothreitol, and 1% Protease Inhibitor Cocktail) was added to the cell powder, followed by sonication three times on ice using a high intensity ultrasonic processor (Ningbo Scientz Biotechnology Co., Ningbo, China). The remaining debris was removed by centrifugation at 20,000 g at 4 • C for 10 min. Finally, the protein was precipitated with cold 20% TCA for 2 h at −20 • C. After centrifugation at 12,000× g 4 • C for 10 min, the supernatant was discarded. The remaining precipitate was washed with cold acetone for three times. The protein was redissolved in 8 M urea and the protein concentration was determined with BCA kit according to the manufacturer's instructions.
For digestion, the protein solution was reduced with 5 mM dithiothreitol for 30 min at 56 • C and alkylated with 11 mM iodoacetamide for 15 min at room temperature in darkness. The protein sample was then diluted by adding 100 mM TEAB to urea concentration less than 2 M. Finally, trypsin was added at 1:50 trypsin-to-protein mass ratio for the first digestion overnight and 1:100 trypsin-to-protein mass ratio for a second 4 h-digestion.

iTRAQ Labeling
After trypsin digestion, peptide was desalted by Strata X C18 SPE column (Phenomenex) and vacuum-dried. Peptide was reconstituted in 0.5 M TEAB and processed according to the manufacturer's protocol for iTRAQ kit. Briefly, one unit of iTRAQ reagent were thawed and reconstituted in acetonitrile. The peptide mixtures were then incubated for 2 h at room temperature and pooled, desalted and dried by vacuum centrifugation. Then, pooled proteins were then labeled with iTRAQ reagents as follows M12CK: 113 and 114 (control), M12W: 115 and 116, M25CK: 117 and 118, and M25W: 119 and 121, and incubated at room temperature for 2 h. All labeled peptides were then combined before being dried in a vacuum concentrator (Eppendorf Concentrator 5301, Hamburg, Germany).

HPLC Fractionation and LC-MS/MS Analysis
The tryptic peptides were fractionated into fractions by high pH reverse-phase HPLC using Agilent 300Extend C18 column (5 µm particles, 4.6 mm ID, 250 mm length). Briefly, peptides were first separated with a gradient of 8% to 32% acetonitrile (pH 9.0) over 60 min into 60 fractions. Then, the peptides were combined into 18 fractions and dried by vacuum centrifuging. The tryptic peptides were dissolved in 0.1% formic acid (solvent A), directly loaded onto a home-made reversed-phase analytical column (15-cm length, 75 µm i.d.). The gradient was comprised of an increase from 6% to 23% solvent B (0.1% formic acid in 98% acetonitrile) over 26 min, 23% to 35% in 8 min and climbing to 80% in 3 min then holding at 80% for the last 3 min, all at a constant flow rate of 400 µm 3 /min on an EASY-nLC 1000 UPLC system (Thermo Fisher Scientific, Bremen, Germany).
The peptides were subjected to NSI source followed by tandem mass spectrometry (MS/MS) in Q Exactive™ Plus (Thermo Fisher Scientific) coupled online to the UPLC. The electrospray voltage applied was 2.0 kV. The m/z scan range was 350 to 1800 for full scan, and intact peptides were detected in the Orbitrap at a resolution of 70,000. Peptides were then selected for MS/MS using NCE setting as 28 and the fragments were detected in the Orbitrap at a resolution of 17,500. A data-dependent procedure that alternated between one MS scan followed by 20 MS/MS scans with 15.0 s dynamic exclusion. Automatic gain control (AGC) was set at 5 × 10 4 . Fixed first mass was set as 100 m/z.

Database Search
The resulting MS/MS data were processed using Maxquant search engine (v.1.5.2.8). Protein identification was performed using Sequest HT engine against the UniprotKB Medicago truncatula database (update to June 9th, 2016, including 70,050 protein sequences) [44]. Trypsin/P was specified as cleavage enzyme allowing up to 2 missing cleavages. The mass tolerance for precursor ions was set as 20 ppm in First search and 5 ppm in Main search, and the mass tolerance for fragment ions was set as 0.02 Da. Carbamidomethyl on Cys was specified as fixed modification and oxidation on Met was specified as variable modifications. FDR was adjusted to <1% and minimum score for peptides was set >40.

Validation of Differentially Expressed Unigenes and Proteins by Quantitative Real-Time PCR (qRT-PCR) Analysis
Total RNA extraction, cDNA synthesis, and qRT-PCRs were performed as Zhang et al. (2019) [45] described. Briefly, 0.1 g fresh tissues were used for total RNA extraction by using Trizol reagent (Invitrogen). RNA quality and integrity were checked by Nanodrop 2000 and 0.8% agarose gel. Then, the first strand cDNA were synthesized from 2 µg of total RNA using oligo(dT) 12-18 primer with the cDNA synthesis kit (Fermentas, Burlington, ON, Canada). Primer sequences were listed in Supplementary File 6. The real-time RT-PCR was conducted in ABI7500 with a final volume of 20 mm 3 . The PCR conditions were as follows: 40 cycles of 95 • C denaturation for 5 s, and 52-55 • C annealing and extension for 20 s. The relative expression level of genes for each sample was calculated relative to a calibrator using the DDCT method as described by Livak and Schmittgen (2001) [46].

Statistical Analysis
Statistical analysis was carried out following the ANOVA analysis of variance using SAS for Windows (SAS Institute, Cary, NC, USA). Comparisons between means were carried out using Duncan's multiple range tests at a significance level of p < 0.05.

Conclusions
The schematic figure summarized the alfalfa response to waterlogging stress with indicated difference between the two cultivars ( Figure 8). Compared dataset analysis of proteomics and transcriptomics revealed that 27 and eight genes displayed jointly up-regulated or down-regulated expression profiles at both mRNA and protein levels in M12_W vs. M12_CK comparison, whereas 30 and 27 genes were found to be co-up-regulated or co-down-regulated in M25_W vs. M25_CK comparison, respectively. The strongly enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways for co-up-regulated genes at mRNA and protein levels in M12_W vs. M12_CK comparison were 'Amino sugar and nucleotide sugar metabolism', 'Arginine and proline metabolism', and 'Starch and sucrose metabolism', whereas co-up-regulated protein-related pathways including 'Arginine and proline metabolism' and 'Valine, leucine and isoleucine degradation' were largely enriched in M25_W vs. M25_CK comparison. Importantly, the identified genes related to beta-amylase, Ethylene response Factor (ERF), Calcineurin B-like (CBL) interacting protein kinases (CIPKs), Glutathione peroxidase (GPX), and Glutathione-S-transferase (GST) may play key roles in conferring alfalfa tolerance to waterlogging stress. The present study may contribute to our understanding the molecular mechanism underlying the responses of alfalfa to waterlogging stress, and also provide important clues for further study and in-depth characterization of waterlogging-resistance breeding candidate genes in alfalfa. degradation' were largely enriched in M25_W vs. M25_CK comparison. Importantly, the identified genes related to beta-amylase, Ethylene response Factor (ERF), Calcineurin B-like (CBL) interacting protein kinases (CIPKs), Glutathione peroxidase (GPX), and Glutathione-S-transferase (GST) may play key roles in conferring alfalfa tolerance to waterlogging stress. The present study may contribute to our understanding the molecular mechanism underlying the responses of alfalfa to waterlogging stress, and also provide important clues for further study and in-depth characterization of waterlogging-resistance breeding candidate genes in alfalfa.