- Key Laboratory of Ruminant Molecular and Cellular Breeding, College of Animal Science and Technology, Ningxia University, Yinchuan, China
Wagyu cattle are well-known for their rich marbling. Qinchuan cattle have slower-depositing marbling than Wagyu cattle. However, because of an increase in the consumer demand for high-quality beef and the increasingly stringent standards of beef quality, improving the marbling grade of Qinchuan cattle has become particularly crucial. Therefore, we here considered castrated crossbred Wagyu cattle (crossed with Qinchuan cattle) as the research subjects. Flavor substances in the longissimus dorsi muscle (LDM) of A1 and A5 grades were detected through headspace-solid-phase microextraction-gas chromatography–mass spectrometry (HS-SPME-GC-MS) and electronic nose (E-nose) analysis. Fat deposition-regulating functional genes in both groups were identified through RNA sequencing (RNA-seq) and Weighted gene co-expression network analysis (WGCNA). The results showed that the intramuscular fat (IMF) was significantly higher in A5-grade beef (32.96 ± 1.88) than in A1-grade beef (10.91 ± 1.07) (p < 0.01). In total, 41 and 39 flavor compounds were detected in A1 and A5 grade beef, respectively. Seven aroma compounds were identified base on odor activity values (OAVs) ≥ 1, namely decanal, hexanal, nonanal, heptanol, 1-octen-3-ol, pentanol, and hexanoic acid-methyl ester. Additionally, FABP4, PLIN1, LIPE, ACACA, and CIDEA were the key genes primarily involved in cholesterol metabolism, sterol metabolism, and the PPAR signaling pathway in the two grades of beef. This study attempted to offer comprehensive information on marbling formation-associated candidate genes and gene-enriched pathways, which provides data for future research in beef cattle breeding and beef quality improvement.
1 Introduction
Beef is a nutritious and flavorful meat, and its quality and taste determine consumer requirements (1). Beef quality traits are evaluated based on shear, cooking loss, backfat thickness, and marbling. The richness of marbling is the key criterion for meat quality evaluation and is closely related to tenderness and flavor (2). The higher intramuscular fat (IMF), the marbling, tenderness, and flavor are also better in beef (3). The catabolism of key volatile compounds significantly increased in high IMF roast beef, which is consistent with a more intense sensory flavor (4). The flavor, sweetness, tenderness, and juiciness of roasted beef increased with increased marbling, whereas its sourness and astringency decreased (5). Therefore, increasing the IMF content of beef is essential for augmenting its flavor.
In flavor studies, sensory is a method for evaluating flavor. At present, the E-nose and HS-SPME-GC-MS are also frequently used to study flavor components and correct subjective judgments made during sensory evaluation. E-nose is a sensing method used for analyzing and discriminating aroma profiles and involves 10 sensors, including W1C, W5S, W3C, W6S, W5C, W1S, WIW, W2S, W2W, and W3S (6, 7). HS-SPME-GC-MS, an analytical method combining SPME and GC-MS, which is often applied in flavor studies of meat (8). In a study, ethyl-acetate, ethyl-propionate, and ethyl-hexanoate were identified in pork by HS-SPME-GC-MS. These volatile compounds were negatively correlated with pork freshness (9). Through HS-SPME-GC-MS, 31 volatile compounds were identified in yak and cattle-yak, and the contents of 2-methyl butanal, 2,3-butanedione, 1-propanediol, and 4-methyl-2-pentanone varied between the two beef type (10). In addition to HS-SPME-GC-MS, aroma thresholds and OAVs are also crucial for evaluating flavor compounds, with OAVs ≥ 1 considered to indicate key aroma compounds (11). For example, 19 key flavor compounds are involved in the grilling process of lambs (12). Furthermore, 3-methyl butyraldehyde was the key aroma compound detected in roasted chicken meat (13).
The flavor of livestock meat is chiefly influenced by sex, breed, and feeding conditions. However, differences in nutrient composition and meat flavor are observed among the same breeds under uniform feeding conditions. Genetic factors such as feed conversion efficiency and fat deposition capacity are the predominant cause of these differences. As sequencing technologies, various sequencing tools have been used to investigate genetic factors that affect meat quality. Of them, Beef eq aids in rapidly and efficiently detecting meat properties such as screening for differentially expressed genes (DEGs) in muscle, fat, and other tissues (14). RNA-seq analysis unveiled that 14 genes were differentially expressed between distinct marbling traits of longissimus dorsi muscle (LDM) samples from 20 Nellore cattle, and marbling formation was found to be strongly correlated with lipid and myoglobin oxidation (15). According to a functional analysis of DEGs and metabolites, ACACA, PLIN1, and FABP4 were significantly upregulated in Pingliang Red cattle and crossbred Wagyu cattle, and the contents of 3-iodo-I-tyrosine, arachidonic acid, and cis-aconitic acid were higher in crossbred Wagyu cattle. These genes and metabolites were critical regulators and intermediate in lipid oxidation, which increased fat deposition and beef tenderness (16).
Wagyu has a high propensity to accumulate IMF (marbling), with even distribution of marbling in the meat (17). Marbling in Wagyu cattle belongs to high heritability (0.38–0.50) (18). The IMF content is lower in Qinchuan cattle than in Wagyu cattle. Crossbreeding improves animal growth, adaptability, and meat quality (19). For example, Angus × Nellore progeny were heavier at birth than Nellore progeny (20). The shear force was lower, and the myofibrillar fragmentation index was higher in Aberdeen Angus × Nelore than in Nelore (21). However, studies exploring meat flavor and molecular markers affecting marbling deposition in Wagyu × Qinchuan cattle are fewer. Therefore, in this study, the F1 generation of Wagyu × Qinchuan cattle crossbreeding was considered as the research target. The LDM of A1 and A5-grade marbling were selected to detect flavor compounds. The key aroma compounds were identified on the base of OAVs ≥1. RNA-seq was performed on the A1 and A5-grade beef samples to screen DEGs affecting marbling formation. Using RNA-, modules correlating with the key aroma compounds were constructed, and IMF deposition-regulating core genes in beef cattle were screened through protein–protein (PPI) analysis. The present study provides new insights for improving beef marbling and meat flavor to satisfy the increasingly demanding market requirements.
2 Materials and methods
2.1 Sample collection
Thirty healthy castrated crossbred Wagyu cattle (crossed with Qinchuan cattle, age: 28–30 months) with similar weights (680.53 kg ± 30.78 kg) were randomly selected from the farms in the Beijing region, which had maintained uniform feeding management conditions. The cattle were electrocuted and humanely slaughtered after being starved for 12 h. The cattle were slaughtered according to the GBT19477-2018 cattle slaughtering procedures. Professionally certified technicians rated marbling in beef portions collected between the 12th and 13th rib of the left half-carcass from 30 crossbred Wagyu cattle. The marbling was rated according to the Japan Meat Grading Association marbling score standard (22). At the end of grading, the LDM tissues of the A1-grade (n = 3) and A5-grade (n = 3) cattle were randomly collected, vacuum-sealed, and stored at −20°C for estimating IMF, dry matter content, and other indices, or were stored at −80°C for RNA-seq.
2.2 Determination of flavor components
2.2.1 E-nose analysis of flavor components
LDM tissues (3 g) were kept in a 50 mL headspace vial and incubated for 40 min at 25°C. Then, flavor components in these tissues were identified using the PNE3.5 E-nose (PEN3.5 Airsense, Schwerin, Germany). Supplementary Table S1 lists the components of the E-nose. The assay procedure is consistent with that of Guan et al. (23).
2.2.2 HS-SPME-GC-MS analysis of flavor components
Aroma compounds were analyzed using a GC-MS system (GC-MS 2010 plus, SHIMADZU; Kyoto, Japan) equipped with a DB-WAX capillary column (30 m × 0.25 mm × 0.25 μm, Agilent Technologies; Santa Clara, CA, USA). The SPME fiber of 50/30 μm DVB/CAR/PDMS (Supelco, PA, USA) should be aged before the aroma compounds are extracted. The A1 and A5-grade LDM samples (2 g) were placed in a 15 mL headspace bottle. Then, 2-dichlorobenzene (4 μL, 6.42 μg/mL) was added to each sample as an internal standard. After the mixture was vortexed, the vessel was sealed with a PTFE membrane and placed in a water bath at 60°C for 20 min. Subsequently, the SPME fibers were inserted into a sealed extraction vial for adsorption and extraction for 20 min and immediately transferred to the gas chromatograph (GC) inlet for 5 min of desorption at 250°C. The GC conditions were as follows: the initial temperature was maintained at 40°C for 3 min, ramped to 200°C at a rate of 5°C/min, and again ramped to 230°C at a rate of 10°C/min. Helium was used as the carrier gas at a 2 mL/min flow rate, and the front inlet temperature was 250°C. The MS conditions were as follows: the ion source was electron impact at 70 eV, and its temperature was 230°C, the mass spectrometry (MS) transfer line temperature was 280°C, the solvent delay was 3 min, and a full-scan mode was adopted across 50 to 500 m/z.
2.2.3 Content analysis of flavor components and OAVs
The flavor components were identified by referring to the NIST14 database, retention index (RI) reference values, and authentic volatile standards (24). The content of each volatile in the varying beef grades was calculated based on the o-dichlorobenzene peak area at a known mass concentration. The OAVs were employed to evaluate the contribution of volatile aroma components to the overall beef flavor. OAVs ≥ 1 were considered to indicate the key aroma components for the A1 and A5-grade beef (25). Supplementary Table S2 lists the calculation methods for flavor compounds and key aroma components.
2.3 Measurements of beef quality traits
Using a straightedge, backfat thickness was measured at the 10th and 11th ribs of the carcass. The eye muscle area was the LDM cross-sectional area at the penultimate 1 and 2 thoracic vertebrae positions of the carcass. The eye muscle area was traced using sulfate paper and the cross-sectional area was calculated as follows: height (cm) × width (cm) × 0.7.
2.3.1 Dry matter content analysis
The weighing flask was placed at 105°C until a constant weight was attained and weighed using an electronic weighing scale. This weight was recorded as W0. Then, the meat sample was completely defrosted and chopped into mince. Subsequently, 3 g of the minced meat sample was weighed into a weighing flask and recorded as W1, after which they were placed at 105°C until a constant weight was attained and weighed. This weight was recorded as W2. The dry matter content was calculated using Carvalho et al.’s method (26), as follows: (W2 – W0)/W1 × 100%.
2.3.2 IMF content analysis
The IMF content in the LDM was extracted using the Soxhlet extraction method (27). The meat sample was dried at 105°C and weighed accurately at approximately 0.7 g. This weight was recorded as W. The sample was then wrapped in a filter paper and baked at 105°C until a constant weight was achieved, cooled, and weighed. This weight was recorded as W1. The filter paper was placed in a Soxhlet extractor and extracted with anhydrous acetaldehyde for 8 h. Then, the paper was removed, dried at 105°C for 2 h, and weighed. This weight was recorded as W2. The crude fat content was calculated as follows: (W1-W2)/W × 100%.
2.3.3 Crude protein analysis
The crude protein content was analyzed using a fully automatic Kjeldahl nitrogen determination instrument (Kjeltec, FOSS). The determination standard was referred to as GB/T 9695.11-2008.
2.3.4 pH24h analysis
After the LDM was refrigerated at 4°C for 24 h, the electrode tip of an acidimeter was inserted into the center of meat samples of different grades and allowed to stand for 5 min, and then, the data were read. The pH24h of each sample was determined three times. The pH meter was calibrated with pH 4.0 and 7.0 buffers at 4°C before use. During measurements, the pH meter was recalibrated using standard buffers that corresponded approximately to the temperature of the muscle to compensate for the effect of temperature on pH readings.
2.3.5 Cooking loss analysis
The cooking loss was analyzed according to Musundire et al.’s method (28). First, 100 g of meat samples of different grades were weighed using an electronic scale, and this weight was recorded as W1. Then, all meat samples were simultaneously kept in a thermostatic water bath at 90°C for 45 min, removed, cooled to room temperature, dried, and weighed. This weight was recorded as W2. The cooking loss was calculated as follows: (W1 – W2)/W1 × 100%.
2.4 mRNA preparation and sequencing
RNA was isolated from A1 (n = 3) and A5 (n = 3) grade LDM by Trizol (Invitrogen, Carlsbad, CA, USA). The integrity of the isolated RNA was examined using a Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA). Magnetic beads with oligo (dT) were used to enrich mRNA after total RNA quality was assessed. Using purified mRNA as a template, single-stranded cDNA was synthesized with a random hexamer primer. Then, DNA polymerase I and RNase H were added to synthesize double-stranded cDNA. Using the Illumina HiSeq2000 platform, the cDNA libraries were bipartite sequenced by Chidio Biotechnology Co., Ltd. (Guangzhou, China).
2.5 Sequencing data analysis
After the data were sequenced, the raw data were filtered using FastQc v0.1. The clean data were mapped to the bovine reference genome (Bos_taurus.ARS_UCD1.2.new.genome.fa) by using Hisat2 v2.2.1 (29). StringTie v2.1.2 was used to assemble the mapped results (30). The gene expression levels of A1 and A5 groups were compared using DESeq2 v1.20, with |Fold Change| ≥ 2 and p < 0.05 established as thresholds of significant difference between the groups. The functional and pathway enrichments of DEGs were analyzed using clusterProfiler v4.0.0 (31). p < 0.05 denoted significant enrichment.
2.6 WGCNA analysis
Coexpression modules were built using WGCNA R v1.69 (32). Thresholds were determined using “PickSoft Threshold,” and genes were clustered using TOM values. Genes having similar expression trends were grouped into modules, each containing at least 50 genes. Similar modules were merged with a threshold of 0.25, and genes in coexpressed modules were then identified. Pearson correlation between gene expression and key flavor substances was further determined. The correlation between gene significance (GS) and module membership (MM) was analyzed, with GS > 0.9 and MM > 0.9 being considered as thresholds for identifying the hub genes in the modules.
2.7 PPI analysis
PPI networks were constructed using STRING v11.5 and visualized using Cytoscape v3.8.0 (33). The core genes were screened using MCODE and CytoHubba plugins in Cytoscape.
2.8 RT-qPCR analysis
In total of 10 DEGs were selected to verify the accuracy of RNA-seq results through RT-qPCR. Total RNAs were extracted and assayed for concentration, and then were reverse transcribed into cDNA according to PrimeScript™ RT reagent Kit (Takara, Kyoto, Japan) with gDNA Eraser to remove genomic DNA. The first step of reverse transcription including 2 μL 5 x gDNA Eraser Buffer, 1 μL gDNA Eraser, 2 μL total RNAs, and 5 μL RNase-Free ddH2O. Then, water bath at 42°C for 2 min. After that, 1 μL PrimeScript RT Enzyme Mix I, 1 μL RT Primer, 4 μL 5 x PrimeScript Buffer 2, and 4 μL RNase Free ddH2O were mixed with the reaction solution in the first step. The reaction procedure was 37°C for 15 min, and then 85°C for 5 s. RT-qPCR was performed according to a previous method of our laboratory (34). Supplementary Table S3 lists all primer sequences. All samples contained 3 biological replicates and 3 technical replicates. All data were expressed as Mean ± SEM.
3 Results
3.1 Meat quality traits analysis
Differences in the eye muscle area and cooking loss between the A1 and A5-grade beef were nonsignificant (p > 0.05). Backfat thickness and IMF were higher in the A5-grade beef than in the A1-grade beef (p < 0.01). Crude fat increased significantly in the A5-grade beef than in the A1-grade beef (p < 0.05). By contrast, dry matter content, crude protein, and pH24h were lower in the A5-grade beef than in the A1-grade beef (p < 0.05) (Table 1).
3.2 E-nose analysis of flavor components between two marbling groups
E-nose radar fingerprinting revealed that the W3S, W1C, W3C, W6S, and W5C sensors displayed lower responses in both A1 and A5-grade beef samples, which indicated that alkanes, ammonia, aromatic compounds, and olefins exert less influence on meat quality. By contrast, the sensors W5S, W1S, W1W, W2S, and W2W displayed increased responses in both A1 and A5-grade beef samples. Furthermore, the responses of W1W, W2W, W1S, and W2S were higher in the A5-grade beef than in the A1-grade beef. This indicated that the content of alcohol compounds was higher in the A5-grade beef (Figure 1a).
The PCA analysis was performed using the E-nose data at 100 s for A1 and A5-grade beef. The results demonstrated that the PC1 and PC2 total contribution rate reached 99.2%, thereby exhibiting good ability to distinguish volatile compounds in the A1 and A5-grade beef. According to the PCA, W1W, W3C, W1C, and W1S were correlated with the A5-grade beef, and W5S was correlated with the A1-grade beef (Figure 1b).
3.3 HS-SPME-GC-MS analysis of flavor components
Through HS-SPME-GC-MS, 55 flavor compounds were evaluated in the A1 and A5-grade beef, including 9 aldehydes, 14 alcohols, 4 ketones, 7 acids, 2 alkenes, 9 alkanes, 6 esters, and 4 other compounds. Of them, 41 flavor compounds were identified in the A1-grade beef and 39 were identified in the A5-grade beef (Table 2). Based on OAVs ≥ 1, both A1 and A5-grade beef were found to contain 7 key aroma compounds, including three aldehydes, three alcohols, and one ester. Additionally, the contents of decanal, hexanal, heptanol, and 1-octen-3-ol were higher in the A5-grade beef than in the A1-grade beef (Table 3).
3.4 Evaluation of RNA-seq data
To determine the genes that affect meat quality, we performed RNA-seq analysis of the LDM of the A1 and A5-grade beef. Supplementary Table S4 lists the specific information of the RNA-seq data. Furthermore, compared with the A1 group, 297 DEGs increased and 62 DEGs decreased in the A5 group (p < 0.05) (Figure 2a; Supplementary Table S5). Brown fat cell differentiation and lipid droplet were the key enriched GO terms (Figure 2b; Supplementary Table S6). Moreover, DEGs were primarily enriched in AMPK and pentose phosphate signaling pathways (Figure 2c; Supplementary Table S7). The PPI results unveiled that FABP4, PLIN1, LIPE, ACACA, LEP, and CIDEA were the hub genes (Figure 2d). To verify the accuracy of the RNA-seq results, 10 DEGs were validated through RT-qPCR. Based on the results, DEGs expression levels were consistent with those in the RNA-seq results (Figure 2e).

Figure 2. Enrichment analysis of DEGs. (a) Volcano map of DEGs (the red dot indicates up-regulated DEGs). (b) GO terms analysis of DEGs (BP indicates a biological process, CC indicates cellular component, MF indicates molecular function). (c) KEGG pathway analysis of DEGs. (d) PPI analysis of DEGs. e: RT-qPCR verification of DEGs. Data were represented as Mean ± SE. n = 3.
3.5 WGCNA analysis
For WGCNA, 0.9 was the correlation coefficient threshold, and 12 was the soft threshold power (Figure 3a). Twenty co-expression modules were constructed. The largest module (light yellow) contained 1962 genes, whereas the smallest module (light cyan) contained 79 genes (Figure 3b). The heatmap displayed that these modules were mutually independent (Figure 3c). The constructed gene co-expression modules were linked to the aroma components. According to the results, the turquoise, blue, and yellow modules were related to nonanal (r = −0.97, p = 0.002), 1-octen-3-ol (r = 0.95, p = 0.004), and hexanal (r = 0.83, p = 0.04), respectively (Figure 3d).

Figure 3. Correlation analysis of modules and traits. (a) Screening of optimal soft thresholds. (b) Co-expressed gene modules clustering tree and modules delineation. (c) Co-expression gene modules correlation heat map. (d) Heatmap of key aroma compounds and module correlation.
The turquoise, blue, and yellow modules, respectively, contained 325, 1938, and 155 genes (Supplementary Table S8). The GO term analysis demonstrated that genes in the turquoise module were significantly enriched in the mitochondrial small ribosomal subunit and mitochondrial matrix (Figure 4a). Furthermore, the 325 genes participated in insulin resistance, apoptosis, and FOXO signaling pathways (Figure 4b). The 1938 DEGs in the blue module were prominently concentrated in the negative regulation of cell–cell adhesion and actin cytoskeleton (Figure 4c) and principally involved in glucagon signaling pathways (Figure 4d). The 155 DEGs in the yellow module were chiefly involved in GO terms such as mitochondrion organization, NADH dehydrogenase complex assembly, oxidative phosphorylation, and actin cytoskeleton (Figure 4e). KEGG revealed that these DEGs were also mainly enriched in oxidative phosphorylation, NOD-like receptor, and FOXO signaling pathways (Figure 4f).

Figure 4. Functional enrichment analysis of DEGs in modules. (a) GO analysis of DEGs in the turquoise module. (b) KEGG analysis of DEGs in the turquoise module. (c) GO analysis of DEGs in the blue module. (d) KEGG analysis of DEGs in the blue module. (e) GO analysis of DEGs in the yellow module. (f) KEGG analysis of DEGs in the yellow module.
3.6 Functional analysis of the hub gene
In this study, hub genes in the turquoise, blue, and yellow modules were screened using GS > 0.9 and MM > 0.9. In total, 8 hub genes were identified in the yellow module. These 8 genes exhibited a correlation of 0.17 (p = 0.035) between GS and MM (Figure 5a). A total of 277 hub genes were identified in the blue module, and the correlation was 0.81 (p = 1e-200) (Figure 5b). In total, 66 hub genes were detected in the turquoise module, and the correlation was 0.78 (p = 1e-67) (Figure 5c). Supplementary Table S9 presents the hub genes in the three modules. According to the enrichment analysis, the hub genes primarily participated in cholesterol and sterol metabolism (Figure 5d). The KEGG results unveiled that the hub genes were predominantly concentrated in lipid synthesis and lipid metabolism signaling pathways, such as MAPK and PPAR signaling pathways (Figure 5e). The PPI analysis revealed that FABP4, PLIN1, PPARA, VASP, MSN, ACTN1, TLN1, and CD34 were the core genes in DEGs (Figure 5f).

Figure 5. Hub gene screening and functional analysis. (a) Genes scatterplot in the yellow module. (b) Genes scatterplot in the blue module. (c) Genes scatterplot in the turquoise module. (d) GO enrichment of hub genes. (e) KEGG analysis of hub genes. (f) PPI analysis of hub genes.
4 Discussion
The beef flavor is determined by its texture and aroma (35), with volatile flavor components contributing to aroma component formation (36). These components are formed through lipid oxidation, Maillard reaction, and thermal degradation (37). The key by-products of lipid oxidation are aldehydes, alcohols, and acidic compounds, including 1-octen-3-ol, nonanal, octanal, and hexanal (38). Pyrazines and furans are mainly produced through Maillard reactions, and thermal degradation and lipid oxidation contribute to flavor compound formation from non-volatile water-soluble lipids (39). Several volatile flavor compounds have been detected in livestock, such as alcohols, aldehydes, acids, hydrocarbons, and heterocyclic compounds (40). The IMF content increased significantly with an increase in beef grade (41, 42). Similarly, the content of volatile flavor compounds increased (43), including alcohols and aldehydes with a fresh flavor, ketones with an oily flavor, and esters with a milky flavor. Most alcohol compounds have specific flavors, such as 1-octen-3-ol has a mushroom flavor, octanol has a lemon flavor, and nonanol has a grassy flavor. Octanal, nonanal, and decanal have a sweet flavor, fatty flavor, and sweet flavor, respectively (44). Arginine, citrate, glucose, propionate, 3-hydroxybutyrate, and lipids are correlated to marbling in crossbred Wagyu cattle (45). Propionic acids were converted to glucose through the TCA cycle. This glucose is then available for fatty acid synthesis, which results in IMF deposition (18). The fatty content in lamb was positively correlated with flavor. Additionally, flavor compounds, such as aldehydes and alcohols, IMF content, juiciness, and tenderness were significantly higher in pork with carcasses of higher quality grades than in pork with carcasses of lower grades (46). In the present study, the IMF content was significantly higher in the A5-grade beef (32.96 ± 1.88) than in the A1-grade beef (10.91 ± 1.07). In addition, the aldehydes, ketones, alcohols, aromatic, acids, and lipid compounds were rich between A5 and A1 grade beef. Further, aldehydes, ketones, alcohols, acids, and lipids were greater in the A5-grade beef. This was concordant with the aforementioned findings, and the difference in the content of these flavor substances was the key reason for the richer marbling and IMF content of the A5-grade beef.
Meat mostly attains its flavor from volatile compounds, which rely on the IMF content (47). Fats were precursors for flavor substance formation. These substances could produce aldehydes, ketones, acids, and alcohol compounds through hydrolysis, pyrolysis, oxidation, and the Maillard reaction (48). In addition to feeding, breeding, and sex, key genes and molecular regulatory pathways are crucial factors influencing IMF deposition. RNA-seq revealed that the core genes affecting marbling and IMF content in Nellore cattle were chiefly PLIN1, CISH, UFM1, and TSHZ1 (49). FABP4, TPI1, ACTA1, and MDH2 were highly expressed in marbling-rich meat in Korean cattle (50).
The present study results were similar and revealed that FABP4, PPARG, ACACA, and PLIN1 were highly expressed in the A5-grade marbling beef of the crossbred Wagyu cattle. PLIN1 is a member of the lipid family that augments lipid formation (51). When adipocytes begin to catabolism, PLIN1 modulates the activity of hydrolytic enzymes in lipid droplets to complete hydrolysis (52). Additionally, PLIN1 was strongly expressed in subcutaneous adipose tissues, with its SNP being related to IMF content and thoracic depth in Qinchuan cattle (53). FABP4 is an intracellular lipid chaperone abundantly expressed in adipocytes and macrophages, which can modulate lipid fluxes, transportation, esterification, and β-oxidation and regulate the lipid signal transduction and metabolism (54). Lipids are oxidized to produce hydroperoxides, which are then further oxidized to volatile lipid products by enzyme catalysis, including 3-hexenal, propionaldehyde, and trans-2-hexenol (55). Furthermore, 3-iodo-L-tyrosine, 2,6-diamino hexanoic acid, cis-aconitate, and arachidonic acid were positively associated with marbling, unsaturated, and tenderness (16). In lipogenesis, FABP4 participated in the PPAR signaling pathway as an up-regulated protein, which regulates lipogenesis in human skeletal muscle cells (56). PLIN1 was highly expressed in porcine adipose tissue and was significantly enriched in the PPAR signaling pathway (57). In summary, these genes and metabolites were pivotal factors involved in marbling formation. They are essential for enhancing fat deposition and marbling richness in beef muscles.
In this study, volatile flavor compound-related three modules were determined through WGCNA. The yellow module was positively and notably correlated with hexanal, and the blue module was significantly and positively correlated with 3-octanol and methyl hexanoate. However, the turquoise module exhibited a significantly negative connection with nonanal and pentanol. The modules were screened by determining gene-to-module correlation and gene significance. The results revealed that 8, 277, and 66 hub genes were obtained in the yellow, blue, and turquoise modules, respectively, and primarily included PPARG, CD34, and FABP4. Hub genes were concentrated in the MAPK, cGMP-PKG, and PPAR signaling pathways, and were involved in cholesterol metabolism. PPAR is a major pathway that regulates lipid metabolism, adipogenesis, energy homeostasis, cell growth, and differentiation (58). PPARG is a core regulatory gene of the PPAR signaling pathway. It was a major regulator of adipocyte differentiation (59), and a key regulator of lipid metabolism in adipocytes (60). Improving LPL, FABP4, and PLIN1 expression activates PPARG, thereby driving fat deposition (61). The MAPK signaling pathway is among the major intracellular pathways for muscle development and adipogenesis (62). In a study, Qiangguyin inhibits adipogenic differentiation through the p38 MAPK signaling pathway, thereby reducing fat accumulation in OVX mice (63). Altogether, the genes FABP4, PLIN1, and PPARG directly or indirectly regulate IMF metabolism through molecular pathways and therefore influence the content of volatile flavor compounds. However, since the samples we collected could only satisfy the minimum number of biological replicates for RNA-seq. The results still need to be further confirmed in subsequent experiments. Meanwhile, the complex mechanism of marbling deposition at the cellular level must also be further verified.
5 Conclusion
In summary, the present study preliminarily elucidated flavor differences in different marbling grade beef and relevant genes affecting IMF deposition. Regarding flavor indicators, the contents of key aroma substances such as heptanol, 1-octen-3-ol, hexanoic acid, and methyl ester were higher in the A5-grade beef than in the A1-grade beef. Moreover, FABP4, PLIN1, PPARG, and ACTN1 were potential candidate genes for regulating IMF deposition. These genes were enriched in cholesterol metabolism, cGMP-PKG, MAPK, and PPAR signaling pathways. The present study unveils the reasons for flavor differences in beef of different grades at the molecular level. The finding provides deeper insights into molecular regulatory mechanisms occurring during beef marbling.
Data availability statement
The data presented in the study are deposited in the NCBI repository, accession number PRJNA1009594.
Ethics statement
Animal experiments were performed under the guidelines of the China Animal Care Committee, and agreed upon by the Animal Protection Society of Ningxia University (Permit No. NXUC20200908). The studies were conducted in accordance with the local legislation and institutional requirements. Written informed consent was obtained from the owners for the participation of their animals in this study.
Author contributions
YD: Formal analysis, Writing – original draft, Methodology, Validation. YZ: Data curation, Software, Writing – review & editing. XZ: Methodology, Validation, Writing – review & editing. CL: Formal analysis, Validation, Writing – review & editing. ZS: Data curation, Formal analysis, Writing – review & editing. JX: Methodology, Software, Writing – review & editing. CQ: Data curation, Formal analysis, Writing – review & editing. YM: Supervision, Funding acquisition, Writing – review & editing. YS: Funding acquisition, Writing – review & editing. XK: Writing – review & editing.
Funding
The author(s) declare that financial support was received for the research and/or publication of this article. This research was funded by the National Science and Technology Major Special Project- Science and Technology Innovation 2030 Special Project (2023ZD0404803-02), the National Natural Science Foundation of China (32160776, 32460812), the earmarked fund for the China Agriculture Research System (CARS-36), the Natural Science Foundation of Ningxia Province, China (Grant No. 2021AAC03027), and the West Light Foundation of the Chinese Academy of Sciences (Grant No. XAB2022YW11).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those 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.
Supplementary material
The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fvets.2025.1501177/full#supplementary-material
References
1. Ortega, DL, Hong, SJ, Wang, HH, and Wu, L. Emerging markets for imported beef in China: results from a consumer choice experiment in Beijing. Meat Sci. (2016) 121:317–23. doi: 10.1016/j.meatsci.2016.06.032
2. Kazala, EC, Lozeman, FJ, Mir, PS, Laroche, A, Bailey, DR, and Weselake, RJ. Relationship of fatty acid composition to intramuscular fat content in beef from crossbred wagyu cattle. J Anim Sci. (1999) 77:1717–25. doi: 10.2527/1999.7771717x
3. Mannen, HJ. Identification and utilization of genes associated with beef qualities. Anim Sci J. (2011) 82:1–7. doi: 10.1111/j.1740-0929.2010.00845.x
4. Frank, D, Kaczmarska, K, Paterson, J, Piyasiri, U, and Warner, R. Effect of marbling on volatile generation, oral breakdown and in mouth flavor release of grilled beef. Meat Sci. (2017) 133:61–8. doi: 10.1016/j.meatsci.2017.06.006
5. Frank, D, Ball, A, Hughes, J, Krishnamurthy, R, Piyasiri, U, Stark, J, et al. Sensory and flavor chemistry characteristics of Australian beef: influence of intramuscular fat, feed, and breed. J Agric Food Chem. (2016) 64:4299–311. doi: 10.1021/acs.jafc.6b00160
6. Qiu, SS, and Wang, J. The prediction of food additives in the fruit juice based on electronic nose with chemometrics. Food Chem. (2017) 230:208–14. doi: 10.1016/j.foodchem.2017.03.011
7. Cui, SQ, Wu, JF, Wang, J, and Wang, XL. Discrimination of American ginseng and Asian ginseng using electronic nose and gas chromatography-mass spectrometry coupled with chemometrics. J Ginseng Res. (2017) 41:85–95. doi: 10.1016/j.jgr.2016.01.002
8. Ahamed, Z, Seo, J-K, Eom, J-U, and Yang, H. Optimization of volatile compound extraction on cooked meat using HS-SPME-GC-MS, and evaluation of diagnosis to meat species using volatile compound by multivariate data analysis. Agric Food Sci. (2023) 188:115374. doi: 10.1016/j.lwt.2023.115374
9. Sun, Y, Fu, M, Li, Z, and Peng, X. Evaluation of freshness in determination of volatile organic compounds released from pork by HS-SPME-GC-MS. Food Anal. Methods. (2018) 11:1321–9. doi: 10.1007/s12161-017-1109-6
10. Zhang, B, Cao, M, Wang, X, Guo, S, Ding, Z, Kang, Y, et al. The combined analysis of GC-IMS and GC-MS reveals the differences in volatile flavor compounds between yak and cattle-yak meat. Foods. (2024) 13:2364. doi: 10.3390/foods13152364
11. Wang, Q, Zhang, Q, Liu, K, An, J, Zhang, S, Chen, Q, et al. Optimization of solid-state fermentation technology and analysis of key aroma components of compound rice wine. Food Sci Technol Res. (2022) 28:35–43. doi: 10.3136/fstr.FSTR-D-21-00094
12. Cheng, K, Liu, T, Yang, C, Yang, H, and Liu, D. Relationship between phospholipid molecules species and volatile compounds in grilled lambs during the heating process. Food Chem X. (2024) 21:101113. doi: 10.1016/j.fochx.2023.101113
13. Yu, Y, Wang, G, Yin, X, Ge, C, and Liao, G. Effects of different cooking methods on free fatty acid profile, water-soluble compounds and flavor compounds in Chinese Piao chicken meat. Food Res Int. (2021) 149:110696. doi: 10.1016/j.foodres.2021.110696
14. Valdés-Hernández, J, Ramayo-Caldas, Y, Passols, M, Sebastià, C, Criado-Mesas, L, Crespo-Piazuelo, D, et al. Global analysis of the association between pig muscle fatty acid composition and gene expression using RNA-Seq. Sci Rep. (2023) 13:535. doi: 10.1038/s41598-022-27016-x
15. Muniz, MMM, Fonseca, LFS, dos Santos Silva, DB, Magalhães, AFB, Ferro, JA, Chardulo, LAL, et al. Small genetic variation affecting mRNA isoforms associated with marbling and meat color in beef cattle. Funct Integr Genomics. (2022) 22:451–66. doi: 10.1007/s10142-022-00844-w
16. Wang, Z, An, X, Yang, Y, Zhang, L, Jiao, T, Zhao, SJ, et al. Comprehensive analysis of the longissimus dorsi transcriptome and metabolome reveals the regulatory mechanism of different varieties of meat quality. J Agric Food Chem. (2023) 71:1234–45. doi: 10.1021/acs.jafc.2c07043
17. Gotoh, T, Nishimura, T, Kuchida, K, and Mannen, H. The Japanese wagyu beef industry: current situation and future prospects – a review. Asian Austral J Anim. (2018) 31:933–50. doi: 10.5713/ajas.18.0333
18. Park, SJ, Beak, SH, Jung, DJS, Kim, SY, Jeong, IH, Piao, MY, et al. Genetic, management, and nutritional factors affecting intramuscular fat deposition in beef cattle – a review. Asian Austral J Anim. (2018) 31:1043–61. doi: 10.5713/ajas.18.0310
19. Hegde, NG. Impact of crossbreeding and upgrading of nondescript cattle and buffaloes on livestock quality and income. Indian J Anim Sci. (2017) 88:606–11. doi: 10.56093/ijans.v88i5.80009
20. Favero, R, Menezes, GO, Torres, R, Silva, L, Bonin, M, Feijó, G, et al. Crossbreeding applied to systems of beef cattle production to improve performance traits and carcass quality. Animal. (2019) 13:2679–86. doi: 10.1017/S1751731119000855
21. Oliveira, P, Oliveira, M, Bonin, M, Avalo, S, Cancio, P, Nascimento, J, et al. Carcass and meat characteristics of feedlot finished nelore cattle and their crossbreeds in the Brazilian Pantanal. Livestock Sci. (2021) 244:104360. doi: 10.1016/j.livsci.2020.104360
22. Motoyama, M, Sasaki, K, and Watanabe, AJMS. Wagyu and the factors contributing to its beef quality: a Japanese industry overview. Meat Sci. (2016) 120:10–8. doi: 10.1016/j.meatsci.2016.04.026
23. Guan, C, Liu, T, Li, Q, Wang, D, and Zhang, Y. Analyzing the effect of baking on the flavor of defatted tiger nut flour by E-tongue, E-nose and HS-SPME-GC-MS. Food Secur. (2022) 11:446. doi: 10.3390/foods11030446
24. Liu, H, Wang, Z, Zhang, D, Shen, Q, Hui, T, and Ma, JJFRI. Generation of key aroma compounds in Beijing roasted duck induced via Maillard reaction and lipid pyrolysis reaction. Food Res Int. (2020) 136:109328. doi: 10.1016/j.foodres.2020.109328
25. Xiao, X, Hou, C, Zhang, D, Li, X, Ren, C, Ijaz, M, et al. Effect of pre-and post-rigor on texture, flavor, heterocyclic aromatic amines and sensory evaluation of roasted lamb. Meat Sci. (2020) 169:108220. doi: 10.1016/j.meatsci.2020.108220
26. Carvalho, R, Boari, C, Villela, S, Pires, A, Mourthé, M, Oliveira, F, et al. Differences between sexes, muscles and aging times on the quality of meat from wagyu× Angus cattle finished in feedlot. Animal Product Sci. (2016) 58:350–7. doi: 10.1071/AN15804
27. Han, F, Li, J, Zhao, R, Liu, L, Li, L, Li, Q, et al. Identification and co-expression analysis of long noncoding RNAs and mRNAs involved in the deposition of intramuscular fat in Aohan fine-wool sheep. BMC Genomics. (2021) 22:1–14. doi: 10.1186/s12864-021-07385-9
28. Musundire, M, Halimani, T, and Chimonyo, MJBPS. Physical and chemical properties of meat from scavenging chickens and helmeted guinea fowls in response to age and sex. Br Poult Sc. (2017) 58:390–6. doi: 10.1080/00071668.2017.1313961
29. Li, C, Li, S, Yang, C, Wang, X, Zhou, X, Shi, Y, et al. Blood transcriptome reveals immune and metabolic-related genes involved in growth of pasteurized colostrum-fed calves. Front Genet. (2023) 14:1075950. doi: 10.3389/fgene.2023.1075950
30. Pertea, M, Pertea, GM, Antonescu, CM, Chang, TC, Mendell, JT, and Salzberg, SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. (2015) 33:290–5. doi: 10.1038/nbt.3122
31. Wu, T, Hu, E, Xu, S, Chen, M, Guo, P, Dai, Z, et al. clusterProfiler 4.0: a universal enrichment tool for interpreting omics data. Innovation. (2021) 2:100141. doi: 10.1016/j.xinn.2021.100141
32. Langfelder, P, and Horvath, S. WGCNA: an R package for weighted correlation network analysis. Bmc Bioinformat. (2008) 9:9. doi: 10.1186/1471-2105-9-559
33. Shannon, P, Markiel, A, Ozier, O, Baliga, NS, Wang, JT, Ramage, D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. (2003) 13:2498–504. doi: 10.1101/gr.1239303
34. Zhao, L, Ding, Y, Yang, C, Wang, P, Zhao, Z, Ma, Y, et al. Identification and characterization of hypothalamic circular RNAs associated with bovine residual feed intake. Gene. (2023) 851:147017. doi: 10.1016/j.gene.2022.147017
35. Tansawat, R, Maughan, CAJ, Ward, RE, Martini, S, and Cornforth, DP. Chemical characterisation of pasture- and grain-fed beef related to meat quality and flavour attributes. Int J Food Sci Tech. (2013) 48:484–95. doi: 10.1111/j.1365-2621.2012.03209.x
36. Legako, JF, Brooks, JC, O'Quinn, TG, Hagan, TDJ, Polkinghorne, R, Farmer, LJ, et al. Consumer palatability scores and volatile beef flavor compounds of five USDA quality grades and four muscles. Meat Sci. (2015) 100:291–300. doi: 10.1016/j.meatsci.2014.10.026
37. Troise, AD, Fogliano, V, Vitaglione, P, and Berton-Carabin, CC. Interrelated routes between the Maillard reaction and lipid oxidation in emulsion systems. J Agr Food Chem. (2020) 68:12107–15. doi: 10.1021/acs.jafc.0c04738
38. Xia, W, and Budge, SM. Techniques for the analysis of minor lipid oxidation products derived from Triacylglycerols: epoxides, alcohols, and ketones. Compr Rev Food Sci F. (2017) 16:735–58. doi: 10.1111/1541-4337.12276
39. Gardner, K, and Legako, JF. Volatile flavor compounds vary by beef product type and degree of doneness. J Anim Sci. (2018) 96:4238–50. doi: 10.1093/jas/sky287
40. Watkins, PJ, Rose, G, Warner, RD, Dunshea, FR, and Pethick, DW. A comparison of solid-phase microextraction (SPME) with simultaneous distillation-extraction (SDE) for the analysis of volatile compounds in heated beef and sheep fats. Meat Sci. (2012) 91:99–107. doi: 10.1016/j.meatsci.2011.12.004
41. Zakharova, IO, and Avrova, NF. The effect of cold stress on ganglioside fatty acid composition and ganglioside-bound sialic acid content of rat brain subcellular fractions. J Therm Biol. (2001) 26:215–22. doi: 10.1016/S0306-4565(00)00045-0
42. Ripoll, G, Alberti, P, Panea, B, Olleta, JL, and Sañudo, C. Near-infrared reflectance spectroscopy for predicting chemical, instrumental and sensory quality of beef. Meat Sci. (2008) 80:697–702. doi: 10.1016/j.meatsci.2008.03.009
43. Machiels, D, and Istasse, L. Evaluation of two commercial solid-phase microextraction fibres for the analysis of target aroma compounds in cooked beef meat. Talanta. (2003) 61:529–37. doi: 10.1016/S0039-9140(03)00319-9
44. O'Sullivan, MG, Byrne, DV, Nielsen, JH, Andersen, HJ, and Martens, M. Sensory and chemical assessment of pork supplemented with iron and vitamin E. Meat Sci. (2003) 64:175–89. doi: 10.1016/S0309-1740(02)00177-8
45. Connolly, S, Dona, A, Hamblin, D, D’Occhio, MJ, and González, LA. Changes in the blood metabolome of wagyu crossbred steers with time in the feedlot and relationships with marbling. Sci Rep. (2020) 10:18987. doi: 10.1038/s41598-020-76101-6
46. Hoa, VB, Seong, P-N, Cho, S-H, Kang, S-M, Kim, Y-S, Moon, S-S, et al. Quality characteristics and flavor compounds of pork meat as a function of carcass quality grade. Asian Australas J Anim Sci. (2019) 32:1448–57. doi: 10.5713/ajas.18.0965
47. Maughan, C, Tansawat, R, Cornforth, D, Ward, R, and Martini, S. Development of a beef flavor lexicon and its application to compare the flavor profile and consumer acceptance of rib steaks from grass-or grain-fed cattle. Meat Sci. (2012) 90:116–21. doi: 10.1016/j.meatsci.2011.06.006
48. Xu, L, Yu, X, Li, M, Chen, J, and Wang, X. Monitoring oxidative stability and changes in key volatile compounds in edible oils during ambient storage through HS-SPME/GC–MS. Int J Food Propert. (2017) 20:S2926–38. doi: 10.1080/10942912.2017.1382510
49. Fonseca, LFS, Silva, DBD, Gimenez, DFJ, Baldi, F, Ferro, JA, Chardulo, LAL, et al. Gene expression profiling and identification of hub genes in Nellore cattle with different marbling score levels. Genomics. (2020) 112:873–9. doi: 10.1016/j.ygeno.2019.06.001
50. Shin, SC, and Chung, ER. Identification of differentially expressed genes between high and low marbling score grades of the longissimus lumborum muscle in Hanwoo (Korean cattle). Meat Sci. (2016) 121:114–8. doi: 10.1016/j.meatsci.2016.05.018
51. Padilla-Benavides, T, Velez-delValle, C, Marsch-Moreno, M, Castro-Muñozledo, F, and Kuri-Harcuch, W. Lipogenic enzymes complexes and cytoplasmic lipid droplet formation during adipogenesis. J Cell Biochem. (2016) 117:2315–26. doi: 10.1002/jcb.25529
52. Sun, Z, Gong, J, Wu, H, Xu, W, Wu, L, Xu, D, et al. Perilipin1 promotes unilocular lipid droplet formation through the activation of Fsp27 in adipocytes. Nat Commun. (2013) 4:1594. doi: 10.1038/ncomms2581
53. Raza, SHA, Shijun, L, Khan, R, Schreurs, NM, Manzari, Z, Abd El-Aziz, AH, et al. Polymorphism of the PLIN1 gene and its association with body measures and ultrasound carcass traits in Qinchuan beef cattle. Genome. (2020) 63:483–92. doi: 10.1139/gen-2019-0184
54. Walenna, NF, Kurihara, Y, Chou, B, Ishii, K, Soejima, T, Itoh, R, et al. Chlamydia pneumoniae exploits adipocyte lipid chaperone FABP4 to facilitate fat mobilization and intracellular growth in murine adipocytes. Biochem Biophys Res Commun. (2018) 495:353–9. doi: 10.1016/j.bbrc.2017.11.005
55. Shahidi, F, and Hossain, AJM. Role of lipids in food flavor generation. Molecules. (2022) 27:5014. doi: 10.3390/molecules27155014
56. Wang, X-W, Sun, Y-J, Chen, X, and Zhang, W-Z. Interleukin-4-induced FABP4 promotes lipogenesis in human skeletal muscle cells by activating the PPAR γ signaling pathway. Cell Biochem Biophys. (2022) 80:355–66. doi: 10.1007/s12013-022-01063-7
57. Li, B, Weng, Q, Dong, C, Zhang, Z, Li, R, Liu, J, et al. A key gene, PLIN1, can affect porcine intramuscular fat content based on transcriptome analysis. Genes (Basel). (2018) 9:194. doi: 10.3390/genes9040194
58. Patsenker, E, Thangapandi, VR, Beil-Wagner, J, Buch, T, vom Berg, J, Hampe, J, et al. Genetic variant PNPLA3 I148M accelerates fat accumulation in livers of mice with ASH/NASH via damping of PPAR alpha and PPAR gamma signalling pathways. Cactus. (2020) 73:S175. doi: 10.1016/S0168-8278(20)30857-6
59. Ma, X, Wei, D, Cheng, G, Li, S, Wang, L, Wang, Y, et al. Bta-miR-130a/b regulates preadipocyte differentiation by targeting PPARG and CYP2U1 in beef cattle. Mol Cell Probes. (2018) 42:10–7. doi: 10.1016/j.mcp.2018.10.002
60. Su, YQ, Liu, XM, Lian, JM, and Deng, C. Epigenetic histone modulations of PPARγ and related pathways contribute to olanzapine-induced metabolic disorders. Pharmacol Res. (2020) 155:104703. doi: 10.1016/j.phrs.2020.104703
61. Gu, H, Zhou, Y, Yang, JZ, Li, JA, Peng, YX, Zhang, X, et al. Targeted overexpression of PPARγ in skeletal muscle by random insertion and CRISPR/Cas9 transgenic pig cloning enhances oxidative fiber formation and intramuscular fat deposition. FASEB J. (2021) 35:e21308. doi: 10.1096/fj.202001812RR
62. Xiao, YT, Wang, J, Yan, WH, Zhou, KJ, Cao, Y, and Cai, W. p38α MAPK antagonizing JNK to control the hepatic fat accumulation in pediatric patients onset intestinal failure. Cell Death Dis. (2017) 8:8. doi: 10.1038/cddis.2017.523
Keywords: crossbred Wagyu cattle, beef grade, flavor components, marbling, transcriptome
Citation: Ding Y, Zhang Y, Zhou X, Li C, Su Z, Xu J, Qu C, Ma Y, Shi Y and Kang X (2025) Transcriptome and HS-SPME-GC-MS analysis of key genes and flavor components associated with beef marbling. Front. Vet. Sci. 12:1501177. doi: 10.3389/fvets.2025.1501177
Edited by:
Peter Dovc, University of Ljubljana, SloveniaReviewed by:
Chaoyun Yang, Xichang College, ChinaXiaodong He, Ding xi Vocational and Technical College, China
Copyright © 2025 Ding, Zhang, Zhou, Li, Su, Xu, Qu, Ma, Shi and Kang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Xiaolong Kang, a2FuZ3hsOTUyN0AxMjYuY29t
†These authors share first authorship