There was no considerable difference between the AGES and the NCBI databases regarding the proportion of false negative results obtained for samples 1–7. At a threshold of 0.05%, false negative results were only obtained for pig in sample 4. At a threshold of 0.1%, none of the combinations of sample–animal species led to false negative results, with the exception of proportions being close to the threshold (0.1%). The data indicate that at a threshold of ≥0.5%, the probability of obtaining false negative results is very low.

### 3.2.3. Probability of Detection (Qualitative Evaluation)

Figure 7 shows the probability of detection for three thresholds (0.1%, 0.5%, and 1%) and three scenarios, namely a laboratory with average performance, a laboratory with positive bias, and a laboratory with negative bias.

<figure><img src="image_65.png" /><figcaption>**Figure 7.** Probability of detection for a laboratory with average performance, a laboratory with positive bias, and a laboratory with negative bias, depending on the threshold and the database.</figcaption></figure>
Figure 7 indicates that a threshold value of 1% seems to be a good compromise, provided that the variance function determined in the ring trial is equal to the actual variance function. A threshold value of 1% guarantees that a laboratory with a positive bias (overestimating the actual proportion) does not complain if the proportion of a certain animal species is 0.5%, whereas even a laboratory with a negative bias (underestimating the actual proportion) will be able to identify proportions >1.5% reliably.

### 3.3. *Negative Control*

Sample 8 was a DNA extract from maize, serving as a negative control. Since the marker system designed for mammals and poultry species does not detect maize, all reads that were obtained for sample 8 had to be regarded as false positive.

Table 8 lists the number of total reads (after pipeline) per laboratory and species. Laboratories 12 and 13 did not submit results for sample 8, laboratory 07 only provided results for one of the two subsamples.

In total, eight animal species were identified in the negative control by alignment against the AGES database. Fourteen further animal species, including the species *Homo sapiens* were identified, when the NCBI database was used. Per laboratory, up to five animal species were only identified with the NCBI database, with the exception of laboratory 06, which even detected eight additional animal species in subsample A.

In most cases, within a laboratory, the number of reads per animal species was similar for both subsamples. When the AGES database was used for alignment, most reads were assigned to beef and pig. Sample multiplexing in general, together with an inappropriate index layout, carries the risk of index misassignment. This is obviously the reason for the over-represented number of reads for pig and beef in the negative controls, as these animal species represent the main quantities in the samples. Although the index kit was used according to the manufacturer's instructions, the number of reads of these animal species could be reduced to the expected level in a supplementary experiment with an alternative index layout. Alignment against the NCBI database also resulted in considerably high numbers of reads (laboratories 06, 07, and 14 > 300, laboratory 07 even > 1000) for *Homo sapiens*, which was not contained in the AGES database.

As mentioned above, the marker system applied does not detect maize. Thus, the high number of reads for maize obtained by laboratory 06 for one subsample seems to be caused by a random error.

In general, the total number of false reads obtained for both subsamples was similar within a laboratory. Larger differences between the total reads of subsamples was observed for laboratories 03 and 14 (AGES and NCBI databases) and laboratory 06 (only NCBI database).
<table><caption>Table 8. Reads obtained for sample 8 (negative control). (A): AGES database, (B): NCBI database.</caption><thead><tr><th rowspan="3">Species</th><th colspan="15">(A)</th></tr><tr><th colspan="15">Laboratory</th></tr><tr><th>01</th><th>02</th><th>03</th><th>04</th><th>05</th><th>06</th><th>07</th><th>08</th><th>14</th><th>15</th><th>20</th></tr><tr><th>A</th><th>B</th><th>A</th><th>B</th><th>A</th><th>B</th><th>A</th><th>B</th><th>A</th><th>B</th><th>A</th><th>B</th><th>A</th><th>B</th></tr></thead><tbody><tr><td>Chicken</td><td>-</td><td>-</td><td>13</td><td>28</td><td>-</td><td>-</td><td>42</td><td>30</td><td>4</td><td>98</td><td>-</td><td>-</td><td>2</td><td>-</td><td>10</td><td>8</td><td>5</td></tr><tr><td>Horse</td><td>-</td><td>-</td><td>7</td><td>6</td><td>-</td><td>2</td><td>57</td><td>51</td><td>-</td><td>2</td><td>-</td><td>-</td><td>-</td><td>-</td><td>1</td><td>-</td><td>7</td><td>5</td></tr><tr><td>Turkey</td><td>-</td><td>-</td><td>5</td><td>11</td><td>-</td><td>6</td><td>51</td><td>52</td><td>3</td><td>7</td><td>1</td><td>-</td><td>-</td><td>3</td><td>232</td><td>6</td><td>15</td><td>13</td><td>9</td></tr><tr><td>Beef</td><td>10</td><td>9</td><td>25</td><td>22</td><td>2</td><td>171</td><td>17</td><td>65</td><td>101</td><td>97</td><td>222</td><td>-</td><td>-</td><td>43</td><td>51</td><td>185</td><td>148</td><td>57</td><td>99</td><td>110</td><td>55</td></tr><tr><td>Pig</td><td>-</td><td>1</td><td>1</td><td>3</td><td>23</td><td>87</td><td>11</td><td>3</td><td>216</td><td>114</td><td>563</td><td>-</td><td>-</td><td>21</td><td>21</td><td>91</td><td>56</td><td>38</td><td>54</td><td>50</td><td>26</td></tr><tr><td>Sheep</td><td>5</td><td>14</td><td>51</td><td>79</td><td>2</td><td>2</td><td>82</td><td>74</td><td>3</td><td>9</td><td>19</td><td>-</td><td>-</td><td>3</td><td>1</td><td>11</td><td>3</td><td>-</td><td>-</td><td>9</td><td>10</td></tr><tr><td>Goat</td><td>1</td><td>1</td><td>-</td><td>1</td><td>-</td><td>-</td><td>-</td><td>1</td><td>1</td><td>1</td><td>2</td><td>-</td><td>-</td><td>-</td><td>1</td><td>-</td><td>1</td><td>-</td><td>-</td><td>1</td><td>1</td></tr><tr><td>Bison bonasus</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>1</td><td>-</td><td>2</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Total reads</td><td>16</td><td>25</td><td>102</td><td>150</td><td>27</td><td>268</td><td>261</td><td>276</td><td>328</td><td>328</td><td>807</td><td>-</td><td>-</td><td>69</td><td>77</td><td>529</td><td>394</td><td>103</td><td>178</td><td>197</td><td>111</td></tr></tbody></table><table><thead><tr><th rowspan="3">Species</th><th colspan="15">(B)</th></tr><tr><th colspan="15">Laboratory</th></tr><tr><th>01</th><th>02</th><th>03</th><th>04</th><th>05</th><th>06</th><th>07</th><th>08</th><th>14</th><th>15</th><th>20</th></tr><tr><th>A</th><th>B</th><th>A</th><th>B</th><th>A</th><th>B</th><th>A</th><th>B</th><th>A</th><th>B</th><th>A</th><th>B</th><th>A</th><th>B</th><th>A</th><th>B</th><th>A</th><th>B</th></tr></thead><tbody><tr><td>Chicken</td><td>-</td><td>-</td><td>13</td><td>28</td><td>-</td><td>-</td><td>42</td><td>30</td><td>4</td><td>98</td><td>-</td><td>-</td><td>2</td><td>-</td><td>10</td><td>177</td><td>2</td><td>10</td><td>8</td><td>5</td></tr><tr><td>Horse</td><td>-</td><td>-</td><td>7</td><td>6</td><td>-</td><td>2</td><td>57</td><td>51</td><td>-</td><td>2</td><td>-</td><td>-</td><td>-</td><td>-</td><td>1</td><td>-</td><td>-</td><td>-</td><td>7</td><td>5</td></tr><tr><td>Turkey</td><td>-</td><td>-</td><td>5</td><td>11</td><td>-</td><td>6</td><td>51</td><td>50</td><td>3</td><td>7</td><td>1</td><td>-</td><td>-</td><td>3</td><td>231</td><td>6</td><td>6</td><td>15</td><td>13</td><td>9</td></tr><tr><td>Beef</td><td>10</td><td>9</td><td>25</td><td>22</td><td>2</td><td>171</td><td>17</td><td>65</td><td>101</td><td>97</td><td>222</td><td>-</td><td>-</td><td>43</td><td>52</td><td>187</td><td>149</td><td>57</td><td>99</td><td>111</td><td>55</td></tr><tr><td>Pig</td><td>-</td><td>1</td><td>1</td><td>3</td><td>23</td><td>99</td><td>11</td><td>3</td><td>226</td><td>114</td><td>564</td><td>-</td><td>-</td><td>21</td><td>21</td><td>91</td><td>57</td><td>38</td><td>54</td><td>50</td><td>27</td></tr><tr><td>Sheep</td><td>6</td><td>14</td><td>50</td><td>76</td><td>2</td><td>2</td><td>80</td><td>74</td><td>3</td><td>9</td><td>19</td><td>-</td><td>-</td><td>3</td><td>1</td><td>7</td><td>3</td><td>-</td><td>-</td><td>9</td><td>10</td></tr><tr><td>Goat</td><td>1</td><td>1</td><td>-</td><td>1</td><td>-</td><td>-</td><td>1</td><td>1</td><td>1</td><td>1</td><td>2</td><td>-</td><td>-</td><td>-</td><td>-</td><td>1</td><td>-</td><td>-</td><td>-</td><td>-</td><td>1</td></tr><tr><td>Bison bonasus</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>1</td><td>-</td><td>1</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Bos mutus</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>1</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>1</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Brachypodium<br>syzygium</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>1</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Coregonus</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>2</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>migratorius</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Equus zebra</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Eukaryotic<br>synthetic</td><td>-</td><td>-</td><td>-</td><td>15</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>
<table><caption>Table 8. Cont.</caption><thead><tr><th rowspan="3">Species</th><th colspan="15">(B)</th></tr><tr><th colspan="15">Laboratory</th></tr><tr><th colspan="2">01</th><th colspan="2">02</th><th colspan="2">03</th><th colspan="2">04</th><th colspan="2">06</th><th colspan="2">07</th><th colspan="2">08</th><th colspan="2">14</th><th colspan="2">15</th><th colspan="2">20</th></tr><tr><th>A</th><th>B</th><th>A</th><th>B</th><th>A</th><th>B</th><th>A</th><th>B</th><th>A</th><th>B</th><th>A</th><th>B</th><th>A</th><th>B</th><th>A</th><th>B</th><th>A</th><th>B</th><th>A</th><th>B</th></tr></thead><tbody><tr><td>Homo sapiens</td><td>-</td><td>-</td><td>58</td><td>46</td><td>30</td><td>188</td><td>82</td><td>215</td><td>302</td><td>784</td><td>1086</td><td>-</td><td>5</td><td>124</td><td>687</td><td>451</td><td>-</td><td>15</td><td>-</td><td>49</td></tr><tr><td>Meleagris ocellata</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>2</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>1</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Oncorhynchus environmental</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>1</td><td>1</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Oncorhynchus mykiss</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>37</td><td>91</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Ovis ammon</td><td>-</td><td>1</td><td>1</td><td>5</td><td>-</td><td>-</td><td>2</td><td>1</td><td>1</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>4</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Ovis vignei</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>1</td><td>1</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Plascolosoma esculenta</td><td>-</td><td>-</td><td>-</td><td>-</td><td>1</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Synthetic construct</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>20</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Zea mays</td><td>-</td><td>-</td><td>2</td><td>-</td><td>-</td><td>-</td><td>3</td><td>851</td><td>3</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>2</td><td>-</td></tr><tr><td>total reads</td><td>17</td><td>26</td><td>160</td><td>215</td><td>57</td><td>469</td><td>344</td><td>496</td><td>1553</td><td>1208</td><td>1894</td><td>-</td><td>74</td><td>202</td><td>1218</td><td>846</td><td>104</td><td>193</td><td>200</td><td>161</td></tr></tbody></table>
4. Conclusions

In summary, evaluation of data from the interlaboratory ring trial indicates that the DNA metabarcoding method performed on an Illumina platform is applicable for determining the proportion of the seven animal species with the given precision. Furthermore, the applicability of the method for testing foodstuff was demonstrated by the correct identification of the ingredients of a model sausage, which also supports the results in our study published recently.

Based on the data of the ring trial, a threshold of 0.5% is suitable to reliably assess whether a certain animal species is contained in a sample. The DNA metabarcoding method turned out to be rather robust and is therefore suitable to be implemented in routine analysis in official food control laboratories. Even laboratories that did not have much experience in NGS were able to provide reliable results. We suggest strictly following the given protocol. The results of the interlaboratory ring trial indicate that even alternative test kits or various sequencing platforms might be applied. However, the impact of any deviations from the experimental conditions has obviously to be tested before implementation in routine analysis.

Correct index recognition is of particular importance for pooled DNA libraries. We recommend frequently changing the index kits or the use of longer index sequences to avoid false positive and/or false negative results.

For taxonomic assignment, we suggest applying a customized database, as the pipeline is completed significantly faster and no nonsense results from erroneous database entries occur. However, if unexpected read losses and non-identifiable reads occur, the additional use of the entire NCBI database or any other appropriate sequence database is recommended.

In order to increase interlaboratory comparability of results obtained by DNA metabarcoding methods, it would be necessary to establish a reference database with verified sequence entries of relevant species. Access to adequate reference material would also facilitate harmonization of the methods used.

In general, determination of the meat content (w/w) from the number of NGS reads or the determined target DNA concentration is a well-known difficulty, especially in the quantification of meat species in processed foods. The result is also influenced by the degree of processing of the sample present and by the type of animal ingredients used. Data from testing reference samples out of proficiency testing schemes confirm the limitations known for DNA quantification in meat products [12]. Quantitative results should therefore serve only as rough estimates for weight ratios of different species in food.

**Supplementary Materials:** The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/foods11081108/s1, Supplementary Table S1: Sequences included into the reference database, Supplementary Table S2: Total number of reads after pipeline (n = 14, samples 1–7, two subsamples each). min: minimal value, max: maximal value, mean: arithmetic mean, RSD: relative standard deviation, Supplementary Table S3: Recovery (%) (total number of reads after pipeline related to the number of raw reads before analysis pipeline). (n = 14, samples 1–7, two subsamples each). min: minimal value, max: maximal value, mean: arithmetic mean, RSD: relative standard deviation.

**Author Contributions:** Conceptualization, R.H.; methodology, S.D. and R.H.; validation, S.U., K.F., A.S., K.N. and K.S.; writing—original draft preparation, M.C.-M., S.D. and R.H.; writing—review and editing, R.H., M.C.-M., S.D., S.U., K.F., A.S.; visualization, M.C.-M., K.F. and A.S.; supervision, R.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** The ring trial was supported by the Federal Office of Consumer Protection and Food Safety, Germany (BVL). BVL funded all reagents for performing the tests and the statistical evaluation by QuoData GmbH.

**Institutional Review Board Statement:** Ethical review and approval was waived for this study as no live animals were used or slaughtered to achieve the aims of the study.
**Data Availability Statement:** The datasets generated during the current study are available from the corresponding authors upon reasonable request.

**Acknowledgments:** We thank the members of the working group “NGS Species identification” within the scope of the official method collection according to §64 of the German Food and Feed Code (LFGB) for their participation in the study.

**Conflicts of Interest:** The authors declare no conflict of interest.

## References

1. Ballin, N.Z.; Vogensen, F.K.; Karlsson, A.H. Species determination—Can we detect and quantify meat adulteration? *Meat Sci.* **2009**, *83*, 165–174. [CrossRef] [PubMed]

2. Ballin, N.Z. Authentication of meat and meat products. *Meat Sci.* **2010**, *86*, 577–587. [CrossRef] [PubMed]

3. Sajali, N.; Wong, S.C.; Abu Bakar, S.; Khairil Mokhtar, N.F.; Manaf, Y.N.; Yuswan, M.H.; Mohd Desa, M.N. Analytical approaches of meat authentication in food. *Int. J. Food Sci. Technol.* **2021**, *56*, 1535–1543. [CrossRef]

4. Kumar, A.; Kumar, R.R.; Sharma, B.D.; Gokulakrishnan, P.; Mendiratta, S.K.; Sharma, D. Identification of species origin of meat and meat products on the DNA basis: A review. *Crit. Rev. Food Sci. Nutr.* **2015**, *55*, 1340–1351. [CrossRef] [PubMed]

5. Köppel, R.; Ruf, J.; Zimmerli, F.; Breitenmoser, A. Multiplex real-time PCR for the detection and quantification of DNA from beef, pork, chicken and turkey. *Eur. Food Res. Technol.* **2008**, *227*, 1199–1203. [CrossRef]

6. Köppel, R.; Ruf, J.; Rentsch, J. Multiplex real-time PCR for the detection and quantification of DNA from beef, pork, horse and sheep. *Eur. Food Res. Technol.* **2011**, *232*, 151–155. [CrossRef]

7. Kaltenbrunner, M.; Hochegger, R.; Cichna-Markl, M. Tetraplex real-time PCR assay for the simultaneous identification and quantification of roe deer, red deer, fallow deer and sika deer for deer meat authentication. *Food Chem.* **2018**, *269*, 486–494. [CrossRef] [PubMed]

8. Staats, M.; Arulandhu, A.J.; Gravendeel, B.; Holst-Jensen, A.; Scholtens, I.; Peelen, T.; Prins, T.W.; Kok, E. Advances in DNA metabarcoding for food and wildlife forensic species identification. *Anal. Bioanal. Chem.* **2016**, *408*, 4615–4630. [CrossRef] [PubMed]

9. Nehal, N.; Choudhary, B.; Nagpure, A.; Gupta, R.K. DNA barcoding: A modern age tool for detection of adulteration in food. *Crit. Rev. Biotechnol.* **2021**, *41*, 767–791. [CrossRef] [PubMed]

10. Galimberti, A.; Casiraghi, M.; Bruni, I.; Guzzetti, L.; Cortis, P.; Berterame, N.M.; Labra, M. From DNA barcoding to personalized nutrition: The evolution of food traceability. *Curr. Opin. Food Sci.* **2019**, *28*, 41–48. [CrossRef]

11. Dobrovolny, S.; Blaschitz, M.; Weinmaier, T.; Pechatschek, J.; Cichna-Markl, M.; Indra, A.; Hufnagl, P.; Hochegger, R. Development of a DNA metabarcoding method for the identification of fifteen mammalian and six poultry species in food. *Food Chem.* **2019**, *272*, 354–361. [CrossRef] [PubMed]

12. Preckel, L.; Brünen-Nieweler, C.; Denay, G.; Petersen, H.; Cichna-Markl, M.; Dobrovolny, S.; Hochegger, R. Identification of mammalian and poultry species in food and pet food samples using 16s rdna metabarcoding. *Foods* **2021**, *10*, 2875. [CrossRef] [PubMed]

13. ISO 21571 I; Foodstuffs—Methods of Analysis for the Detection of Genetically Modified Organisms and Derived Products—Nucleic Acid Extraction; ISO: Geneva, Switzerland, 2005.

14. Meyer, R.; Höfelein, C.; Lüthy, J.; Candrian, U. Polymerase chain reaction-restriction fragment length polymorphism analysis: A simple method for species identification in food. *J. AOAC Int.* **1995**, *78*, 1542–1551. [CrossRef] [PubMed]

15. Palumbi, S.R.; Martin, A.; Romano, S.; McMillan, W.O.; Stice, L.; Grabowski, G. *The Simple Fool's Guide to PCR*; Department of Zoology and Kewalo Marine Laboratory, University of Hawaii: Honolulu, HI, USA, 2002.

16. Martin, M. Cutadapt removes adapter sequences from high-throughput sequencing reads. *EMBnet J.* **2011**, *17*, 10. [CrossRef]

17. Aronesty, E. Comparison of sequencing utility programs. *Open Bioinform. J.* **2013**, *7*, 1–8. [CrossRef]

18. Edgar, R.C. Search and clustering orders of magnitude faster than BLAST. *Bioinformatics* **2010**, *26*, 2460–2461. [CrossRef] [PubMed]

19. Quodata. PROLab Plus Software for Planning, Organizing, Performing and Analyzing Interlaboratory Studies. Available online: www.quadata.de (accessed on 10 December 2021).

20. *Official Collection of Test Methods A. Statistik—Planung und Statistische Auswertung von Ringversuchen zur Methodenvalidierung*; Federal Office of Consumer Protection and Food Safety (BVL), Beuth Verlag: Berlin, Germany, 2006.

21. MacConaill, L.E.; Burns, R.T.; Nag, A.; Coleman, H.A.; Slevin, M.K.; Giorda, K.; Light, M.; Lai, K.; Jarosz, M.; McNeill, M.S.; et al. Unique, dual-indexed sequencing adapters with UMIs effectively eliminate index cross-talk and significantly improve sensitivity of massively parallel sequencing. *BMC Genom.* **2018**, *19*, 30. [CrossRef] [PubMed]
Article

# Development of a Duck Genomic Reference Material by Digital PCR Platforms for the Detection of Meat Adulteration

Xiaoyun Chen<fnref n="1" />, Yi Ji<fnref n="1" />, Kai Li<fnref n="2" />, Xiaofu Wang<fnref n="1" />, Cheng Peng<fnref n="1" />, Xiaoli Xu<fnref n="1" />, Xinwu Pei<fnref n="2" />, Junfeng Xu<fnref n="1" />,*
and Liang Li<fnref n="2" />,†

**Abstract:** Low-cost meat, such as duck, is frequently used to adulterate more expensive foods like lamb or beef in many countries. However, the lack of DNA-based reference materials has limited the quality control and detection of adulterants. Here, we report the development and validation of duck genomic DNA certified reference materials (CRMs) through the detection of the duck *interleukin 2* (IL2) gene by digital PCR (dPCR) for the identification of duck meat in food products. The certified value of IL2 in CRMs was 5.78 ± 0.51 × 10³ copies/µL with extended uncertainty (coverage factor *k* = 2) based on IL2 quantification by eight independent collaborating laboratories. Quantification of the mitochondrial gene *cytb* revealed a concentration of 2.0 × 10⁶ copies/µL, as an information value. The CRMs were also used to determine the limit of detection (LOD) for six commercial testing kits, which confirmed that these kits meet or exceed their claimed sensitivity and are reliable for duck detection.

**Keywords:** certified reference material; digital PCR; duck *interleukin 2*; mitochondrial gene; meat adulteration

**Citation:** Chen, X.; Ji, Y.; Li, K.; Wang, X.; Peng, C.; Xu, X.; Pei, X.; Xu, J.; Li, L. Development of a Duck Genomic Reference Material by Digital PCR Platforms for the Detection of Meat Adulteration. *Foods* **2021**, *10*, 1890. https://doi.org/10.3390/foods10081890

Academic Editors:
Mohammed Gagaoua and
Brigitte Picard

Received: 13 July 2021

Accepted: 13 August 2021

Published: 15 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_66.png" /></figure>

Copyright: © 2021 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article
distributed under the terms and
conditions of the Creative Commons
Attribution (CC BY) license (https://
creativecommons.org/licenses/by/
4.0/).

## 1. Introduction

Meat adulteration has become a major global issue. Fraud by substitution or adulteration with inexpensive raw materials poses an attractive shortcut for unscrupulous food producers, although it is fraught with potential public health risks [1]. In order to definitely show that a product has been adulterated or fraudulently labeled, the product composition must be first determined for comparison of its authenticity with the description provided on its label. This process requires quantitative analysis of characteristic compounds or analytes specific to the ingredient in question, or other evidence that it is present in concentrations at or above levels required by regulatory agencies [2]. Therefore, the development of an accurate and reliable method for determining the content of specific meat components in food products has major economic implications and far-reaching social significance.

Traditional methods for species identification have historically relied on anatomy, histology, sensory judgment, chemistry, electrophoresis, chromatography, and immunology [3]. However, these approaches are each accompanied by limitations in their accuracy and/or sensitivity. To address this issue, recent studies have developed polymerase chain reaction (PCR)-based assays to accommodate the discriminatory capability necessary for identification of specific ingredients or contaminants. The most widely used methods are DNA-based screens [4], including gel-based PCR [5], real-time qPCR [6], multiplex PCR [7], digital PCR [8], isothermal nucleic acid amplification [9–11], and other PCR techniques [12,13]. The recent and extremely rapid expansion in detection methods and the

<footnote n="1">¹ State Key Laboratory for Managing Biotic and Chemical Threats to the Quality and Safety of Agro-Products, Zhejiang Academy of Agricultural Sciences, Hangzhou 310021, China; xiaoyunchen_2016@163.com (X.C.); jymemory12138@163.com (Y.J.); yywxf1981@163.com (X.W.); pc_phm@163.com (C.P.); xuxiaoli@zju.edu.cn (X.X.)</footnote>

<footnote n="2">² Biotechnology Research Institute, Chinese Academy of Agricultural Sciences, Beijing 100081, China; likaij@163.com (K.L.); peixinwu@caas.cn (X.P.)</footnote>

<footnote>* Correspondence: njjfxu@163.com (J.X.); liliang@caas.cn (L.L.)</footnote>
standardization of detection methods for animal-based food products has inadvertently circumvented some steps that are essential to ensure the rigor and quality of data, which is compounded by a lack of reference materials. These factors together can incur a bottle-neck in data quality assurance. To rectify this issue in data quality assurance in food production, many reference materials are urgently needed for the ongoing evaluation of the methods used for quality control.

A reference material (RM) is sufficiently homogeneous and stable with respect to one or more specified properties that have been established to be fit for its intended use in a measurement process. Certified reference materials (CRMs) are reference materials (RMs) characterized using a metrologically valid procedure for one or more specified properties, accompanied by an RM certificate [14], which can be used for calibration of a measurement system, assessment of a measurement procedure, for assigning values to other materials, and quality control. Several DNA-based CRMs have been developed, such as for genetically modified organisms [15–17], cancer diagnosis [18–20], foodborne pathogens [21], and forensic science [22,23].

Due to the increasingly high frequency of adulteration of meat products with inex-
pensive duck, this work aimed to develop a novel DNA reference material for duck meat
through the targeted detection and quantification of the duck *IL2* gene. Here, we describe
the preparation, homogeneity, and short- and long-term stability of the gDNA materials.
The mean *IL2* copy number values of these CRMs for duck were validated through digital
PCR (dPCR) by eight independent collaborating diagnostic laboratories, and the materials
were certified by China's State Administration for Market Regulation. In addition, the
reference materials were used to evaluate the limit of detection for six commercial testing
kits. These CRMs are intended for qualitative and quantitative screening for duck meat
in meat products through detection of the *IL2* gene. These CRMs can also be used for
validation of other quantitative DNA-based methods and laboratory quality control.

## 2. Materials and Methods

2.1. *Extraction and Evaluation of Duck Genomic DNA*

Duck leg meat was provided by the Institute of Animal Husbandry and Veteri-
nary Science, Zhejiang Academy of Agricultural Sciences. Duck leg meat was used,
and a sterile scalpel was used to remove the skin and cut it into small pieces, after
which it was blended in Philips Mixture HR2027 (Hongkong, China), and then trans-
ferred to a Nalgene® 3118-0050 Oak Ridge Centrifuge Tube. Genomic DNA extrac-
tion was performed using a Simgen Animal Tissue DNA Midi Kit (Hangzhou, China).
The integrity of the DNA critically affects the success of the characterization, which
was analyzed by gel electrophoresis using a DYY-6C gel electrophoresis system (Liuyi
Biotechnology, Beijing, China). The quality and purity of gDNA was evaluated at 230 nm,
260 nm, and 280 nm by the UV absorbance method (Nanodrop™ 2000, Thermo Fisher,
Wilmington, DE, USA). The extracted gDNA was quantified using the Quant-iT™ dsDNA
PicoGreen® Kits (Invitrogen, Shanghai, China) according to the manufacturer's instruc-
tions using a lambda DNA standard solution.

## 2.2. Digital PCR Assay

2.2.1. Certified Value

The target gene is a single-copy nuclear gene in the genome. The primers and probe (Figure 1) used for specific quantification were as described previously [24]. The primers and probe were synthesized by Sangon (Shanghai, China). The amplicon size was 212 bp. The primer and probe sequences are shown as follows:

5'-GGAGCACCTCTATCAGAGAAAGACA-3'

5'-GTGTGTAGAGCTCAAGATCAATCCC-3'

5'-FAM-TGGGAACAAGCATGAATGTAAGTGGATGGT-BHQ1-3'
<figure><img src="image_67.png" /><figcaption>**Figure 1.** Sequence of the *interleukin 2* precursor gene (AY821656.1). The primers (red) and probe (green) for the digital PCR are underlined. The amplicon size was 212 bp.</figcaption></figure>

Optimization of the PCR condition is critical for developing a PCR-based method for the specific detection of duck. This study fully optimized the key parameters of real-time qPCR and digital PCR, including specificity, oligonucleotide concentration, annealing temperature, dynamic range, the limit of detection, and the limit of quantification (data not shown). The optimized PCR reaction contains 10 µL of 2X ddPCR Supermix for probes (no dUTP) (Bio-Rad, Hercules, CA, USA), 1 µL of each primer (0.5 µM), 0.5 µL of probe (0.25 µM), 2 µL of DNA template, and 5.5 µL of ddH₂O. After that, the 20 µL of the reaction mixture was then loaded on eight-channel disposable droplet generator cartridges (Bio-Rad). Droplets were generated with 70 µL of droplet generation oil (Bio-Rad) in the droplet generator of the QX200 system (Bio-Rad). The generated droplets were transferred to a 96-well PCR plate (Bio-Rad). The amplification was carried out at a uniform ramp rate of 2.5 °C/s for 95 °C for 10 min, 45 cycles of 95 °C for 15 s followed by 60 °C for 1 min, and a final enzyme deactivation at 98 °C for 10 min. Fluorescent signals from amplified droplets were captured individually in the QX200 droplet reader (Bio-Rad) and analyzed with QuantaSoft 1.6.6.0320 software. The target concentrations were reported as the number of copies per µL of PCR reaction after correction with the Poisson distribution.

### 2.2.2. Information Value

A mitochondrial *cytb* gene TaqMan MGB probe assay [25] was used. For absolute quantification, the primers and probe were fully performed in a QX200 ddPCR system. Primers and probes were purified with high-performance liquid chromatography. The preparation of the reaction mixture and optimized PCR thermal procedure is listed in the Supplementary Materials. The primers and probes are as follows:

5'-GGCACACAAATCCTCACAG-3'

5'-TGTGTTGGCTACTGAGGAGAAA-3'

5'-FAM- CCTACTGGCTATGCACTACACCGCAGAC-BHQ1-3'
## 2.3. Homogeneity Testing

Most CRMs are prepared as batches of bottles or vials. It is important that all units are the same within the stated uncertainty for each property value. ISO 17034 accordingly requires the assessment of the homogeneity of a certified reference material. Homogeneity includes within-unit homogeneity and between-unit homogeneity. It is always necessary to assess the between-unit homogeneity, which was assessed to ensure the equivalence between different units following the Section 7.5.2 of the ISO Guide 35:2017. Furthermore, the within-unit homogeneity is directly reflected in the minimum size of the subsample. The homogeneity of the *IL2* duck gene was evaluated by randomly selecting 15 units from the CRM candidates. Three subsamples from different positions of each unit were taken and measured by dPCR. Measurement results were used for the assessment of homogeneity. Statistical analysis of the data was performed for a 95% confidence level. The values of the copy number concentration for this *IL2* gene of the duck samples were analyzed using the F test and compared with the significance of the calculated F value of the copy number concentration and the critical $F_{0.05(14, 30)}$ value of 2.04.

An experimental assessment must consist of a determination of the minimum sample intake. It is always essential to ensure that the sample intake is sufficient. The between-unit study does not provide such assurance, so the experimental within-unit homogeneity study should be carried out.

## 2.4. Stability Monitoring

Stability is an important parameter for the reference material. The stability of all CRMs should be assessed. Two types of stability are relevant in the production of reference material: the long-term stability and transportation stability (short-term stability). This CRM is the DNA solution, so it needs to be repeatedly frozen and thawed, and thus the freeze-thaw stability is a concern. Stability testing was carried out to evaluate the influence of different storage temperatures and times on the CRMs using digital PCR (dPCR). The uncertainty of stability needs to include the uncertainty of long-term stability and short-term stability.

Transportation stability (short-term stability) is a property of the material referring to stability under expected transport conditions. Extreme high temperatures cannot be ruled out during transportation, so the CRM units were stored at 4 °C, 25 °C, 37 °C, and 60 °C for 1, 3, 5, 7, and 14 days. The samples were stored at -70 °C after sampling. Three tubes were randomly selected for each storage temperature and each tube was sampled 3 times. The digital PCR tested the stability of *IL2* gene of the duck samples. The *t*-test was performed on the chosen dates to evaluate the short-term stability of the CRMs. The *t*-test results showed no significant slope (β1) for the CRMs at 4 °C, 25 °C, and 37 °C at the 95% confidence level during the two weeks.

Long-term stability studies are conducted to assess stability under storage conditions specified for the lifetime of the product. The long-term stability study was extended to 6 months, the *IL2* gene of the duck CRMs was evaluated by analyzing 3 tubes stored at -20 °C for 1, 2, 4, and 6 months. The samples were stored at -70 °C after sampling. Three tubes were randomly selected for each storage temperature and each tube was sampled 3 times by the digital PCR. Because the *IL2* gene of the duck CRMs is the DNA solution, their freeze-thaw stability needs to be implemented. This batch of CRMs was 100 µL per tube and the minimum sample was 2 µL. In the experiment, 3 tubes of the *IL2* gene of the duck CRMs were taken, repeated freezing-thawing was performed 10 times, and samples were taken 3 times per tube by the digital PCR.

## 2.5. Collaborative Characterization

The property value of this batch of CRMs was determined by eight laboratories (Supplementary Table S1) using dPCR. The eight laboratories were engaged in DNA measurement for a long time and have a certain technical authority.
The gDNA samples of the *IL2* gene of ducks and the related ddPCR reagents were mailed to each participating laboratory in a closed box filled with dry ice. Each laboratory received two gDNA samples (label: duck1, duck2), each with a volume of 100 µL. In order to simplify the sampling procedure and reduce the deviation, primers/probes were mixed according to the proportion before sending them to each laboratory. The participating laboratories were requested to measure the copy number concentration by the operation protocol of dPCR experiments. Each sample was repeatedly measured four times, and a total of eight subsamples were required simultaneously on the same PCR plate for each participant. The eight participants need to export data and original files to the organization lab. The returned data were analyzed according to the requirements of ISO Guide 35:2017.

## 2.6. Evaluation of Limit of Detection of Commercial Assay Kits

In this experiment, six commercial duck-derived detection qPCR kits were purchased. There were four mitochondrial gene targets and two nuclear gene kits. Company and product information are summarized in Supplementary Table S2. The CRM was serially diluted and applied to limit of detection (LOD) probit regression analysis for six diagnostic assays [26]. Each of the concentration levels were tested with multiple replicates per concentration according to the manufacturer's instructions with blank controls. Probit regression analysis of 95% hit rates was performed with SPSS 16.0 software (SPSS Inc., Chicago, IL, USA).

# 3. Results

## 3.1. gDNA Extracted from Commercially Obtained Duck Leg Meat

Prior to testing for adulteration, we first sought to confirm the quality of gDNA extracted from leg meat samples of duck obtained from the Institute of Animal Husbandry and Veterinary Science, Zhejiang Academy of Agricultural Sciences. To this end, we first analyzed the extracted gDNA by 1% agarose gel electrophoresis and found neither a visible smear nor an RNA band, which indicated that the gDNA was intact (i.e., not degraded) and that RNA had been removed during the isolation process. In addition, the quantity and quality of each gDNA extraction were confirmed in six technical replicates using a Nanodrop 2000 spectrometer (Thermo Scientific, Wilmington, DE, USA). The A260/A280 ratio averaged 1.98, while the average A260/A230 ratio was 2.28, indicating high purity. The quantification of gDNA concentration by the Qubit 2.0 Fluorometer (Invitrogen) showed that the samples contained approximately 100 ng/µL. DNA solutions were then diluted to concentrations of 6.5 ng/µL with 0.1 × Tris-EDTA (TE) buffer and stored in a -20 °C freezer.

## 3.2. Bottling and Storage/Homogeneity Testing and Minimum Sample Intake

In order to establish protocols for future sample collection and preparation, we next established standardized steps for the proper storage and handling of DNA extractions following methods described before [17]. The gDNA solutions were sterilized using a 0.2 µm nylon filter and 100 µL aliquots were placed in low static, sterilized polypropylene microcentrifuge tubes, with approximately 5 × 10³ copies of duck gDNA per tube. Cardboard freezer boxes accommodating 100 samples each were used for storage, and random samples were taken from each box for testing homogeneity as well as short- and long-term stability. All samples were stored in 4 °C refrigerators in the dark.

In order to establish the protocol for evaluating the duck certified reference material (CRM), we examined the homogeneity of the *IL2* gene in our gDNA extractions using a random stratified sampling method [27] to select 15 units from among 500 CRM candidate samples. Three subsamples were taken from different vial positions for each sample and analyzed by digital PCR (dPCR) to determine whether this technique was sufficiently sensitive to identify differences in *IL2* uniformity within and between samples (for raw data see Supplementary Table S3). Significant differences in the *IL2* copy number among within-vial subsamples and between different samples were then determined using an
F test to compare the Fcalculated with the Fcritical value (2.04) for a 95% confidence level. The results showed that the Fcalculated (1.53) value was lower than the Fcritical (2.04) value (Table 1), which thus indicated high homogeneity among the CRM gDNA samples in this batch.

<table><caption>**Table 1.** Results of the homogeneity analysis.</caption><thead><tr><th>Parameter</th><th>Copy Number/Unit</th></tr></thead><tbody><tr><td>Mean</td><td>5.76 × 10<sup>3</sup></td></tr><tr><td>Q<sub>1</sub></td><td>3.59 × 10<sup>5</sup></td></tr><tr><td>V<sub>1</sub></td><td>14</td></tr><tr><td>S<sub>1</sub><sup>2</sup></td><td>2.57 × 10<sup>4</sup></td></tr><tr><td>Q<sub>2</sub></td><td>5.02 × 10<sup>4</sup></td></tr><tr><td>V<sub>2</sub></td><td>30</td></tr><tr><td>S<sub>2</sub><sup>2</sup></td><td>1.67 × 10<sup>4</sup></td></tr><tr><td>F</td><td>1.53</td></tr><tr><td>F<sub>0.05(14, 30)</sub></td><td>2.04</td></tr><tr><td>Conclusion</td><td>F &lt; F<sub>0.05(14, 30)</sub></td></tr><tr><td>μ<sub>bbrel</sub></td><td>0.0095</td></tr></tbody></table>

### 3.3. Stability Monitoring and Freeze-Thaw Cycles Testing

The stability of the properties of interest (in this case, the duck IL2 copy number) represents an essential feature of any given reference material. It is well-known that long-term stability is related to storage conditions, while short-term stability is related to external factors during sample transportation. Moreover, stability studies can be categorized as either classical stability or synchronous stability studies. In this work, we incorporated data from other international studies [16,19,28] characterizing DNA standard materials in a classical stability assessment of duck IL2 gene stability at −70 °C. Under these storage conditions, we found that the IL2 copy number of the candidate duck CRMs did not change. It warrants mention here that the copy number stability at other storage temperatures requires a synchronous stability study for comparison with storage at −70 °C.

Transportation stability (short-term stability) is a property of the material that describes its stability under expected transport conditions in order to guide proper handling methods prior to evaluation. Since extremely high temperatures cannot be ruled out during transportation, gDNA CRMs were stored at 4 °C, 25 °C, and 60 °C for 1, 3, 5, 7, and 14 days. Three tubes were randomly selected for storage at each temperature, and each tube was tested in three technical replicates by digital PCR. The results of the t-test showed no significant differences in the slope (β1) between the CRMs stored for two weeks at 4 °C or 25 °C at a 95% confidence level (Supplementary Table S4 and Figure 2). However, at 60 °C, the IL2 copy number significantly decreased after a single day of storage. The relative uncertainty in the IL2 copy number due to instability at 25 °C storage was calculated to be 0.035. Collectively, these results of short-term CRM stability analysis indicated that duck gDNA is stable under room temperature (25 °C) storage and transport for up to 14 days, though we advise cold chain transportation to minimize the likelihood of degradation.
<figure><img src="image_68.png" /><figcaption>**Figure 2.** Trends and concentration of stability analysis for the candidate materials: (A) long-term stability at −20 °C and (B–D) short-term stability at 4 °C, 25 °C, and 60 °C. The error bars represent the standard deviation of within-bottle triplicate subsamples.</figcaption></figure>

In order to determine the long-term stability of the CRMs, three randomly selected tubes were stored at −20 °C and tested for the *IL2* copy number by dPCR at 1, 2, 4, and 6 months time points. The results of this analysis showed no significant differences in the *IL2* copy number among any of the time points throughout the 6-month experiment at −20 °C (*p* < 0.05) (Supplementary Table S5 and Figure 2). We then estimated a 0.020 relative uncertainty of the copy number instability caused by long-term storage at −20 °C. Collectively, these results showed that gDNA CRMs could be stably stored at −20 °C for at least 6 months.

Finally, we investigated whether freeze-thaw cycles could adversely affect the *IL2* copy number. Since this batch of CRM samples each contained 100 µL per tube and our above results showed that 2 µL of template was appropriate for analysis, we then randomly selected three tubes of duck CRM gDNA and performed 10 repeated freezing-thaw cycles. Subsequently, dPCR was performed using each sample in three technical replicates, which showed no significant difference in the duck *IL2* gene copy number (*p* > 0.05) among samples. These results indicated that duck gDNA CRMs are stable and robust against damage due to repeated freezing and thawing (data not shown).

### 3.4. Characterization of Certified and Information Value

To validate our results through independent third parties, eight different laboratories each determined the *IL2* copy number for two randomly selected duck CRM samples using four technical replicates per sample in dPCR. Thus, each participating external lab returned eight dPCR results (Supplementary Table S6, Figure 3), and the copy numbers were estimated and analyzed according to ISO GUIDE 35:2017. The statistical analysis indicated that the datasets followed a normal distribution, and none of the datasets contained points outside of the uncertainty range for a 95% confidence level. Furthermore, no significant differences were found in the mean *IL2* copy numbers from each lab, and thus, the certified mean *IL2* copy number in the duck CRMs was estimated to be 5.78 × 10³ copies/µL.
<figure><img src="image_69.png" /><figcaption>Figure 3. Results for dPCR quantification of the IL2 gene copy number by eight independent participating laboratories. Solid line, average of participant results (♦); dotted lines, expanded uncertainty of the certified values of the copy number concentration with a 95% level of confidence.</figcaption></figure>

In addition to IL2, mitochondrial genes can also serve as quantifiable targets for the amplification of duck gDNA in meat-based food products. Our laboratory optimized a digital PCR method for detection of the mitochondrial gene cytb, which revealed a concentration of 2.0 × 10⁶ copies/µL in the CRM samples.

### 3.5. Statistical Estimation of Uncertainty

The uncertainty of the CRMs for the IL2 gene in these duck CRMs was estimated from the contributions according to ISO Guide 35 and consisted of uncertainty components from characterization (u$_{char}$), potential between-unit heterogeneity (u$_{bb}$), and potential instability during long-term (u$_{its}$) and short-term storage (u$_{sts}$).

The relative standard uncertainty of characterization (u$_{char,rel}$) is estimated by Equation (1). The result is 0.012.

$$u_{\text{char,rel}} = \sqrt{u_{\text{A,rel}}^2 + u_{\text{B,rel}}^2} \qquad (1)$$

where u$_{A,rel}$ is type A uncertainty (random error) and u$_{B,rel}$ is type 2 uncertainty (systematic error).

The relative standard uncertainty of potential between-unit heterogeneity (u$_{bb,rel}$) is 0.0095, estimated as described in Table 1.

The relative standard uncertainty of the potential degradation during transport (u$_{sts,rel}$), and long-term storage (u$_{its,rel}$) is estimated by Equations (2) and (3). The results are 0.035 and 0.020, respectively.

$$u_{\text{sts}} = s(\beta_1) \times X/\bar{x} \qquad (2)$$

where s(β₁) is the standard deviation of all results of the transport time stability study (Supplementary Table S4), X is the chosen transport time (14 days at 25 °C), and x̄ is the mean value.

$$u_{\text{its}} = s(\beta_1) \times X/\bar{x} \qquad (3)$$

where s(β₁) is the standard deviation of all the results of the long-time stability study (Table S5), X is the chosen transport time (six months at −20 °C), and x̄ is the mean value.
The expanded relative uncertainty (U<sub>CRM,rel</sub>) for the IL2 gene in these duck CRMs was estimated by Equation (4), and the expanded uncertainty of the copy number concentration (U<sub>CRM</sub>) was estimated to be 0.51 × 103 (coverage factor k = 2, approximate 95% confidence interval).

$$U_{\text{CRM,rel}} = k \times \sqrt{u_{\text{char,rel}}^2 + u_{\text{bb,rel}}^2 + u_{\text{its,rel}}^2 + u_{\text{sts,rel}}^2} \quad (4)$$

where u<sub>char</sub> is the measurement of the certified value, u<sub>bb</sub> is the potential between-unit heterogeneity, u<sub>its</sub> and u<sub>sts</sub> is the uncertainties in the long-term and short-term stability.

### 3.6. Results of Assay Kit Evaluations

We then investigated whether commercially available kits were sufficiently sensitive to be used for extraction of duck gDNA for comparison with CRMs characterized here. For this purpose, we used the externally validated CRM samples and accompanying data to evaluate the LODs of six diagnostic assays (Figure 4). The lot numbers and LODs claimed by each commercial diagnostic assay are listed in Supplementary Table S2. Then, using CRM samples, we tested each of these six kits and compared the results of the LOD determination using a probit regression analysis of 95% hit rates with the LOD in the manufacturer’s instructions for each kit. The results showed that all of the assays met or exceeded their claimed sensitivities and are reliable for detection of IL2 in gDNA extracted from duck meat. Moreover, the results also show that these duck gDNA CRMs are suitable for the assessment of commercial kits and allow for comparable LOD studies of various detection methodologies.

<figure><img src="image_70.png" /><figcaption>**Figure 4.** Evaluation of LODs of commercial assays. (A) Comparison of LODs claimed by the kit manufacturer using six duck DNA detection kits with the LODs determined with the CRMs. (B–G) Probit regression analysis of the six assays for the detection of duck IL2 (SPSS 16.0). The probit versus IL2 concentration ((B–E) mitochondrial gene, (F,G) nuclear gene) was obtained from 21 replicates of serial dilutions and an additional 10 replicates of a blank sample.</figcaption></figure>
4. Discussion

In recent years, DNA detection methods have been an essential technique for molecu-
lar diagnosis. RM and CRM are the quality assurance of testing data, which can evaluate the
measurement method and monitor the measurement process. DNA CRMs originated from
the detection of genetically modified ingredients because quantitative testing is needed
in this field [15]. Subsequently, a large number of molecular diagnostic standards have
been developed in the field of medical testing. The international metrology field also pays
great attention to DNA measurement, especially traceability to SI units [28]. As an absolute
quantitative method, the digital PCR method is currently a potential primary reference
measurement method for DNA target measurement. As the carrier of the measurement
value, the reference material plays a vital role in the measurement traceability and the
guarantee of the data quality of the measurement results.

In this study, we chose the nuclear gene *IL2* as the target gene. Compared with
mitochondrial genes, nuclear genes are single-copy genes, genetically stable, and there is
no difference between different tissues [6], so it is easy to perform absolute quantification.
We fully optimized the *IL2* gene method on qPCR and dPCR platforms. Although the
results of this study characterizing a CRM for duck meat are primarily based on dPCR,
standard laboratory testing typically relies on qPCR-based analysis. Therefore, it was
also necessary to investigate whether the gDNA samples contained inhibitors that could
interfere with the qPCR reaction. We first performed six serial dilutions of each gDNA
sample, ranging from 100 ng/µL to 0.02 ng/µL, and generated a standard curve of Ct
values by qPCR. The linear regression equation between Ct (y) and log10 starting copy
number (X) was Y = -3.415X + 42.632. The slope of the curve was -3.415, indicating a
96.3% amplification efficiency and R² = 0.998. The linear regression analysis thus confirmed
that the duck gDNA samples were suitably pure for qPCR and met the performance criteria
defined by the European Network of GMO Laboratories (ENGL). We next determined
that the appropriate linear range for the copy number for subsequent dPCR analysis was
between 20 and 14,000 copies per reaction to obtain a good correlation (i.e., R² = 1) with the
expected value. Digital PCR commonly uses 20 µL reaction volumes containing 2 µL of the
DNA template. The results of the dPCR showed that the *IL2* copy number in these duck
CRM samples averaged 5780 copies per microliter, within the linear range, which suggested
that 2 µL of the gDNA template was appropriate for analysis. We then determined the
standard uncertainty of heterogeneity among technical replicates of each sample to account
for variation between samples. The relative uncertainty for the *IL2* copy number was 0.0095,
which corresponded with a 95% confidence level. Taken together, these results show that
the samples were sufficiently homogenous for accurate evaluation as a standard CRM.

The CRM has two values to meet the detection needs of different targets. The certified
value is expressed as (5.78 ± 0.51) × 10³ copies/µL with extended uncertainty (a coverage
factor *k* = 2) based on the quantification of *IL2* with a collaborative characterization of eight
participants. An information value is 2.0 × 10⁶ copies/µL based on the quantification of a
mitochondrial gene. Our laboratory determined this value as an informative supplementary
value for the quantification of duck meat, which may be of value to users of the CRMs.
However, the replication data are currently insufficient to assess the uncertainty associated
with this value rigorously.

Lots of novel detection methods of meat adulteration have been developed [29,30],
while the literature on reference materials is rarely reported. Some commercial companies
have developed RM, but not CRM. Therefore, it is urgent to develop animal DNA standard
materials according to the ISO standard system. The identification of meat adulteration
is an important part of the field of food safety. DNA-based testing is almost qualitative,
and it needs to be further developed into quantitative testing. At the same time, CRM
should be developed, international proficiency testing should be organized, and assay kit
evaluations should be carried out to ensure accurate and consistent measurement results
in the industry.
**5. Conclusions**

Here, we developed gDNA CRMs (GBW(E) 091060) certified by China's State Administration for Market Regulation for the analysis of duck adulteration in meat products. Eight independent diagnostic laboratories contributed to the validation and certification of these duck gDNA CRMs by digital PCR, confirming a mean of 5.78 ± 0.51 × 10³ copies/µL for the IL2 gene. Homogeneity and stability testing demonstrated that the CRMs were homogenous and stable for at least 6 months at -20 °C storage and for 14 days in cold chain delivery conditions. This batch of duck gDNA CRMs can serve as an essential tool for method validation and proficiency testing in the analysis of duck content in meat products. This study also provides a technical basis for development of other animal-derived reference materials.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10.3390/foods10081890/s1, Table S1: Laboratories participating in the characterization of the IL2 gene from duck genomic DNA. Table S2: List of commercial assays assessed by the IL2 certified reference standard. Table S3: Raw data of the homogeneity analysis. Table S4: Results of the short-term stability study. Table S5: Results of the long-term stability study. Table S6: Results of independent laboratory validation.

**Author Contributions:** Conceptualization, L.L. and J.X.; methodology, X.C. and X.X.; software, X.C.; validation, Y.J. and K.L.; formal analysis, Y.J. and K.L.; investigation, X.W. and C.P.; resources, L.L. and J.X.; writing—original draft preparation, X.C.; writing—review and editing, L.L. and J.X.; visualization, X.W.; supervision, C.P. and X.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Key Research and Development Program of Zhejiang (2021C02059); the Agricultural Science and Technology Innovation Program of CAAS (CAAS-ZDRW202009).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data are available upon request.

**Acknowledgments:** The authors gratefully acknowledge the financial support from the Key Research and Development Program of Zhejiang (2021C02059) and the Agricultural Science and Technology Innovation Program of CAAS (CAAS-ZDRW202009).

**Conflicts of Interest:** The authors declare no conflict of interest.

**References**

1. Köppel, R.; Daniels, M.; Felderer, N.; Brünen-Nieweler, C. Multiplex real-time PCR for the detection and quantification of DNA from duck, goose, chicken, turkey and pork. *Eur. Food Res. Technol.* **2013**, *236*, 1093–1098. [CrossRef]

2. Burns, M.; Wiseman, G.; Knight, A.; Bramley, P.; Foster, L.; Rollinson, S.; Damant, A.; Primrose, S. Measurement issues associated with quantitative molecular biology analysis of complex food matrices for the detection of food fraud. *Analyst* **2015**, *141*, 45–61. [CrossRef]

3. Kumar, A.; Kumar, R.R.; Sharma, B.D.; Gokulakrishnan, P.; Mendiratta, S.K.; Sharma, D. Identification of Species Origin of Meat and Meat Products on the DNA Basis: A Review. *Crit. Rev. Food Sci. Nutr.* **2013**, *55*, 1340–1351. [CrossRef] [PubMed]

4. Böhme, K.; Calo-Mata, P.; Barros-Velázquez, J.; Ortea, I. Review of Recent DNA-Based Methods for Main Food-Authentication Topics. *J. Agric. Food Chem.* **2019**, *67*, 3854–3864. [CrossRef] [PubMed]

5. Kitpipit, T.; Sittichan, K.; Thanakiatkrai, P. Are these food products fraudulent? Rapid and novel triplex-direct PCR assay for meat identification. *Forensic Sci. Int. Genet. Suppl. Ser.* **2013**, *4*, e33–e34. [CrossRef]

6. Li, T.; Wang, J.; Wang, Z.; Qiao, L.; Liu, R.; Li, S.; Chen, A. Quantitative determination of mutton adulteration with single-copy nuclear genes by real-time PCR. *Food Chem.* **2021**, *344*, 128622. [CrossRef] [PubMed]

7. Xu, R.; Wei, S.; Zhou, G.; Ren, J.; Liu, Z.; Tang, S.; Cheung, P.C.; Wu, X. Multiplex TaqMan locked nucleic acid real-time PCR for the differential identification of various meat and meat products. *Meat Sci.* **2018**, *137*, 41–46. [CrossRef]

8. Floren, C.; Wiedemann, I.; Brenig, B.; Schütz, E.; Beck, J. Species identification and quantification in meat and meat products using droplet digital PCR (ddPCR). *Food Chem.* **2015**, *173*, 1054–1058. [CrossRef]
9. Lin, L.; Zheng, Y.; Huang, H.; Zhuang, F.; Chen, H.; Zha, G.; Yang, P.; Wang, Z.; Kong, M.; Wei, H.; et al. A visual method to detect meat adulteration by recombinase polymerase amplification combined with lateral flow dipstick. *Food Chem.* **2021**, *354*, 129526. [CrossRef]

10. Cao, Y.; Zheng, K.; Jiang, J.; Wu, J.; Shi, F.; Song, X.; Jiang, Y. A novel method to detect meat adulteration by recombinase polymerase amplification and SYBR green I. *Food Chem.* **2018**, *266*, 73-78. [CrossRef] [PubMed]

11. Wang, J.; Wan, Y.; Chen, G.; Liang, H.; Ding, S.; Shang, K.; Li, M.; Dong, J.; Du, F.; Cui, X.; et al. Colorimetric Detection of Horse Meat Based on Loop-Mediated Isothermal Amplification (LAMP). *Food Anal. Methods* **2019**, *12*, 2535-2541. [CrossRef]

12. Hossain, M.A.M.; Ali, E.; Hamid, S.B.A.; Asing; Mustafa, S.; Desa, M.N.M.; Zaidul, I.S.M. Double Gene Targeting Multiplex Polymerase Chain Reaction-Restriction Fragment Length Polymorphism Assay Discriminates Beef, Buffalo, and Pork Substitution in Frankfurter Products. *J. Agric. Food Chem.* **2016**, *64*, 6343-6354. [CrossRef]

13. López-Oceja, A.; Nuñez, C.; Baeta, M.; Gamarra, D.; de Pancorbo, M. Species identification in meat products: A new screening method based on high resolution melting analysis of cyt b gene. *Food Chem.* **2017**, *237*, 701-706. [CrossRef]

14. ISO/TC334. Reference Materials—Selected Terms and Definitions; ISO Guide 30:2015; ISO: Geneva, Switzerland, 2015.

15. Yang, Y.; Li, L.; Yang, H.; Li, X.; Zhang, X.; Xu, J.; Zhang, D.; Jin, W.; Yang, L. Development of Certified Matrix-Based Reference Material as a Calibrator for Genetically Modified Rice G6H1 Analysis. *J. Agric. Food Chem.* **2018**, *66*, 3708-3715. [CrossRef] [PubMed]

16. Li, J.; Zhang, L.; Li, L.; Li, X.; Zhang, X.; Zhai, S.; Gao, H.; Li, Y.; Wu, G.; Wu, Y. Development of Genomic DNA Certified Reference Materials for Genetically Modified Rice Kefeng 6. *ACS Omega* **2020**, *5*, 21602-21609. [CrossRef]

17. Li, J.; Li, L.; Zhang, L.; Zhang, X.; Li, X.; Zhai, S.; Gao, H.; Li, Y.; Wu, G.; Wu, Y. Development of a certified genomic DNA reference material for detection and quantification of genetically modified rice KMD. *Anal. Bioanal. Chem.* **2020**, *412*, 7007-7016. [CrossRef]

18. Xu, J.; Qu, S.; Sun, N.; Zhang, W.; Zhang, J.; Song, Q.; Lin, M.; Gao, W.; Zheng, Q.; Han, M.; et al. Construction of a reference material panel for detecting KRAS/NRAS/EGFR/BRAF/MET mutations in plasma ctDNA. *J. Clin. Pathol.* **2021**, *74*, 314-320. [CrossRef] [PubMed]

19. Dong, L.; Wang, X.; Wang, S.; Du, M.; Niu, C.; Yang, J.; Li, L.; Zhang, G.; Fu, B.; Gao, Y.; et al. Interlaboratory assessment of droplet digital PCR for quantification of BRAF V600E mutation using a novel DNA reference material. *Talanta* **2020**, *207*, 120293. [CrossRef]

20. He, H.-J.; Das, B.; Cleveland, M.H.; Chen, L.; Camalier, C.E.; Liu, L.-C.; Norman, K.L.; Fellowes, A.P.; McEvoy, C.R.; Lund, S.P.; et al. Development and interlaboratory evaluation of a NIST Reference Material RM 8366 for EGFR and MET gene copy number measurements. *Clin. Chem. Lab. Med.* **2019**, *57*, 1142-1152. [CrossRef] [PubMed]

21. Vallejo, C.V.; Tere, C.P.; Calderon, M.N.; Arias, M.M.; Leguizamon, J.E. Development of a genomic DNA reference material for Salmonella enteritidis detection using polymerase chain reaction. *Mol. Cell. Probes* **2021**, *55*, 101690. [CrossRef] [PubMed]

22. Børsting, C.; Mas, C.T.; Morling, N. SNP typing of the reference materials SRM 2391b 1-10, K562, XY1, XX74, and 007 with the SNPforID multiplex. *Forensic Sci. Int. Genet.* **2011**, *5*, e81-e82. [CrossRef] [PubMed]

23. Steffen, C.R.; Kiesler, K.M.; Borsuk, L.A.; Vallone, P.M. Beyond the STRs: A comprehensive view of current forensic DNA markers characterized in the PCR-based DNA profiling standard SRM 2391D. *Forensic Sci. Int. Genet. Suppl. Ser.* **2017**, *6*, e426-e427. [CrossRef]

24. Cheng, X.; He, W.; Huang, M. Establishment and Application of Fluorescent Real-Time PCR for the Detection of Duck Meat in Foods. *Food Sci.* **2013**, *34*, 92-96. [CrossRef]

25. SAC/TC387. Identification of Animal Ingredients from Common Livestock and Poultry Real-Time PCR; GB/T38164; AQSIQ of China: Beijing, China, 2019.

26. Tholen, D.W.; Kallner, A.; Kennedy, J.W.; Krouver, J.S.; Meier, K. Evaluation of Precision Performance of Quantitative Measurement Methods; Approved Guideline—Second Edition. NCCLS, 2004. Available online: https://yeec.com/uploadimages1/forum/2013-2/201321114424796191.pdf (accessed on 28 May 2021).

27. ISO/Guide 35:2017(en). Reference Material-Guidance for Characterization and Assessment of Homogeneity and Stability. 2017. Available online: https://www.iso.org/obp/ui/#iso:std:iso:guide:35-ed-4:v1:en (accessed on 20 May 2021).

28. Dong, L.; Meng, Y.; Wang, J.; Liu, Y. Evaluation of droplet digital PCR for characterizing plasmid reference material used for quantifying ammonia oxidizers and denitrifiers. *Anal. Bioanal. Chem.* **2014**, *406*, 1701-1712. [CrossRef] [PubMed]

29. Kim, M.-J.; Lee, Y.-M.; Suh, S.-M.; Kim, H.-Y. Species Identification of Red Deer (Cervus elaphus), Roe Deer (Capreolus capreolus), and Water Deer (*Hydropotes inermis*) Using Capillary Electrophoresis-Based Multiplex PCR. *Foods* **2020**, *9*, 982. [CrossRef]

30. Seddaoui, N.; Amine, A. Smartphone-based competitive immunoassay for quantitative on-site detection of meat adulteration. *Talanta* **2021**, *230*, 122346. [CrossRef]
Article

# Gas Chromatography–Mass Spectrometry-Based Metabolite Profiling for the Assessment of Freshness in Gilthead Sea Bream (*Sparus aurata*)

Athanasios Mallouchos *, Theano Mikrou and Chrysavgi Gardeli

Department of Food Science and Human Nutrition, Agricultural University of Athens, Iera Odos 75, 118 55 Athens, Greece; theano.mikrou@gmail.com (T.M.); agardeli@aua.gr (C.G.)

* Correspondence: amallouchos@aua.gr; Tel.: +30-210-529-4581

Received: 10 March 2020; Accepted: 6 April 2020; Published: 9 April 2020

**Abstract:** Gilthead sea bream (*Sparus aurata*) is one of the most important farmed Mediterranean fish species, and there is considerable interest for the development of suitable methods to assess its freshness. In the present work, gas chromatography–mass spectrometry-based metabolomics was employed to monitor the hydrophilic metabolites of sea bream during storage on ice for 19 days. Additionally, the quality changes were evaluated using two conventional methods: sensory evaluation according to European Union’s grading scheme and K-value, the most widely used chemical index of fish spoilage. With the application of chemometrics, the fish samples were successfully classified in the freshness categories, and a partial least squares regression model was built to predict K-value. A list of differential metabolites were found, which were distinguished according to their evolution profile as potential biomarkers of freshness and spoilage. Therefore, the results support the suitability of the proposed methodology to gain information on seafood quality.

**Keywords:** sea bream; fish; spoilage; metabolomics; multivariate analysis; biomarkers

## 1. Introduction

Gilthead sea bream (*Sparus aurata*) is farmed intensively in Greece and accounts for over half of all production in Europe. In 2018, the volume of production reached 61,000 tons, with a value of EUR 276 million. Greece in particular is expected to double its production by 2030 in order to meet the growing demand and maintain its market position globally [1].

Fish quality is objectively the most important characteristic that affects acceptance by the consumer, and it is dependent on a wide range of factors [2]. Freshness (or degree of spoilage) is a decisive factor in assessing fish quality. Its deterioration begins immediately after slaughter and takes place through biochemical, physicochemical, and microbiological mechanisms [3]. Post-mortem changes depend on species, age, diet, slaughter method, processing, and conditions during transportation and storage, such as temperature, which is the most important factor affecting the commercial life of the product [4]. Preservation on ice is the most common method of maintaining fresh fish, which limits microbial growth, the main cause of spoilage.

The European Community has established common marketing standards for fishery products to assess freshness through organoleptic examination [5]. Thus, the fishing industry must grade the products in three freshness categories defined as Extra, A, and B. Fish not classified in any of these grades are considered unacceptable. Although organoleptic examination is still the most satisfactory way of assessing fish freshness, issues of objectivity and convenience can be claimed if compared with instrumental methods.

From an analytical point of view, several methods have been recommended in order to evaluate fish quality, which rely on the determination of chemical, microbiological, and physical parameters [2].
The K-value, one of the most widely used chemical indexes to monitor fish quality, is based on the measurement of adenosine triphosphate (ATP) and its degradation products, namely, adenosine diphosphate (ADP), adenosine monophosphate (AMP), inosine phosphate (IMP), inosine (INO), and hypoxanthine (Hx) [6]. However, the K-value is subject to large inter- and intra-species variations and is dependent on many factors [7].

The term “metabolomics” refers to the systematic study of low molecular mass metabolites, which vary under a given set of conditions in the cell, tissue, or organism [8]. In recent years, the metabolomics studies on seafood products have been steadily increasing and have focused mainly on the nutritional status of fish [9], differentiation between wild and farmed fish [10–12], and classification according to aquaculture system [13–15], but also some steps have been taken towards seafood freshness. More specifically, the changes in metabolic profiles during cold storage have been investigated on yellowtail [16], bogue [17], mussels [18], and salmon [19]. With regards to sea bream, Picone et al. [20] investigated the molecular profiles using ¹H-NMR only at the beginning and the end of iced storage of fish produced with different aquaculture systems. Heude and co-workers [21] proposed a method based on NMR spectroscopy for the rapid determination of K-value and trimethylamine nitrogen content on sea bream, among other fish.

In the present study, gas chromatography–mass spectrometry (GC–MS) was used to monitor the changes in the polar metabolite fraction of sea bream during storage on ice, in order to identify potential markers of freshness and spoilage. Multivariate data analysis was applied to classify fish samples in freshness categories according to EU sensory scheme, and a partial least squares regression (PLS-R) model was built to predict K-value.

## 2. Materials and Methods

### 2.1. *Fish Provision, Storage, and Sampling*

Gilthead sea bream samples (400–600 g, 25–30 cm) were obtained directly from a Greek fish processing plant (PLAGTON S.A., Mitikas, Aitoloakarnania, Greece). Fish were farmed in cages in the geographical area designated as Food and Agriculture Organization (FAO) 37.2.2 (Ionian Sea) and were slaughtered by immersion in ice cold water (hypothermia), packed with flaked ice into self-draining polystyrene boxes, and delivered to the laboratory within 3–4 h of harvesting. Two fish batches, each consisting of 30 ungutted whole fish, were used in the course of two independent storage trials. The fish batches were harvested in April and July of the same year. The fish samples were stored in a refrigerator, and fresh ice was added daily. The harvesting day was considered as day 0 of storage period. The sampling began the next day of storage (day 1), and afterwards continued every 2 days for a total period of 19 days. At each sampling point, three randomly chosen fish were removed from the batch and used for the subsequent analyses.

### 2.2. *Sample Preparation for GC–MS Metabolomics*

Fish were treated as described in Association of Official Agricultural Chemists (AOAC) Official Method 937.07. The heads, scales, tails, fins, guts, and inedible bones were removed and discarded. Then, fish were filleted to obtain all flesh and skin from head to tail and from top of back to belly on left side only. Each fillet (white muscle with skin) was cut quickly in small cubes and snap-frozen in liquid nitrogen to quench the metabolism. Tissue grinding was performed in a pre-cooled A11 analytical mill (IKA, Wilmington, NC, USA) to obtain a fine frozen powder. The mill was operated in pulse mode for 10–15 s per grinding batch in order to prevent the thawing of the sample. Aliquots (50 mg) of each powdered sample were accurately weighed (± 0.1 mg) into 2 mL Eppendorf tubes with O-ring screw caps (Sarstedt, Germany) and transferred to –80 °C for storage. The remaining quantity of each sample powder was stored at –80 °C in sealed bags and used for the determination of K-value. Furthermore, from each sampling point, a suitable quantity (10 g) of fish powder was pooled to obtain
a single quality control (QC) sample, which was further processed similarly to unknown samples, as described below.

Tissue disruption and subsequent metabolite extraction was undertaken using a Tissuelyser LT (Qiagen, Germantown, MD, USA) according to a modified Bligh and Dyer method [22]. Pre-chilled and degassed homogenization solvent (525 µL methanol/water, 2:0.625 v/v, HPLC grade), internal standard (50 µL glycine-d5, 0.2 mg/mL in 0.1 M HCl), and two stainless-steel balls (2.5 mm diameter) were added to each Eppendorf tube and, subsequently, the fish powder was homogenized for 2 min at 20 Hz. Then, 200 µL of chloroform was added, and the homogenization was repeated for 1 min. Then, another 200 µL of chloroform was added to each tube, and the contents were mixed for 10 min using a cell shaker. During this process, the samples were always kept on ice. Finally, 200 µL HPLC-grade water was added, and the samples were vortex mixed for 15 s. To initiate phase separation, the samples were centrifuged for 2 min at 12,000 rpm. A total of 100 µL of the aqueous fraction was transferred in new Eppendorf tubes with pre-punctured screw caps and lyophilized overnight (12 h). After replacing the caps with new ones, the sample pellets were stored at -80 °C until required for analysis.

Prior to GC-MS analysis, a two-stage chemical derivatization process was carried out to impart volatility to non-volatile metabolites, while also enabling thermal stability [23,24]. The lyophilized samples were left to reach room temperature for 15 min and then 40 µL of 20 mg/mL methoxyamine solution in pyridine (Acros Organics, Geel, Belgium) was added, and they were then incubated at 30 °C for 90 min in an orbital heating block. Subsequently, 70 µL MSTFA (N-methyl-N-trimethylsilyltrifluoroacetamide—Acros Organics) was added, and the samples were incubated at 37 °C for 90 min. After cooling the samples for 5 min, 20 µL of retention index solution (n-alkanes C10–C24, 0.6 mg/mL in pyridine—Sigma Aldrich, Darmstadt, Germany) was added and the contents were transferred to 200 µL conical insert placed in 2 mL vial with screw cap for further GC-MS analysis.

## 2.3. GC-MS Analysis

GC-MS analysis was carried out using a Shimadzu GCMS QP-2010 Ultra operated with the accompanied GCMS Solution software. Helium was used as a carrier gas at a constant linear velocity of 36 cm/s. Sample injections (1 µL) were performed with AOC 20 s autosampler in split mode (split ratio 1/25). The temperature of the injection port, interface, and ion source was set at 250, 290, and 230 °C, respectively. Separation of compounds was carried out in a MEGA-5HT fused silica capillary column (30 m × 0.25 mm, 0.25 µm film thickness, MEGA S.r.l., Legnano, Italy). Oven temperature was maintained initially at 60 °C for 1 min, then programmed at 10 °C/min to 325 °C, and held for 5 min. The mass spectrometer was operated in electron ionization mode with the electron energy set at 70 eV and a scan range of 70–600 m/z. The samples (QC and blanks included) were analyzed in a predetermined order [24].

## 2.4. Data Processing Procedure

Raw data were processed with MS-DIAL software, which is freely available at the PRIMe website ([http://prime.psc.riken.jp/](http://prime.psc.riken.jp/)) [25]. Metabolite identification was performed according to the Metabolomics Standards Initiative at four levels [26]:

* MSI level 1 (identified compounds): based on similarity of retention index (RI) and mass spectrum relative to an authentic compound analyzed under identical experimental conditions.

* MSI level 2 (putatively annotated compounds): agreement of retention index (ΔRI < 20) and mass spectrum (match > 850) coming from the publicly available libraries at PRIMe. Amdis (v. 2.72) and NIST MS Search software (v. 2.2) including NIST 14 mass spectral library were used complimentarily.

* MSI level 3 (putatively characterized compound classes): agreement of RI or mass spectrum to known compounds of a chemical class.

* MSI level 4 (unknown compounds).
The resulting output from this procedure was a retention index vs. sample data matrix with related metabolite IDs and peak heights linked to each sample injection. Subsequently, manual data curation was performed, which included the removal of metabolic features detected in < 50% of QC samples; the combination of metabolite rows that had two or more identical peaks, such as sugars; and normalization to sample mass used in extraction. Finally, the data were normalized to the QC samples using a low-order nonlinear locally estimated smoothing function (LOESS) [24], in order to correct for the signal drift within and between analytical blocks. Afterwards, metabolites with relative standard deviation (RSD) > 30% within pooled QCs were removed. The final data matrix was further processed statistically in MetaboAnalyst 4.0 web-based tool suite [27]. This included multivariate and univariate testing as detailed in the Results and Discussion section. Before statistical processing, the data were log-transformed and mean centered. Partial least squares regression (PLS-R) was performed using The Unscrambler X ver. 10.4 (CAMO Software AS, Oslo, Norway).

## 2.5. Freshness Assessment

The freshness rating of raw fish was performed by a panel of three trained assessors according to the European Union’s grading system [5]. This system distinguishes between three freshness categories (Extra, A, B) corresponding to various levels of spoilage. Category E corresponds to the highest quality level, followed by categories A and B, whereas fish graded below B is considered unacceptable for consumption. In order to rate this evaluation, a 0–3 score scale was used (rating of categories: Extra ≥ 2.7, 2 ≤ A < 2.7, 1 ≤ B < 2, unacceptable < 1) according to [28].

## 2.6. ATP Breakdown Products

ATP and its degradation products (ADP, AMP, IMP, Ino, and Hx) were isolated from fish tissue according to Ryder [29]. Chromatography was performed using a JASCO HPLC system (JASCO International Co., Ltd., Tokyo, Japan) consisting of a quaternary pump (PU-2089 Plus), an autosampler (AS-1555), and a photodiode array detector (MD-910). The separation was accomplished with a Luna C18 column (250 mm × 4.6 mm i.d., 5 µm; Phenomenex, Torrance, CA, USA) using gradient elution. Mobile phase A was a 0.05 M phosphate buffer (pH 7) and mobile phase B was acetonitrile (Sigma Aldrich, Louis, MI, USA). The elution program was as follows: 0 min, 100% A; 9 min, 97% A; 15 min, 85% A; 17 min, 60% A. Final conditions were kept for 7 min and the column was equilibrated for 15 min at initial conditions. The flow rate was set at 1 mL/min and the injection volume was 20 µL. The monitoring wavelength was set at 254 nm and the molar concentration of ATP breakdown products were calculated from their corresponding calibration curves using the external standard method. The K-value (%) was calculated from Equation (1):

$$K(\%) = \frac{(Ino + Hx)}{(ATP + ADP + AMP + IMP + Ino + Hx)} \times 100\% \quad (1)$$

# 3. Results and Discussion

## 3.1. Freshness Assessment Using Classical Methods

Quality deterioration of fish during storage on ice was monitored using a sensory method (EU grading system), and a chemical one (K-value), which is based on the measurement of ATP breakdown products. The changes in sensory score and K-value of sea bream during 19 storage days on ice are shown in Figure 1. As expected, the sensory score decreased linearly ($y = -0.1404x + 3.0057$), showing high negative correlation with storage time ($r = -0.9880, p < 0.001$). Until day 1 of storage, the freshness rating of fish was evaluated as Extra. From day 3 to day 7, the freshness of fish was rated A, whereas the category B was assigned to fish stored between 9–15 days. The limit of acceptability of raw sea bream stored on ice was about 16–17 days. As the sensory quality of fish decreased, the K-value increased linearly ($y = 2.5033x + 3.3305$), showing a high positive correlation with storage time
(r = 0.9974, *p* < 0.001) and a negative correlation with sensory score (r = -0.9809, *p* < 0.001). When the fish was considered unacceptable (day 17), the K-value was 45%. Similar findings have been reported by others authors [6,30,31]. Small variations could be attributed to the different rearing area and farming method among others.

<figure><img src="image_71.png" /><figcaption>**Figure 1.** Changes in sensory score (■-▲-) and K-value (■●-●) of sea bream stored in ice. Each point represent the mean value of six replicate measurements (three fish samples x two storage experiments at each sampling point). Error bars denote standard deviation. The labels of sensory scores denote the freshness category according to EU grading system.</figcaption></figure>

## 3.2. Freshness Assessment Using GC-MS Metabolomics

Before starting the real data elaboration, method performance was evaluated using several quality control criteria. The primary requirement was to check the relative abundance of amino acids and sugar trimethylsilyl (TMS) derivatives in QCs. A detailed description is provided elsewhere [23]. After passing the above criteria, the internal standard performance was checked. Deuterated glycine-d5 was added in every sample (including QCs and blanks) to monitor the extraction procedure. For this reason, the RSD% of peak height was calculated and the value obtained was 5.4% for QCs (n = 15) and 15.3% for samples (n = 54). As a last check, a preliminary principal component analysis (PCA) was obtained with the peak heights of the entire dataset. A tight clustering of the QCs, as well as blank samples, was observed in the score plot (Figure S1), which is a further confirmation of the robustness of the analytical procedure, not only for the internal standard but for the whole fingerprint of fish.

After data curation, the samples were grouped in 10 classes (from 0 to 9) according to the sampling sequence during fish storage on ice (class 0 represents the first day of storage, class 1 the third day, etc.). Principal component analysis (PCA) of data exhibited a significant ability to separate samples according to storage time (Figure 2a). It is evident that very fresh samples (class 0, day 1) were clearly separated from the other classes, and most distinctively from those at the later stages of storage (class 8, 9). This shows that there is quantitative information in these data, as PC1, with an explained variance of 62.5%, is the component that describes the evolution of fish spoilage. This was depicted even more clearly in the plot of PC1 scores values vs. storage time (Figure 2b). Thus, PC1 can condensate all information of the metabolite features and give a measure of the molecular quality of fish.
<figure><img src="image_72.png" /><figcaption>**Figure 2.** (a) Principal component analysis (PCA) score plot derived from the hydrophilic metabolites of sea bream during storage on ice. The legend indicates the sampling sequence (0-1-2-3-4-5-6-7-8-9) that corresponds to storage day (1-3-5-7-9-11-13-15-17-19), respectively; (b) evolution of sea bream spoilage as described by PC1 scores values vs. storage time on ice.</figcaption></figure>

Studying the multivariate loadings values revealed various metabolites that either increased or decreased with fish storage. The most important loadings were established and confirmed by Kruskal-Wallis ANOVA, as well as by performing Spearman's correlation analysis. Figure 3 shows the top 25 highly correlating and significant metabolites ($p < 0.05$) that either increased (shown in red; positive $R$) or decreased during storage (shown in blue; negative $R$).

<figure><img src="image_73.png" /><figcaption>**Figure 3.** Pattern recognition—Spearman's correlation analysis showing the top 25 metabolite features correlated significantly with sampling sequence (0-1-2-3-4-5-6-7-8-9 is equivalent to storage day 1-3-5-7-9-11-13-15-17-19). Each row represents the most significant metabolite identified from the test ($p < 0.05$). The x-axis shows correlation score, whereas the y-axis corresponds to gas chromatography-mass spectrometry (GC-MS) peak number from peak index (see Table 1).</figcaption></figure>
After these encouraging results, the data were grouped into four classes (0 to 3) representing the EU freshness grades (i.e., grade Extra: class 0; grade A: class 1; grade B: class 2; unacceptable: class 3), as described in EC no. 2396/1996. Partial least squares discriminant analysis (PLS-DA) was carried out in order to find a discriminant index of freshness. It is evident that the supervised model (Figure 4a) can clearly classify the samples into the correct freshness grade. Although there was some overlap of the confidence ellipses of grade A (class 1) and B (class 2), the grade Extra (class 0) was further apart from unacceptable samples (class 3). Similarly to the aforementioned PCA results, it seems that this separation was described mainly by PC1, which accounted for the 65.5% of model variance. The optimal number of components, as calculated by 10-fold cross-validation, was 3 (Figure S2). The predictive ability of the model (Q²), accuracy, and coefficient of determination (R²) were satisfactory (0.94, 0.97, and 0.76, respectively). The significance of class discrimination was verified by performing a permutation test (*p* < 0.001; 0/1000), and the performance was measured using group separation distance (B/W ratio) [32] (Figure S3). The variable importance in projection (VIP) scores were also calculated from the PLS-DA, and Figure 4b highlights the top 15 highly significant metabolites that were identified for each freshness grade.

<figure><img src="image_74.png" /><figcaption>**Figure 4.** (a) Partial least squares discriminant analysis (PLS-DA) scores plotted for freshness grades of sea bream stored on ice. The legend indicates the four EU grades: Extra (0), A (1), B (2), unacceptable (3). (b) Top 15 metabolite features based on variable importance in projection (VIP) scores from PLS-DA. The x-axis shows the scores whereas the y-axis corresponds to GC-MS peak number from peak index (see Table 1). Color bars show median intensity of metabolite feature in the respective group.</figcaption></figure>

The combined result of PCA, PLS-DA, and Kruskal-Wallis ANOVA is summarized in Table 1, which shows a list of metabolites significantly correlated with fish storage on ice. They can be distinguished in two groups—in the first group were metabolites whose relative content increased significantly during storage, whereas the second group comprised metabolites with a decreasing trend. Thus, we can infer that the first group of compounds constitute potential markers of spoilage, whereas the second group could be markers of freshness. Their respective evolution pattern during storage is summarized in Tables S1 and S2. The relative concentration of six amino acids (leucine, isoleucine, valine, phenylalanine, tyrosine, methionine) increased during storage on ice as a result of autolysis and bacterial spoilage. Similar findings were observed in bogue [17] and salmon [19]. On the contrary, the observed decrease of glycine and glutamic acid probably indicates that their degradation by microorganisms occurred in a higher rate than their release from muscle proteins. This is in contrast to the aforementioned studies, but can be rationalized by the different spoilage
microorganisms developing in each fish species. The amount of inosine increased during storage, as expected, and was confirmed also by HPLC analysis of ATP breakdown products. The levels of succinic, malic, and fumaric acid, involved in the Krebs cycle, decreased during storage, thus indicating a preferential consumption for bacterial growth. In fact, organic acids or amino acids rather than glucose are the preferred carbon sources for *Pseudomonas* [33], the dominating spoilage genus in sea bream at low temperatures [34]. The increasing trend of some sugars, such as ribose and galactose, has been observed previously in other aquatic products, such as mussels and yellowtail fish [16,18]. On the contrary, the sugar phosphates, such as ribulose 5-phosphate, which is the end product of pentose phosphate pathway (oxidative branch), decreased significantly with storage, and thus may represent freshness markers. It should be noted here that the interpretation of the biological significance of metabolomics data is not always straightforward. The main difficulty arises from the nature of freshness loss, which is a process described primarily by two different phenomena—the autolysis from endogenous enzymes and the spoilage due to microbial growth. In addition, the whole picture is complicated by the evolution of microbial diversity that leads to shifts of metabolite profiles.

<table><caption>**Table 1.** Metabolites that either increased or decreased significantly<fnref n="1" /> during storage of sea bream on ice.</caption><thead><tr><th>Peak Number</th><th>Significant Metabolites</th><th>MSI Level</th><th>Identifier from Relevant Database</th></tr></thead><tbody><tr><td colspan="4" style="text-align: left; font-style: italic;">Increasing trend</td></tr><tr><td>41</td><td>Gluconic acid</td><td>2</td><td>HMDB0000625</td></tr><tr><td>16</td><td>Glyceric acid</td><td>2</td><td>CHEBI:32398</td></tr><tr><td>29</td><td>Ribose</td><td>2</td><td>CHEBI:47014</td></tr><tr><td>37</td><td>Galactose</td><td>1</td><td>CHEBI:4139</td></tr><tr><td>8</td><td>Ethanolamine</td><td>2</td><td>CHEBI:16000</td></tr><tr><td>47</td><td>Inosine</td><td>2</td><td>CHEBI:17596</td></tr><tr><td>9</td><td>Leucine</td><td>1</td><td>CHEBI:25017</td></tr><tr><td>6</td><td>Valine</td><td>1</td><td>CHEBI:16414</td></tr><tr><td>27</td><td>Phenylalanine</td><td>1</td><td>CHEBI:17295</td></tr><tr><td>35</td><td>Tyrosine</td><td>2</td><td>CHEBI:17895</td></tr><tr><td>12</td><td>Isoleucine</td><td>2</td><td>CHEBI:17191</td></tr><tr><td>22</td><td>Methionine</td><td>1</td><td>CHEBI:16643</td></tr><tr><td>11</td><td>Glycerol</td><td>1</td><td>CHEBI:17754</td></tr><tr><td>30</td><td>Ribitol</td><td>2</td><td>CHEBI:15963</td></tr><tr><td colspan="4" style="text-align: left; font-style: italic;">Decreasing trend</td></tr><tr><td>45</td><td>Ribulose 5-phosphate</td><td>2</td><td>CHEBI:17363</td></tr><tr><td>44</td><td>Arabinose-5-phosphate</td><td>2</td><td>CHEBI:16241</td></tr><tr><td>46</td><td>Sugar phosphate_2381<sup>2</sup></td><td>3</td><td>N/A</td></tr><tr><td>32</td><td>Glycerol 3-phosphate</td><td>2</td><td>CHEBI:15978</td></tr><tr><td>15</td><td>Succinic acid</td><td>1</td><td>CHEBI:15741</td></tr><tr><td>17</td><td>Fumaric acid</td><td>1</td><td>CHEBI:18012</td></tr><tr><td>21</td><td>Malic acid</td><td>1</td><td>HMDB0000744</td></tr><tr><td>2</td><td>Lactic acid</td><td>1</td><td>CHEBI:78320</td></tr><tr><td>5</td><td>a-Aminobutyric acid</td><td>2</td><td>CHEBI:35621</td></tr><tr><td>24</td><td>Creatinine</td><td>2</td><td>CHEBI:16737</td></tr><tr><td>1</td><td>Methylamine</td><td>2</td><td>CHEBI:16830</td></tr><tr><td>14</td><td>Glycine</td><td>1</td><td>CHEBI:15428</td></tr><tr><td>26</td><td>Glutamic acid</td><td>2</td><td>HMDB00148</td></tr></tbody></table>

The analytes of Table 1 were used further in PLS-R as input variables (predictors, X) in order to predict K-value (output variable, Y). Segmented cross validation was employed using 18 segments with three replicate samples each. Thus, each segment corresponded to a sampling point (18 segments = 9 sampling points × 2 fish batches). External validation using an independent test set was not

<footnote n="1">¹ According to the combined results of Kruskal-Wallis ANOVA, Spearman's correlation analysis, and VIP scores from PLS-DA. ² The number denotes the *n*-alkane retention index in MEGA HT-5 column.</footnote>
performed due to the relatively small number of samples (n = 54) under study. The optimal number of factors was 3 and explained the 95% of total variance. The performance metrics and the regression line of the model are presented in Figure 5. The root mean square error (RMSE) and the coefficient of determination (R²) at the validation stage suggested good prediction performance, with their values being 3.4710 (K-value %) and 0.9473, respectively. According to the slope of the regression line (0.9546), there was an almost perfect linear relationship between the predicted and the measured K-values.

<figure><img src="image_75.png" /><figcaption>**Figure 5.** Comparison between the observed and predicted K-values (%) by the partial least squares regression (PLS-R) model based on important metabolites listed in Table 1. The shape of symbols denotes the calibration (●-) and the validation (-▲-) set (black line: the ideal y = x line; dotted lines: ± 3.5% K-value).</figcaption></figure>

**4. Conclusions**

The present study demonstrated for the first time that GC-MS-based metabolomics is an efficient tool to monitor the quality loss of sea bream during storage on ice. We have clearly presented a panel of hydrophilic metabolites linked directly to storage time of fish that could be used as potential markers of freshness and spoilage. Additionally, with the application of multivariate data analysis, the samples were successfully classified to freshness grades, and the K-value was predicted using a PLS-R model. Therefore, our results support the suitability of the proposed methodology to gain information on seafood quality. However, we should note that an approach based on the hydrophilic fraction of metabolites solely, as described here, is not enough for the complete elucidation of the post-mortem changes occurring at the molecular level. Hence, future applications should include the investigation of the lipophilic metabolites using larger-scale storage experiments in combination with microbiological and sensorial data.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2304-8158/9/4/464/s1, Figure S1: PCA score plot for the entire data set before data curation. The legend indicates the type of samples (BL: blank, QC: quality control, S: sample). Figure S2: PLS-DA classification using different number of components. The red star indicates the best classifier based on Q2. Figure S3: PLS-DA model validation by permutation tests based on group separation distance (B/W). Table S1: Most significant metabolites with an increasing trend during storage of sea bream on ice. The y-axis represents the relative amount before and after variable transformation, respectively. The x-axis shows sampling sequence that is equivalent to storage time. Table S2: Most significant metabolites with a decreasing trend during storage of sea bream on ice. The y-axis represents the relative amount before and after variable transformation, respectively, for the bar plot and the box plot. The x-axis shows sampling sequence that is equivalent to storage time.
**Author Contributions:** Conceptualization, A.M.; methodology, A.M.; formal analysis, A.M.; investigation, T.M. and C.G.; data curation, A.M.; writing—original draft preparation, A.M., T.M., and C.G.; writing—review and editing, A.M., T.M., and C.G.; supervision, A.M. and C.G.; project administration, C.G. All authors have read and agreed to the published version of the manuscript.

**Acknowledgments:** The authors would like to thank PLAGTON S.A. (Aitoloakarnania, Greece) for the provision of fish.

**Conflicts of Interest:** The authors declare no conflict of interest.

## References

1. Federation of Greek Maricultures. Available online: https://fgm.com.gr/uploads/file/FGM_19_ENG_WEB_spreads.pdf (accessed on 7 December 2019).

2. Botta, J.R. *Evaluation of Seafood Freshness Quality*, 1st ed.; VCH Publishers, Inc.: New York, NY, USA, 1995; ISBN 978-0-471-18580-2.

3. Dalgaard, P. FISH Spoilage of Seafood. In *Encyclopedia of Food Sciences and Nutrition*; Caballero, B., Trugo, L., Finglas, P.M., Eds.; Academic Press: Cambridge, UK, 2003; pp. 2462–2471. ISBN 978-0-12-227055-0.

4. Olafsdottir, G.; Nesvadba, P.; Di Natale, C.; Careche, M.; Oehlenschläger, J.; Tryggvadóttir, S.V.; Schubring, R.; Kroeger, M.; Heia, K.; Esaiassen, M.; et al. Multisensor for fish quality determination. *Trends Food Sci. Technol.* **2004**, *15*, 86–93. [CrossRef]

5. Council Regulation (EC) No 2406/96 of 26 November 1996. Laying down common marketing standards for certain fishery products. *Off. J. Eur. Communities* **1996**, *L334*, 1–15.

6. Alasalvar, C.; Taylor, K.D.A.; Öksüz, A.; Garthwaite, T.; Alexis, M.N.; Grigorakis, K. Freshness assessment of cultured sea bream (Sparus aurata) by chemical, physical and sensory methods. *Food Chem.* **2001**, *72*, 33–40. [CrossRef]

7. Huidobro, A.; Pastor, A.; Tejada, M. Adenosine Triphosphate and Derivatives as Freshness Indicators of Gilthead Sea Bream (Sparus aurata). *Food Sci. Technol. Int.* **2001**, *7*, 23–30.

8. Dunn, W.B.; Ellis, D.I. Metabolomics: Current analytical platforms and methodologies. *TrAC Trends Anal. Chem.* **2005**, *24*, 285–294.

9. Gil-Solsona, R.; Nácher-Mestre, J.; Lacalle-Bergeron, L.; Sancho, J.V.; Calduch-Giner, J.A.; Hernández, F.; Pérez-Sánchez, J. Untargeted metabolomics approach for unraveling robust biomarkers of nutritional status in fasted gilthead sea bream (*Sparus aurata*). *PeerJ* **2017**, *5*, e2920. [CrossRef]

10. Melis, R.; Cappuccinelli, R.; Roggio, T.; Anedda, R. Addressing marketplace gilthead sea bream (Sparus aurata L.) differentiation by 1H NMR-based lipid fingerprinting. *Food Res. Int.* **2014**, *63*, 258–264. [CrossRef]

11. Rezzi, S.; Giani, I.; Héberger, K.; Axelson, D.E.; Moretti, V.M.; Reniero, F.; Guillou, C. Classification of Gilthead Sea Bream (*Sparus aurata*) from ¹H NMR Lipid Profiling Combined with Principal Component and Linear Discriminant Analysis. *J. Agric. Food Chem.* **2007**, *55*, 9963–9968. [CrossRef]

12. Mannina, L.; Sobolev, A.P.; Capitani, D.; Iaffaldano, N.; Rosato, M.P.; Ragni, P.; Reale, A.; Sorrentino, E.; D'Amico, I.; Coppola, R. NMR metabolic profiling of organic and aqueous sea bass extracts: Implications in the discrimination of wild and cultured sea bass. *Talanta* **2008**, *77*, 433–444. [CrossRef]

13. Savorani, F.; Picone, G.; Badiani, A.; Fagioli, P.; Capozzi, F.; Engelsen, S.B. Metabolic profiling and aquaculture differentiation of gilthead sea bream by 1H NMR metabonomics. *Food Chem.* **2010**, *120*, 907–914. [CrossRef]

14. Melis, R.; Anedda, R. Biometric and metabolic profiles associated to different rearing conditions in offshore farmed gilthead sea bream (*Sparus aurata L.): General. Electrophoresis* **2014**, *35*, 1590–1598. [CrossRef] [PubMed]

15. Del Coco, L.; Papadia, P.; De Pascali, S.; Bressani, G.; Storelli, C.; Zonno, V.; Fanizzi, F.P. Comparison among Different Gilthead Sea Bream (Sparus aurata) Farming Systems: Activity of Intestinal and Hepatic Enzymes and 13C-NMR Analysis of Lipids. *Nutrients* **2009**, *1*, 291–301. [CrossRef] [PubMed]

16. Mabuchi, R.; Adachi, M.; Ishimaru, A.; Zhao, H.; Kikutani, H.; Tanimoto, S. Changes in Metabolic Profiles of Yellowtail (Seriola quinqueradiata) Muscle during Cold Storage as a Freshness Evaluation Tool Based on GC-MS Metabolomics. *Foods* **2019**, *8*, 511. [CrossRef] [PubMed]

17. Ciampa, A.; Picone, G.; Laghi, L.; Nikzad, H.; Capozzi, F. Changes in the Amino Acid Composition of Bogue (Boops boops) Fish during Storage at Different Temperatures by 1H-NMR Spectroscopy. *Nutrients* **2012**, *4*, 542–553. [CrossRef]
18. Aru, V.; Pisano, M.B.; Savorani, F.; Engelsen, S.B.; Cosentino, S.; Cesare Marincola, F. Metabolomics analysis of shucked mussels' freshness. *Food Chem.* **2016**, *205*, 58–65. [CrossRef]

19. Shumilina, E.; Ciampa, A.; Capozzi, F.; Rustad, T.; Dikiy, A. NMR approach for monitoring post-mortem changes in Atlantic salmon fillets stored at 0 and 4°C. *Food Chem.* **2015**, *184*, 12–22. [CrossRef]

20. Picone, G.; Balling Engelsen, S.; Savorani, F.; Testi, S.; Badiani, A.; Capozzi, F. Metabolomics as a Powerful Tool for Molecular Quality Assessment of the Fish Sparus aurata. *Nutrients* **2011**, *3*, 212–227. [CrossRef]

21. Heude, C.; Lemasson, E.; Elbayed, K.; Piotto, M. Rapid Assessment of Fish Freshness and Quality by 1H HR-MAS NMR Spectroscopy. *Food Anal. Methods* **2015**, *8*, 907–915. [CrossRef]

22. Wu, H.; Southam, A.D.; Hines, A.; Viant, M.R. High-throughput tissue extraction protocol for NMR- and MS-based metabolomics. *Anal. Biochem.* **2008**, *372*, 204–212. [CrossRef]

23. Fiehn, O. Metabolomics by gas chromatography-mass spectrometry: Combined targeted and untargeted profiling. *Curr. Protoc. Mol. Biol.* **2016**, *114*, 30–34. [CrossRef]

24. The Human Serum Metabolome (HUSERMET) Consortium; Dunn, W.B.; Broadhurst, D.; Begley, P.; Zelena, E.; Francis-McIntyre, S.; Anderson, N.; Brown, M.; Knowles, J.D.; Halsall, A.; et al. Procedures for large-scale metabolic profiling of serum and plasma using gas chromatography and liquid chromatography coupled to mass spectrometry. *Nat. Protoc.* **2011**, *6*, 1060–1083. [CrossRef] [PubMed]

25. Lai, Z.; Tsugawa, H.; Wohlgemuth, G.; Mehta, S.; Mueller, M.; Zheng, Y.; Ogiwara, A.; Meissen, J.; Showalter, M.; Takeuchi, K.; et al. Identifying metabolites by integrating metabolome databases with mass spectrometry cheminformatics. *Nat. Methods* **2018**, *15*, 53–56. [CrossRef] [PubMed]

26. Sumner, L.W.; Amberg, A.; Barrett, D.; Beale, M.H.; Beger, R.; Daykin, C.A.; Fan, T.W.-M.; Fiehn, O.; Goodacre, R.; Griffin, J.L.; et al. Proposed minimum reporting standards for chemical analysis: Chemical Analysis Working Group (CAWG) Metabolomics Standards Initiative (MSI). *Metabolomics* **2007**, *3*, 211–221. [CrossRef]

27. Chong, J.; Wishart, D.S.; Xia, J. Using MetaboAnalyst 4.0 for Comprehensive and Integrative Metabolomics Data Analysis. *Curr. Protoc. Bioinforma.* **2019**, *68*, e86. [CrossRef] [PubMed]

28. Council Regulation (EEC) No 103/76 of 19 January 1976. Laying down common marketing standards for certain fresh or chilled fish. *Off. J. Eur. Communities* **1976**, *L20*, 29–34.

29. Ryder, J.M. Determination of adenosine triphosphate and its breakdown products in fish muscle by high-performance liquid chromatography. *J. Agric. Food Chem.* **1985**, *33*, 678–680. [CrossRef]

30. Lougovois, V.P.; Kyranas, E.R.; Kyrana, V.R. Comparison of selected methods of assessing freshness quality and remaining storage life of iced gilthead sea bream (Sparus aurata). *Food Res. Int.* **2003**, *36*, 551–560. [CrossRef]

31. Šimat, V.; Bogdanović, T.; Krželj, M.; Soldo, A.; Maršić-Lučić, J. Differences in chemical, physical and sensory properties during shelf life assessment of wild and farmed gilthead sea bream (Sparus aurata, L.): Shelf life parameters of wild and farmed gilthead sea bream. *J. Appl. Ichthyol.* **2012**, *28*, 95–101. [CrossRef]

32. Bijlsma, S.; Bobeldijk, I.; Verheij, E.R.; Ramaker, R.; Kochhar, S.; Macdonald, I.A.; van Ommen, B.; Smilde, A.K. Large-Scale Human Metabolomics Studies: A Strategy for Data (Pre-) Processing and Validation. *Anal. Chem.* **2006**, *78*, 567–574. [CrossRef]

33. Rojo, F. Carbon catabolite repression in Pseudomonas: Optimizing metabolic versatility and interactions with the environment. *FEMS Microbiol. Rev.* **2010**, *34*, 658–684. [CrossRef]

34. Parlapani, F.F.; Mallouchos, A.; Haroutounian, S.A.; Boziaris, I.S. Microbiological spoilage and investigation of volatile profile during storage of sea bream fillets under various conditions. *Int. J. Food Microbiol.* **2014**, *189*, 153–163. [CrossRef] [PubMed]

<figure><img src="image_76.png" /></figure>

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
<figure><img src="image_77.png" /></figure>
Article

# Upstream Regulator Analysis of Wooden Breast Myopathy Proteomics in Commercial Broilers and Comparison to Feed Efficiency Proteomics in Pedigree Male Broilers

Walter G. Bottje <fnref n="1" />/*, Kentu R. Lassiter <fnref n="1" />, Vivek A. Kuttappan <fnref n="2" />, Nicholas J. Hudson <fnref n="3" />, Casey M. Owens <fnref n="1" />, Behnam Abasht <fnref n="4" />, Sami Dridi <fnref n="1" /> and Byungwhi C. Kong <fnref n="1" />

**Citation:** Bottje, W.G.; Lassiter, K.R.; Kuttappan, V.A.; Hudson, N.J.; Owens, C.M.; Abasht, B.; Dridi, S.; Kong, B.C. Upstream Regulator Analysis of Wooden Breast Myopathy Proteomics in Commercial Broilers and Comparison to Feed Efficiency Proteomics in Pedigree Male Broilers.
*Foods* **2021**, *10*, 104. https://doi.org/10.3390/foods10010104

Received: 2 November 2020
Accepted: 23 December 2020
Published: 6 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_78.png" /></figure>

Copyright: © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

**Abstract:** In an effort to understand the apparent trade-off between the continual push for growth performance and the recent emergence of muscle pathologies, shotgun proteomics was conducted on breast muscle obtained at ~8 weeks from commercial broilers with wooden breast (WB) myopathy and compared with that in pedigree male (PedM) broilers exhibiting high feed efficiency (FE). Comparison of the two proteomic datasets was facilitated using the overlay function of Ingenuity Pathway Analysis (IPA) (Qiagen, CA, USA). We focused on upstream regulator analysis and disease-function analysis that provides predictions of activation or inhibition of molecules based on (a) expression of downstream target molecules, (b) the IPA scientific citation database. Angiopoepsilon 2 (ANGPT2) exhibited the highest predicted activation Z-score of all molecules in the WB dataset, suggesting that the proteomic landscape of WB myopathy would promote vascularization. Overlaying the FE proteomics data on the WB ANGPT2 upstream regulator network presented no commonality of protein expression and no prediction of ANGPT2 activation. Peroxisome proliferator coactivator 1 alpha (PGC1α) was predicted to be inhibited, suggesting that mitochondrial biogenesis was suppressed in WB. PGC1α was predicted to be activated in high FE pedigree male broilers. Whereas RICTOR (rapamycin independent companion of mammalian target of rapamycin) was predicted to be inhibited in both WB and FE datasets, the predictions were based on different downstream molecules. Other transcription factors predicted to be activated in WB muscle included epidermal growth factor (EGFR), X box binding protein (XBP1), transforming growth factor beta 1 (TGFB1) and nuclear factor (erythroid-derived 2)-like 2 (NFE2L2). Inhibitions of aryl hydrocarbon receptor (AHR), AHR nuclear translocator (ARNT) and estrogen related receptor gamma (ESRRG) were also predicted in the WB muscle. These findings indicate that there are considerable differences in upstream regulators based on downstream protein expression observed in WB myopathy and in high FE PedM broilers that may provide additional insight into the etiology of WB myopathy.

**Keywords:** wooden breast; myopathy; proteomics; feed efficiency; upstream regulator analysis

## 1. Introduction

The incidence of wooden breast (WB) muscle myopathy has increased in the last several years, possibly in response to genetic selection for growth performance and a shift to heavier market weight birds [1–3]. There have been several reports recently that have investigated WB myopathy to identify genetic signatures and histological defects in the progression of the disease (e.g., [3–9]). Collectively, these studies point toward;

<footnote n="1">¹ Department of Poultry Science & The Center of Excellence for Poultry Science, University of Arkansas, 1260 W. Maple, Fayetteville, AR 72701, USA; klassit@uark.edu (K.R.L.); cmowens@uark.edu (C.M.O.); dridi@uark.edu (S.D.); bkong@uark.edu (B.C.K.)</footnote>

<footnote n="2">² Novus International, 20 Research Park, St. Charles, MO 63304, USA; Vivek.Kuttappan@novusint.com</footnote>

<footnote n="3">³ School of Agriculture and Food Sciences, The University of Queensland, Gatton, Brisbane, QLD 4072, Australia; n.hudson@uq.edu.au</footnote>

<footnote n="4">⁴ Department of Animal and Food Sciences, University of Delaware, Newark, DE 19716, USA; abasht@Udel.edu</footnote>

<footnote>* Correspondence: wbottje@uark.edu</footnote>
(a) phlebitis and inflammation, (b) oxidative stress and metabolic dysfunction, (c) myofiber
degeneration, (d) lipid deposition and (e) development of hard pectoral muscle, particularly
in the cranial region. All of these factors contribute to the poor muscle quality and lower
shelf life of processed breast muscle fillets [10].

We have conducted a shotgun proteomic analysis on breast muscle tissue from pedi-
gree male (PedM) broilers phenotyped for feed efficiency (FE) [11] and from commercial
broilers with and without the presence of severe WB muscle myopathy [2]. It has been
hypothesized that the increased incidence of WB might be related to the increased age
at market weight as well as faster growth rates in broilers. Thus, there is a concern that
increased selection for growth performance might be contributing to the increased inci-
dence of WB myopathy. While WB is an area with unknowns, damaged mitochondria,
low capillarity, lactic acidosis and hypoxia are all suggestive of an aerobic energy supply
issue. Ingenuity Pathway Analysis (IPA) software was used to facilitate the organization
and interpretation of the proteomic data in the studies by Kong et al. [11] and Kuttappan
et al. [2]. A powerful feature of the IPA program is the ability to compare datasets using an
overlay function in which data from one study can be projected (overlaid) onto another
dataset to reveal similarities and differences between the two datasets. It can also be used to
connect upstream regulators to diseases and functions. This overlay function can produce
hypotheses for future studies and insight into fundamental mechanisms between these
global expression datasets. Therefore, the intent of the current study is, firstly, to conduct
upstream regulator analysis of WB data not previously reported in Kuttappan et al. [2],
and secondly, to conduct comparisons of the WB data with the proteomic data associated
with feed efficiency reported by Kong et al. [11].

**2. Materials and Methods**

2.1. *Ethics Statement*

The present study was conducted in accordance with the recommendations in the
guide for the care and use of laboratory animals of the National Institutes of Health. All
procedures for animal care complied with the University of Arkansas Institutional Animal
Care and Use Committee (IACUC): Protocol Nos. 14,012 and 17,080.

2.2. *Wooden Breast Myopathy Samples in Commercial Broilers*

Breast muscle samples for wooden breast (WB) proteomics were obtained from male broilers (n = 24) that were randomly selected from a floor pen study at 52 d in a previous study [2]. After weighing, the birds were euthanized using carbon dioxide. The skin over the *Pectoralis major* muscle was excised and the muscle scored for woody breast [12]. Breast muscle samples (~5 g from right cranial region of the *Pectoralis major*) were collected, immediately snap frozen, and stored at -80 °C. Based on the lesions scores from histological analysis conducted in the same region, the muscle samples were categorized for severity of muscle lesions. From the 24 muscle samples, five samples were chosen that exhibited normal (NORM) histology and five others were chosen with severe (SEV) amounts of lesions to be subjected to proteomic analysis [2].

2.3. *Breast Muscle Samples in Pedigree Male Broilers*

Shotgun proteomics was conducted on breast muscle samples obtained from Pedigree
Male (PedM) broilers individually phenotyped for high or low feed efficiency (FE) (n = 4
per group) [11]. These samples had been obtained from the right *Pectoralis muscle* in PedM
broilers between 8–9 weeks that were humanely killed as part of a larger study [13].

**2.4. *Protein Extraction***

In Kong et al. [11], muscle was homogenized in 1.5 mL of 20 mM potassium phosphate
buffer at pH 7.4 using a hand-held Tissue-Tearor (Biospec Products Inc., Bartlesville, OK,
USA) at speeds varying from 5000 rpm to 32,000 rpm. Following homogenization, samples
were centrifuged at 10,000× g, and the supernatant was collected. In Kuttappan et al. [2],
protein samples were precipitated using TCA-acetone with 8 M urea, 100 mM Tris HCl (pH 8.5) with 5 mM Tris 2-carboxyethyl phosphine at room temperature. Following dissolution and reduction, a 1/20th volume of 200 mM iodoacetamide was added and alkylation carried out for 15 min in the dark at room temperature. The sample was then diluted with four volumes of 100 mM TrisHCl, and digested with trypsin overnight at 37 °C. Protein concentrations in the supernatant in both studies were determined using the Bicinchoninic Acid Protein Assay (Sigma Aldrich, St. Louis, MO, USA). NORM and SEV muscle samples were then diluted to a protein concentration of 20 µg/150 µL in phosphate buffer, and the samples were stored at -80 °C until further analysis.

## 2.5. Shotgun Proteomics

Individual extracted proteins were used in shotgun proteomics analysis with in-gel trypsin digestion followed by tandem mass spectrometry (MS/MS) conducted at the University of Arkansas for Medical Science proteomics core lab (UAMS, Little Rock, AR, USA) [11]. Raw data were analyzed by database searching using Masco (Matrix Science, Boston, MA, USA) and UniProtKB database (http://www.uniprot.org/help/uniprotkb) with the result compiled using the Scaffold Program (Proteome Software, Portland, OR, USA).

In Kuttappan et al. [2], digested sample homogenates were acidified with 1% formic acid, and purified by reversed phase chromatography using C18 affinity media (Omix-Agilent, Santa Clara, CA, USA). Each sample was subjected to three replicate analyses for LC-MS/MS using a hybrid-Orbitrap XL mass spectrometer (ThermoFisher Scientific, Waltham, MA, USA) according to Voruganti et al. [14]. Mass spectrometry analysis was conducted in the DNA/Protein Resource Facility (Oklahoma State University, Stillwater, OK, USA).

## 2.6. Upstream Regulator Analysis

Upstream regulator analysis by IPA is based on: (a) the number and degree of differential expression of downstream target molecules in the existing dataset and (b) prior knowledge of relationships between upstream transcriptional regulators and their downstream target molecules obtained in published literature citations that have been curated and stored in the IPA program. Upstream regulator analysis determines how many known targets or regulators are within the user's dataset, and compares each differentially expressed molecule to the reported relationship in the literature. If the observed direction of change is mostly consistent with either activation or inhibition of the transcriptional regulator, then a prediction is made and an activation z score generated that is also based on literature-derived regulation direction (i.e., "activating" or "inhibiting"). Activation z scores >2.0 indicate that a molecule is activated whereas activation z scores of <−2.0 indicate that a target molecule is inhibited. Qualified predictions can also be made with activation z scores between 2.0 and −2.0. The *p*-value of overlap measures whether there is a statistically significant overlap between the dataset molecules and those regulated by an upstream regulator is calculated using Fisher's Exact Test, and significance is attributed to *p*-values < 0.05.

# 3. Results and Discussion

The main focus of this study is on upstream regulator analysis of severe-fulminant WB myopathy. However, as parts of the results and discussion below are concerned with comparing the WB myopathy proteomic data to a dataset obtained in PedM broilers phenotyped for high and low FE, a description of proteomics methods and of the animals from which samples were obtained for the two datasets is warranted.

## 3.1. *Phenotypic Data for Birds in Wooden Breast and Feed Efficiency Studies*

Body weights, wooden breast scores for commercial broilers and growth performance during phenotyping for feed efficiency in PedM broilers from which muscle samples
were obtained for proteomic studies are presented in Table 1. It can be seen that there were significant differences in body weight and wooden breast (WB) myopathy scores in commercial broilers in severe vs. normal muscle. Additionally, although there were no differences in body weights in PedM broilers, there was a significant difference in feed efficiency due to lower gain while consuming the same amount of feed in the low FE PedM phenotype.

<table><caption>**Table 1.** Body weights (BW) and wooden breast (WB) scores in commercial broilers with and without WB, and BW, weight gain (Gain), feed intake (FI) and feed efficiency (FE, gain to feed) in Pedigree Male broilers from which shotgun proteomics was conducted on breast muscle tissue <fnref n="1" />.</caption><thead><tr><th colspan="3">Commercial Broilers <sup>2</sup></th><th colspan="3">Pedigree Male Broilers <sup>3</sup></th></tr><tr><th></th><th>Normal</th><th>Severe WB</th><th></th><th>High FE</th><th>Low FE</th><th></th></tr></thead><tbody><tr><td>Variable</td><td>(n = 5)</td><td>(n = 5)</td><td>p</td><td></td><td></td><td></td></tr><tr><td>BW Kg (52 days)</td><td>3.2 + 0.2</td><td>3.8 + 0.2</td><td>0.031</td><td>BW (49 days)</td><td>3.13 + 0.08</td><td>3.18 + 0.03</td></tr><tr><td>WB Score</td><td>0.2 + 0.2</td><td>1.4 + 0.2</td><td>0.003</td><td>Gain (kg/7 days)</td><td>0.64 + 0.04</td><td>0.47 + 0.04</td></tr><tr><td></td><td></td><td></td><td></td><td>FI (kg/7 days)</td><td>0.99 + 0.06</td><td>1.03 + 0.09</td></tr><tr><td></td><td></td><td></td><td></td><td>FE (G:F)</td><td>0.65 + 0.01</td><td>0.46 + 0.01</td></tr></tbody></table>

## 3.2. Proteomic Dataset Comparison

Proteomics conducted on high and low FE PedM phenotypes in Kong et al. [11] had been obtained eight years earlier in Bottje et al. [13] than for the WB myopathy study by Kuttappan et al. [2]. The timing of sample collection is important for at least two reasons; (a) it was well before WB myopathy had been observed in commercial broilers and (b) the genetic impact of selection of the PedM broiler line would have had sufficient time in the genetic pipeline to reach the level of commercial broilers in the study by Kuttappan et al. [2]. Breast muscle samples were obtained at 52 day in the WB study and between 56 and 63 day (after FE phenotyping 6–7 weeks) in the PedM FE study. It is worth noting that the number of samples in each study are small (n = 4 or 5), and therefore, this paper may run the risk of not being fully representative of the two populations in each study. However, it can also be noted that the difference in phenotypes were highly significant (p = 0.003 and p < 0.0001 for the WB and FE phenotypes, respectively), suggesting that the relative differences in global expression patterns could be considered representative of a larger phenotype population.

Proteomic analysis of the two studies were conducted at different laboratories. For the commercial broiler study [2], LC-MS/MS analysis was conducted using a hybrid LTQ-OrbitrapXL mass spectrometer (ThermoFisher Scientific, Waltham, MA, USA) as described previously [14], but using 40-cm C18 columns developed over a 2-h period with 0 to 40% acetonitrile. In the PedM broiler study [11], extracted individual proteins were subjected to shotgun proteomics analysis by in-gel trypsin digestion and tandem mass spectrometry (MS/MS) at the University of Arkansas Medical Science (UAMS) Proteomics Core Lab (Little Rock, AR). Both studies used the *Gallus* reference proteome downloaded from the UniProtKB (http://www.uniprot.org/help/uniprotkb) database to identify individual proteins from the spectrometric data. The annotation of proteins upstream regulator analysis and downstream functions were conducted using Ingenuity Pathway Analysis (IPA, Qiagen, Redwood City, CA USA).

## 3.3. Upstream Regulators and Functional Analysis

Figure 1 provides a visual guide to assist in understanding subsequent figures and tables in this study. The predictions of activation (orange background, positive activation Z score) or inhibition (blue background, negative activation Z score) of upstream regulators (Table 2) and functions (Table 3) in WB muscle were determined using the Ingenuity

<footnote n="1">¹ Values represent the mean ± SE. ² Kuttappan et al. [2]; ³ Bottje et al. [13]. Visual scoring values of WB (ranging 0–3) for samples selected for proteomic analysis based on histological assessment indicating normal tissue and severe WB myopathy.</footnote>
Pathway Analysis program. These predictions (activation Z scores) were based on: (1) differential protein expression in the dataset and (2) relationships of upstream regulators and downstream target molecules in the IPA scientific literature database.

<figure><img src="image_79.png" /><figcaption>**Figure 1.** A prediction legend for the interpretation of protein expression and upstream regulator and function analysis in tables and figures.</figcaption></figure>

Supplementary Table S1 contains all of the differentially expressed proteins in the WB myopathy datasets that will appear in tables and figures below. Proteins highlighted in green were down-regulated and those in red were upregulated in the WB myopathy muscle, respectively.

In Tables 2 and 3, proteins in red and green (underlined) lettering denote up- and down-regulation, respectively, in WB myopathy. The protein abbreviations in Tables 2 and 3 are defined in Supplementary Table S1 along with the *p* value and fold difference in expression for each protein. Using an overlay function in the IPA program, it is possible to compare the expression of proteins and predictions of upstream regulators and functions in different datasets. Thus, a qualitative comparison of the WB myopathy proteomics data in the present study can be made to the high vs. low FE proteomics data from Kong et al. [11].

Examples of functional networks in Table 2 are presented in Figure 2. In the vasculogenesis function network (Figure 2A), all up-regulated proteins, with the exception of CALR (calreticulin), contributed to the prediction of activation (activation Z score = 2.43, $p = 4.46 \times 10^{-4}$). The yellow dashed inhibitory line (indicating that the relationships of up-regulation of CALR and predicted activation of vasculogenesis are inconsistent with literature citations) for CALR indicates that enhanced vasculogenesis is typically associated with down-regulation of CALR. The down-regulation of ACP1 would contribute to the predicted activation of vasculogenesis because it would reduce the inhibitory effect on this function. On the other hand, the effects of PKM and LDHA were inconsistent with their down-regulation, indicating that increased expression of these proteins would normally be associated with increased vasculogenesis.
<table><caption>**Table 2.** Functions predicted to be activated (orange) or inhibited (blue) in wooden breast (WB) myopathy compared to normal muscle based on differentially expressed (DE) target molecules in commercial broilers.</caption><thead><tr><th>Functions</th><th>Activation z-Score</th><th>p-Value of Overlap</th><th>Differentially Expressed Proteins <sup>1</sup></th></tr></thead><tbody><tr><td>Development of vasculature</td><td>2.45</td><td>1.21 × 10<sup>-3</sup></td><td>ACPL,ANXA2,C3,CALR, CAPNS1,CDH13,CRYAB,CTSBS,FGA,FLNA,HSPA5,IGHM, LDHA, PGK1,PKM,PPIA,RHOA,VIM,YWHAZ</td></tr><tr><td>Vasculogenesis</td><td>2.43</td><td>4.46 × 10<sup>-4</sup></td><td>ACPL,ANXA2,C3,CALR, CAPNS1,CDH13,CRYAB,CTSBS,FLNA,HSPA5,LDHA,PKM,PPIA,RHOA,VIM,YWHAZ</td></tr><tr><td>angiogenesis</td><td>2.30</td><td>1.67 × 10<sup>-3</sup></td><td>ACPL,ANXA2,C3,CALR, CAPNS1,CDH13,CRYAB,CTSBS,FLNA,HSPA5, LDHA,PGK1,PKM,PPIA,RHOA,VIM,YWHAZ</td></tr><tr><td>Necrosis</td><td>-3.33</td><td>2.68 × 10<sup>-11</sup></td><td>ACPL,ANXA2,C3,CALR, CAPNS1,CBR1,CRYAB, CS,CTSBS,DDX3X,DYNC1H1,EEF1A1,EIF3L,FDP5,FGA,FLNA, GAPDH,GLO1,HINT1,HK1,HNRNPH1,HSP90AA1,HSP90B1,HSPA5,IGHM,KLHL40, LDHA,P4HB, PARK7,PDCD6IP, PDIA3, PEBP1,PKM,POSTN,PPIA, PRDX6,PSMC2,RACK1,RAN,RHOA,RPL11,RPL6,RPLP0,RPS11,RPS14,RPS24,RPS3,Rps3a1,RPS6,RPS7,TAGLN2, UBE2V2,VIM,YWHAE,YWHAZ</td></tr><tr><td>Cell death</td><td>-3.05</td><td>2.01 × 10<sup>-9</sup></td><td>ACPL,AK1,ALDOC,ANXA2,C3,CALR, CAPNS1,CBR1,CRYAB, CS,CTSBS,DDX3X,DES,DYNC1H1,EEF1A1,EIF3L, FDP5,FGA,FLNA, GAPDH,GLO1,HINT1,HK1,HNRNPH1,HSP90AA1,HSP90B1,HSPA2,HSPA5,IGHM,KLHL40, LDHA,P4HB, PARK7,PDCD6IP, PDIA3, PEBP1,PKM,POSTN,PPIA, PRDX6,PSMC2,RAB1A,RACK1,RAN,RHOA,RPL11,RPL6,RPLP0,RPS11,RPS14,RPS24,RPS3,Rps3a1,RPS6,RPS7,TAGLN2, UBE2V2,VIM,YWHAE,YWHAZ</td></tr><tr><td>Apoptosis</td><td>-1.98</td><td>4.61 × 10<sup>-6</sup></td><td>ALDOC,ANXA2,C3,CALR, CAPNS1,CRYAB, CS,CTSBS,DDX3X,DES,DYNC1H1,EEF1A1,FLNA, GAPDH,GLO1,HINT1,HK1,HNRNPH1,HSP90AA1,HSP90B1,HSPA2,HSPA5,IGHM,KLHL40, LDHA,P4HB, PARK7,PDCD6IP, PDIA3, PEBP1,PKM,PPIA, PRDX6,RACK1,RHOA,RPLP0,RPS24,RPS3,RPS6,TAGLN2, UBE2V2,VIM,YWHAE,YWHAZ</td></tr></tbody></table>

<footnote n="1">¹ Target proteins in plain red type were up-regulated (p < 0.05, >1.3-fold difference) in WB myopathy compared to normal tissue. Target molecules in bold green and underlined type were down-regulated (p < 0.05, < -1.3-fold difference) in WB myopathy compared to normal tissue.</footnote>
<table><caption>**Table 3.** Upstream regulators predicted to be activated (orange) or inhibited (blue) in wooden breast (WB) myopathy compared to normal muscle based on differentially expressed downstream target molecules in commercial broilers.</caption><thead><tr><th>Upstream Regulator</th><th>Protein Name</th><th>Activation z-Score</th><th>p-Value of Overlap</th><th>Differentially Expressed Downstream Target Proteins<sup>1</sup></th></tr></thead><tbody><tr><td>ANGPT2</td><td>Angiopoecitin 2</td><td>3.08</td><td>7.41 × 10<sup>-9</sup></td><td>CALR,CRYAB,HSP90AA1,HSPA2,HSPA5,PA4HB,PDIA3,PDIA6,POSTN,RPS6</td></tr><tr><td>EGFR</td><td>Epidermal growth factor</td><td>2.62</td><td>1.38 × 10<sup>-3</sup></td><td>CRYAB,HNRNPH1,HSP90B1,HSPA5,POSTN,PPIA,VIM</td></tr><tr><td>XBPI</td><td>X box binding protein 1</td><td>2.41</td><td>8.79 × 10<sup>-3</sup></td><td>CALR,HSP90B1,HSPA5,PDIA3,PDIA4,PDIA6</td></tr><tr><td>TGFB1</td><td>Transforming growth factor beta 1</td><td>2.21</td><td>1.95 × 10<sup>-6</sup></td><td>AK1, ANXA2, CAPNS1,CO17A1, CTSB, DES,EEF1A1,FLNA,FSCN1, GATD3A,GATD3B,HNT1</td></tr><tr><td>NFE2L2</td><td>Nuclear factor (erythroid-derived 2)-like 2</td><td>2.18</td><td>4.3 × 10<sup>-9</sup></td><td>HNRNPH1,HSP90A1,HAPA5,JARS,IGHM,LDHA,PGK1, PLS3,POSTN,RAB1A,RACK1,RHOA,VIM,YWHAE ARF1,CBR1,DSTN,HSP90AA1,HSP90B1, ME1,PDIA3,PDIA4,PDIA6,PGD,RACK1,RAN,RPL18,RPLP0,TPJ1</td></tr><tr><td>RICTOR</td><td>Regulatory-associated protein Independent of mTOR complex 2</td><td>-3.61</td><td>2.06 × 10<sup>-9</sup></td><td>PSMC2,RPL11,RPL18,RPL4,RPL6,RPLP0,RPS11,RPS2,RPS24,RPS3</td></tr><tr><td>PPARGC1A</td><td>Peroxisome proliferator activator receptor gamma-coactivator 1 α</td><td>-2.35</td><td>1.92 × 10<sup>-4</sup></td><td>AK1, C3, CS,GOT2,LDHA,ME1,PGAM1</td></tr><tr><td>AHR</td><td>Aryl hydrocarbon receptor</td><td>-2.22</td><td>2.02 × 10<sup>-3</sup></td><td>ABC6,CBR1,GOT1, HSP90AA1,HSPA5,IGHM,VIM</td></tr><tr><td>ARNT</td><td>Aryl hydrocarbon receptor nuclear translocator</td><td>-2.20</td><td>3.35 × 10<sup>-6</sup></td><td>ALDOC,GAPDH,IGHM,LDHA,PGK1,TPJ1,VIM</td></tr><tr><td>ESRRG (NR3B3)</td><td>Estrogen related receptor gamma</td><td>-2.17</td><td>1.89 × 10<sup>-6</sup></td><td>ALDOC,GAPDH,LDHA,PKM,TPJ1</td></tr></tbody></table>

<footnote n="1">¹ Target molecules in plain red type were up-regulated (p < 0.05) >1.3-fold difference) in WB myopathy compared to normal tissue. Target molecules in bold green and underlined type were down-regulated (p < 0.05, < -1.3-fold difference) in WB myopathy compared to normal tissue. ANGPT2: Angiopoecitin 2; EGFR: epidermal growth factor; XBPI: X box binding protein; TGFB1: transforming growth factor beta 1; NFE2L2: nuclear factor (erythroid-derived 2)-like 2.</footnote>
<figure><img src="image_80.png" /><figcaption>**Figure 2.** Network of differentially expressed proteins in the wooden breast (WB) myopathy dataset used in calculating the activation Z score for vasculogenesis (A). Overlay of the feed efficiency (FE) proteomic dataset from Kong et al. [11]; (B) indicates no commonality of differentially expressed proteins in WB vasculogenesis network. Protein abbreviations are defined and differential expression of proteins are provided in Supplementary Table S1.</figcaption></figure>

It is worth noting that predictions for ANGPT2 (described in more detail below) or vasculogenesis were not made in the high vs. low FE proteomic dataset; cell death of muscle cells, necrosis of muscle and apoptosis were all predicted to be inhibited in the high FE phenotype. The activation state and *p* value of overlap were all significant for necrosis of muscle cells (−2.19, *p* = 1.52 × 10⁻⁵), cell death of muscle cells (−2.59, *p* = 1.54 × 10⁻⁶), and apoptosis of muscles (−2.81, *p* = 1.29 × 10⁻⁴) for the high FE phenotype, which suggests that a commonality of mechanisms were presented in the high FE PedM broiler and for broilers exhibiting WB myopathy.
Using the overlay function of the IPA program, it is possible to project one dataset onto another to see commonalities and differences between the datasets. In Figure 2B, the FE proteomic dataset (from Kong et al. [11]) was projected (overlaid) onto the WB myopathy vasculogenesis functional network. The result shows all of the proteins in WB proteomics dataset were not differentially expressed (fold difference, 1.3, *p* < 0.05) in the FE dataset, i.e., vasculogenesis was not predicted to be a significant function in the high vs. low FE PedM phenotypes.

In the necrosis functional network (Figure 3), all up- and down-regulated proteins in WB myopathy contributed to the prediction of inhibition of necrosis; there were no inconsistencies in the relationships between the proteins and the downstream effect of predicted inhibition of necrosis. With the large number of predictions that were obtained in the WB myopathy dataset, the discussion of each upstream regulator may not be sufficient in detail and critical literature citations may be missed. Nonetheless, we hope it will spark hypotheses to be tested that will provide insight into this disease that can lead to preventive measures to ameliorate significant economic losses to the poultry industry.

<figure><img src="image_81.png" /><figcaption>**Figure 3.** The network of differentially expressed proteins that was used to predict inhibition of necrosis in wooden breast (WB) myopathy. Protein abbreviations are defined and differential expression of proteins are provided in Supplementary Table S1.</figcaption></figure>

Although there was no commonality between WB myopathy and FE datasets using the overlay function in IPA, three functions—muscle cell death, apoptosis of muscle cells, and necrosis of muscle cells—were predicted to be inhibited in the high FE vs. low FE PedM dataset (Supplementary Table S2). Thus, there was a commonality between WB myopathy and high FE with respect to these functions. This will be discussed in a little more detail in the summary section towards the end of this manuscript.
3.4. Upstream Regulators

3.4.1. Angiopoeitin 2

The top predicted activated upstream regulator in the WB myopathy was ANGPT2 (Table 2). Murtyn et al. [15] indicated that hypoxia and circulatory insufficiency may play a role in the development of WB myopathy. Thus, the predicted activation of ANGPT2 by IPA could be a tissue-mediated signal to increase blood vessel development in damaged tissue to enhance oxygen delivery. This hypothesis is supported by the predicted activation of functions of development of vasculature, vasculogenesis and angiogenesis shown in Table 3. Moffarri and Hossein [16] reported that ANGPT2 expression in primary skeletal muscle myocytes was responsive to hydrogen peroxide-induced oxidative stress, but not to proinflammatory cytokines. As ANGPT2 promotes skeletal myoblast survival and repair, Moffarri and Hossein [16] hypothesized that muscle-derived ANGPT2 production plays a positive role in muscle fiber repair. Conversely, ANGPT2 activation could help in the delivery of oxygen and nutrients to aid growth and muscle development as birds with WB myopathy were heavier than birds without WB myopathy (Table 1). In line with this general connection of WB to impaired aerobic energy supply, we have recently found that heavier birds (known to be pre-disposed to WB) possess a lower skeletal muscle mitochondrial content, which would reduce the capacity for aerobic ATP synthesis [17].

The ANGPT2 network of downstream molecules presented in Figure 4A indicates that all of the downstream molecules in WB myopathy were up-regulated compared to normal muscle. Each of these molecules are indicated to be up-regulated indirectly (dashed orange line and arrow) by ANGPT2, with the exception of ribosomal protein 6 (RP6), whose up-regulation was not predicted (black dashed arrow line). This relationship, i.e., an unpredicted effect, occurs when there are too few literature citations in the IPA database to establish a clear relationship between the upstream regulator and the downstream target molecule. When the FE proteomics dataset was overlaid (projected) onto the WB myopathy ANGPT2 network (Figure 4B), there were no commonalities in protein expression and therefore no prediction could be made of ANGPT2 activity.

Using the 'Grow' function of the IPA program, it is possible to link upstream regulators to downstream functions, as shown in Figure 4C. In this figure, predictions of activation of vasculogenesis and angiogenesis in WB myopathy are the result of increased VIM, CRYAB, HSPA5 and HSP90AA1 expression and the predicted activation of ANGPT2. All of these are linked by dashed orange lines with arrows. The dashed lines indicate that at least one additional step was between the upstream molecule and the downstream target. Apoptosis and necrosis were predicted to be inhibited in WB myopathy due to increased expression of CRYAB, HSPA5, HSP90AA1, CALR and HSPA2 and by the predicted activation of ANGPT2. These relationships to necrosis and apoptosis functions are indicated by dashed blue lines with a short perpendicular line at the end indicating inhibition. The lighter blue color for necrosis indicates that the prediction made for this function, based on proteins in this network, was qualitatively less than for apoptosis. The darker color (stronger qualitative prediction) for apoptosis is due in part to more proteins in the regulatory network (five up-regulated proteins plus ANGT2) compared to four proteins in the necrosis network.
<figure><img src="image_82.png" /><figcaption>**Figure 4.** The upstream regulatory network of proteins displayed in hierarchical format for Angiopoeitin 2 (ANGPT2) used in the prediction of activation of this molecule (activation Z score and *p* value of overlap as shown in (A)). Overlay of the feed efficiency (FE) proteomics dataset from Kong et al. [11] reveals no commonality of differentially expressed proteins between the wooden breast (WB) and FE data (B). A regulatory network for ANGPT2 indicates target proteins in the dataset that contributed to the predictions of active vasculogenesis and angiogenesis, and inhibited necrosis and apoptosis is provided in (C). Protein abbreviations are defined and differential expression of proteins are provided in Supplementary Table S1.</figcaption></figure>
### 3.4.2. Epidermal Growth Factor Receptor

The upstream regulator network of epidermal growth factor receptor (EGFR) in WB myopathy is presented in Figure 5. EGFR was predicted to be activated in WB myopathy, but the overlay of the FE proteomics dataset resulted in no prediction, indicating no commonality with this upstream regulatory network between WB myopathy and high FE. There are thousands of literature citations for EGFR involvement in carcinoma, but here we will focus on a few citations that are pertinent to EGFR in skeletal muscle development. Olwin and Hauschka [18] reported that EGFR and fibroblast growth factor receptor were permanently lost in terminal differentiation of adult mouse skeletal muscle in vitro. It was also reported that whereas EGF increases in skeletal muscle as pigs age, mRNA expression of the receptor declines [19]. Finally, Leroy et al. [20] reported that down-regulation of EGFR triggers differentiation of human myoblasts in cell culture, indicating that EGFR expression normally declines with age. Therefore, the prediction of activation of EGFR in WB myopathy might indicate an abnormal process that contributes to the pathophysiology in WB myopathy.

<figure><img src="image_83.png" /><figcaption>**Figure 5.** The regulatory network of epithelial growth factor receptor (EGFR) showing downstream target proteins contributing to the prediction of activation of vasculogenesis and inhibition of necrosis and apoptosis in wooden breast (WB) myopathy. Protein abbreviations are defined and differential expression of proteins are provided in Supplementary Table S1.</figcaption></figure>

### 3.4.3. Transforming Growth Factor Beta 1 (TGFB1)

TGFB1 was predicted to be activated in WB myopathy proteomics (Table 2, Figure 6A) as well as in the high FE PedM phenotype (Figure 6B). While the TGFB1 networks in the two studies did not share any differentially expressed proteins, the overlay of the FE proteomics data on the WB dataset (Figure 6C) as well as the overlay of the WB proteomics data on the FE dataset (Figure 6D) resulted in predictions of the activation of TGFB1. TGFB1 is a cytokine that is involved in controlling the growth, proliferation and differentiation of cells and able to serve in an autocrine manner to regulate its own expression [21]. This autocrine stimulation is indicated by the semi-circular arrow in each of the networks in Figure 6. In the study conducted by Hubert et al. [22], TGFB1 was also identified as one of the upstream regulators of the WB myopathy. The role of TGFB1 in this myopathy may be explained in part by its effects on cell differentiation and mitochondrial function. Under normal conditions, this cytokine is involved in the regulation of cell differentiation,
increased mitochondrial activity and increased oxidative phosphorylation. However, studies have shown that increased TGFB1 expression is detrimental and can produce mitochondrial dysfunction with increased oxygen radical production that compromises the mitochondrial antioxidant system [23–25]. The resulting mitochondrial perturbations have been illustrated to lead to the inhibition of myogenic differentiation and the formation of tissue fibrosis [26,27]. The potential for diminished mitochondrial biogenesis in WB myopathy is further supported by predicted inhibition of PGC1α, which is discussed in more detail below. Feed efficiency studies in cattle have also identified TGFB1 as being an upstream regulator [28,29]. This provides further evidence of the of TGFB1 involvement in genetic differences in the various feed efficiency phenotypes.

<figure><img src="image_84.png" /><figcaption>**Figure 6.** The upstream regulatory network of transforming growth factor beta 1 (TGFB1) for wooden breast (WB) (A) and feed efficiency (FE) (B) proteomic datasets. Downstream target proteins that were differentially expressed are shown as up-regulated in pink or red or down-regulated in WB myopathy and high FE PedM broilers. Overlay of the FE proteomic dataset on the WB myopathy network (C) and of the FE data on the WB network (D) resulted in predicted activations of TGFB1 as well. Therefore, despite dissimilar proteomic landscapes of the two datasets, TGFB1 is predicted to be active in both WB myopathy and high FE.</figcaption></figure>
Protein abbreviations for WB myopathy data (A, C) are defined and differential expression of proteins are provided in Supplementary Table S1. Protein abbreviations for FE proteomics (B and D) are as follows: ACTA1 (actin, alpha 1, skeletal muscle), ACTC1 (actin, alpha, cardiac muscle 1), ALB (albumin), ARF (ADP-ribosylation factor 1), ARF4 (ADP-ribosylation factor 4), C4BPA (complement component 4 binding protein, alpha), CAV1 (caveolin 1), CCT2 (cell division cycle 42), CKM (creatinine kinase, muscle), CSPG4 (chondroitin sulfate proteoglycan 4), CTSC (cathepsin C), EEF1A (eukaryotic translation elongation factor 1 alpha 1), ENO2 (enolase 2, gamma neuronal), GPI (glucose 6 phosphate isomerase), IRS1 (insulin receptor S1), OARD1 (O-acyl-ADP-ribose deacylase), PRPS1 (phosphoribosyl pyrophosphate synthetase 1), PSMD1 (proteosome 26S subunit, non-ATPase, 1), SLC25A4 (solute carrier family 25 (mitochondrial carrier; adenine nucleotide translocator), member 4), VDAC2 (voltage-dependent anion channel 2).

Linkage of TGFB1 activation to downstream functions through the differentially expressed proteins is provided in Figure 7. For simplicity, only relationships that clearly predicted activation (angiogenesis, vasculogenesis) or inhibition (apoptosis, necrosis) are presented in this figure. The expression of eight proteins in this regulatory network (PGK1, CTSB, VIM, HSP90AA1, FLNA, HSPA5, ANAXA2 and RHOA) were involved with predictions of increased angiogenesis and vasculogenesis in WB myopathy. The up-regulation of all the proteins, with the exception of PGK1, which was down-regulated, were downstream targets of TGFB1 that contributed to the predictions of enhanced blood vessel-capillary formation. Eight up-regulated proteins (HSP90AAA1, FLNA, HSPA5, ANXA2, RHOA, RACK1, NFRNPH1 and DES) and one down-regulated protein (HINT1) were involved in the predictions of reduced apoptosis and necrosis in wooden breast myopathy.

<figure><img src="image_85.png" /><figcaption>**Figure 7.** The regulatory network of TGFB1 in wooden breast (WB) myopathy leading to predicted activations of angiogenesis and vasculogenesis and inhibition of apoptosis and necrosis through downstream target proteins in WB myopathy. Protein abbreviations are defined and differential expression of proteins are provided in Supplementary Table S1.</figcaption></figure>

### 3.4.4. Nuclear Factor (Erythroid-Derived 2)-Like 2 (NFE2L2)

Nuclear factor (erythroid-derived 2)-like 2 (NFE2L2) coordinates cellular response to oxidative stress [30,31]. NFE2L2 was predicted to be activated in both WB proteomics and in the FE proteomics dataset in PedM broilers (Figure 8A,B), but the predictions were based on different downstream molecules. Overlaying of each dataset on the other upstream regulatory network (Figure 6C,D) resulted in predicted activation of NFE2L2.
A

<figure><img src="image_86.png" /></figure>

WB Proteomics
Z Score = 2.18, *p* = 4.3 × 10⁻⁹

C

<figure><img src="image_87.png" /></figure>

FE Proteomics Overlay of
WB Proteomics NFE2L2 Network

B

<figure><img src="image_88.png" /></figure>

FE Proteomics
Z Score = 2.38, *p* = 1.54 × 10⁻⁶

D

<figure><img src="image_89.png" /></figure>

WB Proteomics Overlay of FE
NFE2L2 Network

**Figure 8.** The upstream regulator network for NFE2L2 (nuclear factor (erythroid-derived 2)-like 2) for wooden breast (WB) myopathy (activation Z score = 2.18, 4.3 × 10⁻⁹) compared to normal muscle (A) and for the high vs. low PedM broiler feed efficiency (FE) proteomic data (activation Z score = 2.38, *p* = 1.54 × 10⁻⁶) (B). Overlay of the FE data on the WB myopathy data (C) and vice versa (D) both resulted in predicted activations of NFE2L2. Protein abbreviations in B and D in the FE data set are as follows: ARF1 (ADP-ribosylation factor 1), COX4I1 (cytochrome oxidase subunit IV isoform 1), EIF4G2 (eukaryotic translation initiation factor 4 gamma 2), GPX1 (glutathione peroxidase 1), PPIB (peptidylprolyl isomerase B (cyclophilin B), PSMB1 (proteosome 26S subunit beta type, 1) PSMD1 (proteosome 26S subunit non-ATPase, 1), RYR3 (ryanodine receptor 3), SOD1 (superoxide dismutase 1), UGDH (UDP-glucose 6-dehydrogenase).

Under normal (non-oxidant stress) conditions, the NFE2L2 protein remains bound to Cullin3 (CUL3) and Kelch like-ECH protein 1 (KEAP1), and is then rapidly directed (within
15-20 min) to proteasomes for degradation [32]. However, in response to oxidative stress, NFE2L2 is transported to the nucleus, where it stimulates antioxidant gene expression after binding to the antioxidant response element. Recently, NFE2L2 was predicted to be activated in high feed efficiency (HFE) compared to low feed efficiency (LFE) broiler phenotypes based on downstream target molecules [11,33]. It has been shown that NFE2L2 undergoes post-transcriptional and translational modifications that can affect the amount of NFE2L2 protein produced, which, along with rapid activation or degradation, makes it an extremely labile protein [32,34]. Avian NFE2L2 appears to function the same as mammalian NFE2L2, but since it shares only 67% homology with mammalian NFE2L2 [35], it may have other roles in avian tissues besides coordination of antioxidant response. In a recent review, NFE2L2 was observed to affect mitochondrial function in many ways, including efficiency of oxidative phosphorylation, mitochondrial biogenesis and mitochondrial integrity [36].

### 3.4.5. X-Box Binding Protein 1

X-box binding protein 1 (XBP1) was predicted to be activated in WB myopathy (Table 2). There was no commonality of protein expression of this upstream regulator in the FE proteomic dataset. Myoblasts that were forced to differentiate in response to XBP1 overexpression produced myotubes that were shorter and less mature indicative of impaired differentiation relative to control cells [37]. Of particular potential relevance to WB myopathy is the following from this study: "Our observations suggest that skeletal muscle tissue is particularly sensitive to physiological stressors that trigger the unfolded protein response, namely glucose deprivation, anabolic stimulation, and perhaps other endoplasmic reticulum stresses, such as hypoxia, imbalances in calcium homeostasis, and ischemia. These perturbations could be triggered by normal or pathological states." The presence of hypoxia in WB myopathy has been reported [38] and will be discussed in greater detail below.

### 3.4.6. Rapamycin-Insensitive Companion of Mammalian Target of Rapamycin (RICTOR)

RICTOR plays an important role in actin-cytoskeletal development, as formation of this structure was impaired in a RICTOR knockout mouse model [39–41]. Ablation of RICTOR (mTORC2) did not adversely affect muscle function in mice, results that contrast dramatically with the loss of function and development of muscle dystrophy when RAPTOR (mTORC1) was removed [42].

RICTOR was predicted to be inhibited in WB myopathy with activation z scores of -3.61 (Table 2). RICTOR was also predicted to be inhibited in PedM broilers exhibiting high FE [11]. The upstream regulator networks of RICTOR in the current WB myopathy study and in the FE proteomics study are presented in Figure 9A,B, respectively. Although both studies yielded predictions of RICTOR inhibition, it is apparent that the predictions were based on different downstream target molecules in the respective datasets. The overlay of the datasets on each other yielded no prediction giving indication of no commonality in the RICTOR networks. Thus, this appears to indicate that the predictions of inhibition of RICTOR in WB myopathy and FE breast muscle proteomics are due to separate mechanisms. The connection between inhibition of RICTOR to inhibition of necrosis shown in Figure 9C is shown to be mediated by up-regulation of several ribosomal proteins. Apoptosis was predicted to be weakly activated in this network through combined inhibition of RICTOR and up-regulation of RPS3.
A

<figure><img src="image_90.png" /></figure>

B

<figure><img src="image_91.png" /></figure>

C

<figure><img src="image_92.png" /></figure>

**Figure 9.** The upstream regulatory network for RICTOR in wooden breast (WB) myopathy (A) and for feed efficiency (FE) proteomics (B). All downstream proteins were up-regulated in both datasets, but there were no common differentially expressed proteins. The regulatory network for RICTOR linking downstream molecules to the predicted inhibition of necrosis and slight activation of apoptosis in WB myopathy (C). For (A,C), protein abbreviations are defined and differential expression of proteins in WB myopathy are provided in Supplementary Table S1. Proteins in (B) are as follows: COX4I1 (cytochrome oxidase subunit IV isoform 1), NDUFB5, B8, S7 (NADH dehydrogenase (ubiquinone) Complex I assembly factor 5, factor 8, subunit 7), PPA2 (pyrophosphate (inorganic) 2, PSMB1, D1, D2 (proteosome 26S subunit, non-ATPase beta type-1) PSMD1, PSMD2 (proteosome 26S subunit, non-ATPase, 1 & non-ATPase 2), RPL12, RPL30 (ribosomal protein L12 & L30), UQCRC2 (ubiquinol-cytochrome c-reductase core protein II).

3.4.7. Peroxisome Proliferator-Activating Gamma Coactivator 1 Alpha (PPARGC1α or PGC1α)

Because PGC1α has been described as the master regulator of mitochondrial biogenesis (Nisoli et al. [43,44]), the predicted inhibition of PGC1α by upstream regulator analysis (Figure 10A, activation Z score = -2.35, *p* = 1.92 × 10⁻⁵ suggests that mitochondrial biogenesis, and possibly mitochondrial function, would be impaired in WB myopathy. In contrast, the upstream regulator analysis of PGC1α was predicted to be activated in the high FE PedM broiler phenotype (Figure 10B, activation Z score = 2.45, *p* = 1.84 × 10⁻⁴). There was no commonality of protein expression between the WB and FE datasets. The mitochondrial canonical pathways in WB and FE proteomics datasets are shown in Figure 10C,D, respec-
tively. Whereas Complex I of the mitochondrial electron transport chain was predicted to be inhibited in WB myopathy (Figure 10C), Complex I, III and IV were predicted to be activated based on expression of proteins in the FE dataset (Figure 10D). Complex I (NADH dehydrogenase), which accepts NADH-linked energy substrates, is the largest complex in the electron transport chain, consisting of approximately 40 proteins. Thus, inhibition of this complex could have a dramatic influence on the energy production in WB tissues. The potential for diminished mitochondrial function in WB myopathy is supported by a recent report showing increased mitochondrial damage and mitophagy in histological analysis of WB myopathy [45]. Several studies provide evidence of enhanced mitochondrial function and complex activities in the high FE PedM phenotype (e.g., see review by [46]) that may be attributed to increased activity of PGC1α in the high FE phenotype shown in Figure 10B.

<figure><img src="image_93.png" /><figcaption>**Figure 10.** The upstream regulator network for PPARGC1α for WB myopathy vs. control breast muscle in commercial broilers (A) and high vs. low FE PedM (B) proteomic datasets [11]. The predicted inhibition of PGC1α in WB myopathy was associated with the predicted inhibition of Complex I of the mitochondrial electron transport chain in WB myopathy (C), whereas the predicted activation of PGC1α in high FE was associated with predicted activations of Complex I, III and IV of the mitochondrial electron transport chain (D). The regulatory network for PPARGC1A showing the downstream target molecules involved in the predictions of inhibition of vasculogenesis and oxidation of fatty acid and activation of reactive oxygen species (ROS) production and quantity of adipose tissue (E). Protein abbreviations in (A,E) are defined and differential expression of proteins are provided in Supplementary Table S1. Protein abbreviations for (D) are as follows: COX4I1 (cytochrome oxidase subunit IV isoform 1), CTSC (cathepsin C), GPX1 (glutathione peroxidase 1), NDUFB5 (NADH dehydrogenase (ubiquinone) Complex I assembly factor 5), NDUFB8 (NADH dehydrogenase (ubiquinone) beta subcomplex, 8), OMM (Outer Mitochondrial Membrane), IMM (Inner Mitochondrial Membrane).</figcaption></figure>
The regulatory network for PGC1α presented in Figure 10E indicates that, unlike the overall prediction of vasculogenesis activation (Table 3), vasculogenesis was predicted to be inhibited by a combination of down-regulation of LDHA and predicted inhibition of PGC1α. Reactive oxygen species (ROS) production was predicted to be enhanced in WB myopathy due to a combination of LDHA down-regulation and up-regulation of C3. Enhanced ROS production may contribute to tissue damage reported in WB myopathy. In this network, oxidation of fatty acid was predicted to be inhibited in WB myopathy in the present study due to the up-regulation of C3, down-regulation of GOT2 and predicted inhibition of PGC1α and the quantity of adipose tissue was predicted to be activated in WB myopathy. These findings agree with many studies, indicating that lipid accumulation in muscle is a characteristic in WB myopathy [6,47,48].

### 3.4.8. Aryl Hydrocarbon Receptor Nuclear Translocator (ARNT) and Aryl Hydrocarbon Receptor Nuclear (AHR)

ARNT and AHR were predicted to be inhibited in WB myopathy (Table 2). The regulatory networks shown in Figure 11A,B, respectively, for ARNT indicate that fibrosis would be predicted to be enhanced, and glycolysis, vasculogenesis as well as apoptosis would be predicted to be inhibited based on expression of proteins within this regulatory network. In contrast, angiogenesis and necrosis were predicted to be activated in the AHR regulatory network.

<figure><img src="image_94.png" /></figure>

<figure><img src="image_95.png" /></figure>

**Figure 11.** The upstream regulatory networks for ARNT (A) and AHR (B). The predicted inhibition of ARNT was based on the downstream expression of proteins that are shown and was associated with predicted inhibition of apoptosis, vasculogenesis and glycolysis and prediction of activation of fibrosis. Conversely, the expression of downstream target molecules of involved in predicting inhibition of AHR were associated with predictions of enhanced angiogenesis and necrosis. The discrepancy in the prediction of active necrosis in the AHR and prediction of inhibition of necrosis in the entire dataset shown in Table 2 is based on the smaller subset of expression data for AHR compared to the overall dataset. Protein abbreviations are defined and differential expression of proteins are provided in Supplementary Table S1.

Murtyn et al. [15] provided a sound argument for the presence of hypoxic conditions playing a role in WB myopathy that was based on the differential expression of several downstream molecules that are responsive to hypoxia inducible factor 1 (HIF1 and HIF-dependent expression). However, they also indicated that transcripts of the two subunits of HIF-1, HIF-1α and ARNT were not differentially expressed between normal and WB myopathy breast muscle. Since both AHR and ARNT are involved in other cellular processes (e.g., xenobiotic metabolism), the predictions of inhibition in the present study
(and lack of difference in gene expression in Murtyn et al. [15]) may have nothing to do with tissue hypoxia, but rather, are playing roles in other cellular processes. For example, AHR knockout in aortic smooth muscle cells resulted in the deregulation of components of TGF-β signaling and could directly affect the metabolism of toxic substances [49]. Recent reviews indicated that there are numerous functions that AHR may play in the cell [50] and immunity [51]. Greene et al. [38] measured lower oxygen levels in breast muscle tissue exhibiting WB myopathy compared to normal tissue. Furthermore, several genes associated with hemoglobin in red blood cells and HIF-1α in breast muscle were down-regulated in WB compared to normal breast muscle tissue [38].

Some discussion is warranted with regard to the predicted inhibition of AHR in the present study. First, predictions of activation or inhibition of an upstream molecule generated in the IPA program are based on the expression of downstream molecules specifically associated with the dataset being examined. In the case of the WB proteomic dataset in the present study, there were far fewer differentially expressed proteins (~130) that were detected compared to approximately 1500 differentially expressed transcripts detected by Murtyn et al. (2015) [15]. Secondly, denaturing gels used in separating proteins prior to protein expression analysis can result in a loss of many proteins with large hydrophobic regions (e.g., membrane-bound or membrane associated proteins), which in turn could influence the predictions generated by the IPA program. Third, the lack of differential expression of ARNT and HIF-1α m-RNA reported by Murtyn et al. [15] does not mean that the proteins would also not be differentially expressed. Furthermore, these proteins may also undergo post-translational modifications that would be required for activation of downstream target molecules. The predicted activation of ANGPT2 in WB myopathy could be a response to tissue hypoxia to increase vascularization caused by tissue damage and/or hypoxia. Thus, further examination AHR and ARNT signaling in WB myopathy is warranted.

Although no predictions were made regarding AHR and ARNT in the FE proteomics study by Kong et al. [11], ESR1 (estrogen receptor 1) was predicted to be activated (activation Z score = 1.97, 3.93 × 10⁻³). In the IPA gene view in the Qiagen software, ESR1 is indicated to be a member of AHR-aryl hydrocarbon-Arnt-ESR1. Thus, this would appear a major difference between the WB and FE proteomics datasets.

### 3.4.9. Estrogen Regulator Receptor Gamma (ERSSγ)

ERSSγ was predicted to be inhibited in WB myopathy (Table 2). The regulatory network of ERSSγ shown in Figure 12 indicates that vasculogenesis would be inhibited by the down-regulation of LDHA and PKM, whereas down-regulation of GAPDH caused a prediction of enhanced necroptosis in WB myopathy. Again, both of these predictions are opposite to functions outlined in Table 3, with the discrepancy attributed to the larger number of proteins involved in the overall predictions shown in Table 3.

Fan et al. [52] reported that ESRRγ helps in reducing muscle damage and improves muscle function in a double PGC1α/β knockout mouse model. The loss of PGC1α/β resulted in a decrease in expression of myoglobin and a lighter color of muscle consistent with a change from oxidative to glycolytic fibers, which were largely restored when ESRRγ transgenic mice were crossed with the PGC1α/β knockout mice [52]. Defects in mitochondrial function from the PGC1α/β knockout were also restored by crossing with the ESRRγ transgenic mice. Furthermore, ESRRγ also increased the expression of genes associated with angiogenesis (e.g., vascular endothelial growth factor A and fibroblast growth factor 1) and vascular density (CD31).
<figure><img src="image_96.png" /><figcaption>**Figure 12.** ESRRGamma (ESRRG or ESSRγ) was predicted to be inhibited in WB myopathy based on the expression of downstream molecules shown in this regulatory network. The regulatory network resulted in the prediction of the inhibition of vasculogenesis and activation of necroptosis. As in Figure 9, discrepancies in the predictions of these functions and those in Table 2 are due to a smaller subset of proteins in the ERRG upstream regulatory network, compared to the larger set of proteins listed for functions in Table 2. Protein abbreviations are defined and differential expression of proteins are provided in Supplementary Table S1.</figcaption></figure>

### 3.5. Summary and Synopsis

There were two major goals of this study; (1) to present an upstream regulatory factor analysis associated with wooden breast myopathy to expand our understanding of this disease and (2) to determine where contributions to WB myopathy might be linked to selection for performance (high FE) in PedM broilers. Table 4 summarizes the results of comparison of upstream regulator analysis between the WB proteomic data to the FE proteomic data reported by [11]. With the exception of TGFB1 and NFE2L2, there were no commonalities for upstream regulators between the WB and FE proteomics data. While PPARGC1α was predicted to be inhibited in WB myopathy, this upstream regulator was predicted to be activated (based on a completely different set of differentially expressed proteins) in muscle of the high FE compared to the low FE PedM broiler. From these findings, we hypothesize that mitochondrial biogenesis would be inhibited in WB myopathy and activated in the high FE PedM broiler [11].

<table><caption>**Table 4.** Upstream regulators that were predicted to be activated (+) or inhibited (−) in wooden breast myopathy (WB) proteomics (from Kuttappan et al. [2]) with the resulting prediction when the FE proteomics dataset from high and low feed efficiency (FE) Pedigree Broilers males [11] was overlaid on the WB upstream regulator network. The designation ‘no’ indicates that no prediction was made with the FE overlay (projection of FE dataset onto WB myopathy dataset).</caption><thead><tr><th>Upstream Regulator</th><th>WB</th><th>FE Overlay</th><th>Upstream Regulator</th><th>WB</th><th>FE Overlay</th></tr></thead><tbody><tr><td>ANGPT2</td><td>+</td><td>no</td><td>RICTOR</td><td>-</td><td>no</td></tr><tr><td>EGFR</td><td>+</td><td>no</td><td>PPARGC1α</td><td>-</td><td>no</td></tr><tr><td>XBP1</td><td>+</td><td>no</td><td>AHR</td><td>-</td><td>no</td></tr><tr><td>TGFB1</td><td>+</td><td>+</td><td>ARNT</td><td>-</td><td>no</td></tr><tr><td>NFE2L2</td><td>+</td><td>+</td><td>ESRRG</td><td>-</td><td>no</td></tr></tbody></table>

Figure 13 provides a graphical presentation of the upstream regulators that were predicted to be activated (orange) or inhibited (blue) in WB myopathy (13A) or FE proteomics (13B) and predicted interactions with other upstream regulators described in this study.
This figure was generated by uploading the upstream regulator networks in Table 2, and then removing differentially expressed downstream target proteins to reveal predicted interactions between the upstream regulators. Major functions/processes, including the prediction of the activation of angiogenesis and the inhibition of cell death (apoptosis and necrosis), were discussed previously. The predicted activation of TGFB1 is indicated as having an indirect (through one or more mediators) inhibitory effect on PPARGC1A and predicted inhibition of PPARGC1A would directly inhibit EGFR. Predictions of reduced activity of AHR and ARNT both contribute to the prediction of lower ESRRG activity. Similarly, predicted inhibition of ARNT and AHR would inhibit ESRRG as well. Because PPARGC1A and ESRRG enhance mitochondrial biogenesis, their predicted inhibition could be hypothesized to lower mitochondrial biogenesis. Elevated NFE2L2 activity is predicted to result in enhanced activation of TGFB1. The predicted inhibition of AHR would function to enhance TGFB1 activity, since AHR has been shown to inhibit TGFB1. NFE2L2 activation would contribute to the predicted activation of TGFB1. Effects of ANGPT2 and TGFB1 would enhance EGFR and their combination would lead to the prediction of enhanced angiogenesis in WB myopathy.

<figure><img src="image_97.png" /><figcaption>**Figure 13.** (A) Upstream regulator interactions contributing to the overall prediction of increased vasculogenesis and decreased cell death in wooden breast muscle myopathy in commercial broilers. The combination of the predicted activities of these upstream regulators would support the prediction of increased angiogenesis and inhibition of cell death in wooden breast myopathy. (B) Arrangement of similar upstream regulators in the high vs. low feed efficiency PedM phenotype dataset revealed a scenario in which angiogenesis was not predicted to be active or inhibited while predictions were made of inhibition cell death, which is similar to that for WB myopathy.</figcaption></figure>
Figure 13B presents upstream regulators from the FE proteomics data. No prediction was made concerning ANGPT2, and unlike WB myopathy, EGFR was predicted to be inhibited in the high FE phenotype. No prediction was made in high vs. low FE data concerning angiogenesis/vasculogenesis. TGFB1 was predicted to be activated in high FE, as were NFE2L2 and PPARGC1A, but as indicated in the text, these were based on different downstream target molecules compared to the WB myopathy data. No predictions were made for AHR, ESSRG or ARNT, but estrogen receptor 1 (ESR1), a member or the ARN family, was predicted to be activated in the high FE phenotype. The combination of increased PPARGC1A and ESR1 would act to enhanced mitochondrial biogenesis in the high FE phenotype.

All of the predictions presented in this study are just that—predictions. Each of these predictions of upstream regulator activity represents hypotheses and mechanistic studies that need to be conducted in order to determine if these upstream regulators contribute to the development of WB muscle myopathy. We hope that they will help shed light on the understanding of WB myopathy that will lead to methods, whether nutritional or chemical, that can be used to ameliorate WB muscle myopathy in commercial broilers.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/2304-8158/10/1/104/s1, Supplementary Table S1. Abbreviations and names with P value and fold change of proteins presented in Tables 2 and 3 listed alphabetically. Fold Differences in green (negative values) were down-regulated in wooden breast (WB) myopathy whereas fold differences in red (positive values) were up-regulated in WB myopathy breast muscle relative to values expressed in normal breast muscle tissue. Supplementary Table S2. Functions predicted to be inhibited (blue) in breast muscle of Pedigree Male Broilers exhibiting a high compared to low feed efficiency phenotype (previously unpublished data from Kong et al. [11]).

**Author Contributions:** Design of studies: W.G.B., V.A.K., C.M.O. and B.C.K.; bioinformatics: V.A.K. and B.C.K.; interpretation of data—pathway analysis: W.G.B., B.C.K., N.J.H. and V.A.K.; drafting and critical review of manuscript: W.G.B., K.R.L., N.J.H., V.A.K., B.A., S.D., C.M.O. and B.C.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was provided by the USDA-NIFA (#2013-01953) to W.B., B.K. and N.J.H., USDA-NIFA SAS (#2019 69012-29905) to W.B., S.D., and B.K. and the Arkansas Biosciences Institute to W.B., B.K. and S.D.

**Institutional Review Board Statement:** The present study was conducted in accordance with the recommendations in the guide for the care and use of laboratory animals of the National Institutes of Health. All procedures for animal care complied with the University of Arkansas Institutional Animal Care and Use Committee (IACUC): Protocol Nos. 14,012 and 17,080.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** No new data were created in this study. Data sharing is not applicable to this article. Support for this research was provided by the USDA-NIFA (#2013-01953) to W.B., B.K. and N.H., USDA-NIFA SAS (#2019 69012-29905) to W.B., S.D., B.K. and N.H. and the Arkansas.

**Acknowledgments:** Support was provided by USDA-NIFA (#2013-01953), Chancellors Challenge Fund (University of Arkansas), USDA-NIFA Sustainable Agriculture Systems (#2019-69012-29905), and Division of Agriculture, University of Arkansas.

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

## References

1. Velleman, S.G. Relationship of Skeletal Muscle Development and Growth to Breast Muscle Myopathies: A Review. *Avian Dis.* **2015**, *59*, 525–531. [CrossRef]

2. Kuttappan, V.A.; Bottje, W.; Ramnathan, R.; Hartson, S.D.; Coon, C.N.; Kong, B.-W.; Owens, C.M.; Vazquez-Añon, M.; Hargis, B.M. Proteomic analysis reveals changes in carbohydrate and protein metabolism associated with broiler breast myopathy. *Poult. Sci.* **2017**, *96*, 2992–2999. [CrossRef] [PubMed]
3. Lake, J.A.; Abasht, B. Glucolipotoxicity: A Proposed Etiology for Wooden Breast and Related Myopathies in Commercial Broiler Chickens. *Front. Physiol.* **2020**, *11*, 1–13. [CrossRef] [PubMed]

4. Velleman, S.G.; Clark, K.L. Histopatholoic and myogenic gene expression changes associated with wooden breast in broiler breast muscle. *Avian Dis.* **2015**, *59*, 410–418. [CrossRef] [PubMed]

5. Abasht, B.; Mutryn, M.F.; Mickalek, R.D.; Lee, W.R. Oxidative stress and metabolic perturbations in wooden breast disorder in chickens. *PLoS ONE* **2016**, *11*, e0153750. [CrossRef]

6. Papah, M.B.; Brannick, E.M.; Schmidt, C.J.; Abasht, B. Evidence and role of phlebitis and lipid infiltration in the onset and pathogenesis of Wooden Breast Disease in modern broiler chickens. *Avian Pathol.* **2017**, *46*, 623–643. [CrossRef]

7. Sihvo, H.-K.; Lindén, J.; Airas, N.; Immonen, K.; Valaja, J.; Puolanne, E. Wooden Breast Myodegeneration of Pectoralis Major Muscle over the Growth Period in Broilers. *Vet. Pathol.* **2017**, *54*, 119–128. [CrossRef]

8. Brothers, B.; Zhuo, Z.; Papah, M.B.; Abasht, B. RNA-Seq Analysis Reveals Spatial and Sex Differences in Pectoralis Major Muscle of Broiler Chickens Contributing to Difference in Susceptibility to Wooden Breast Disease. *Front. Physiol.* **2019**, *10*, 764. [CrossRef]

9. Petracci, M.; Soglia, F.; Madruga, M.; Carvalho, L.; Ida, E.; Estévez, M. Wooden-Breast, White Striping, and Spaghetti Meat: Causes, Consequences and Consumer Perception of Emerging Broiler Meat Abnormalities. *Compr. Rev. Food Sci. Food Saf.* **2019**, *18*, 565–583. [CrossRef]

10. Gratta, F.; Fasolato, L.; Birolo, M.; Zomeño, C.; Novelli, E.; Massimiliano, P.; Pascual, A.; Xiccato, G.; Trocino, A. Effect of breast myopathies on quality and microbial shelf life of broiler meat. *Poult. Sci.* **2019**, *98*, 2641–2651. [CrossRef]

11. Kong, B.; Lassiter, K.; Piekarski-Welsher, A.; Dridi, S.; Reverter-Gomez, A.; Hudson, N.J.; Bottje, W.G. Proteomics of breast muscle tissue associated with the phenotypic expression of feed efficiency within a pedigree male broiler line: I. Highlight on mitochondria. *PLoS ONE* **2016**, *11*, e0155679. [CrossRef]

12. Tijare, V.; Yang, V.; Kuttappan, V.; Alvarado, C.; Coon, C.; Owens, C. Meat quality of broiler breast fillets with white string and woody breast muscle myopathies. *Poult. Sci.* **2016**, *95*, 2167–2173. [CrossRef] [PubMed]

13. Bottje, W.; Brand, M.D.; Ojano-Dirain, C.; Lassiter, K.; Toyomizu, M.; Wing, T. Mitochondrial proton leak kinetics and relationship with feed efficiency within a single genetic line of male broilers. *Poult. Sci.* **2009**, *88*, 1683–1693. [CrossRef] [PubMed]

14. Voruganti, S.; Lacroix, J.C.; Rogers, C.N.; Rogers, J.; Matts, R.L.; Hartson, S. The Anticancer Drug AUY922 Generates a Proteomics Fingerprint That Is Highly Conserved among Structurally Diverse Hsp90 Inhibitors. *J. Proteome Res.* **2013**, *12*, 3697–3706. [CrossRef] [PubMed]

15. Mutryn, M.F.; Brannick, E.M.; Fu, W.; Lee, W.R.; Abasht, B. Characterization of a novel chicken muscle disorder through differential gene expression and pathway analysis using RNA-sequencing. *BMC Genom.* **2015**, *16*, 1–19. [CrossRef] [PubMed]

16. Mofarrahi, M.; Husain, S.N.A. Expression and Functional Roles of Angiopoietin-2 in Skeletal Muscles. *PLoS ONE* **2011**, *6*, e22882. [CrossRef]

17. Reverter, A.; Okimoto, R.; Sapp, R.; Bottje, W.G.; Hawken, R.; Hudson, N.J. Chicken muscle mitochondrial content appears co-ordinately regulated and is associated with performance phenotypes. *Biol. Open* **2017**, *6*, 50–58. [CrossRef]

18. Olwin, B.B.; Haushka, S.D. Cell surface fibroblast growth factor and epidermal growth factor receptors are permanently lost during skeletal muscle terminal differentiation in culture. *J. Cell Biol.* **1988**, *107*, 761–769. [CrossRef]

19. Peng, M.; Palin, M.-F.; Véronneau, S.; Lebel, D.; Pelletier, G. Ontogeny of epidermal growth factor (EGF), EGF receptor (EGFR) and basic fibroblast growth factor (bFGF) mRNA levels in pancreas, liver, kidney, and skeletal muscle of pig. *Domest. Anim. Endocrinol.* **1997**, *14*, 286–294. [CrossRef]

20. Leroy, M.C.; Perroud, J.; Darbellay, B.; Bernheim, L.; Konig, S. Epidermal Growth Factor Receptor Down-Regulation Triggers Human Myoblast Differentiation. *PLoS ONE* **2013**, *8*, e71770. [CrossRef]

21. Van Obberghen-Schilling, E.; Roche, N.S.; Flanders, K.C.; Sporn, M.B.; Roberts, A.B. Transforming growth factor B1 positively regulates its own expression in normal and transformed cells. *J. Biol. Chem.* **1988**, *263*, 7741–7746. [CrossRef]

22. Hubert, S.M.; Williams, T.J.; Athrey, G. Insights into the molecular basis of wooden breast based on comparative analysis of fast-and slow-growth broilers. *bioRxiv* **2018**, 356683. [CrossRef]

23. Casalena, G.; Daehn, I.; Bottinger, E.P. Transforming Growth Factor-β, Bioenergetics, and Mitochondria in Renal Disease. *Semin. Nephrol.* **2012**, *32*, 295–303. [CrossRef] [PubMed]

24. Abe, Y.; Sakairi, T.; Beeson, C.; Kopp, J.B. TGF-beta1 stimulates mitochondrial oxidative phosphorylation and generation of reactive oxygen species in cultured mouse podocytes, mediated in part by the mTOR pathway. *Am. J. Physiol. Ren. Physiol.* **2013**, *305*, F1477–F1490. [CrossRef] [PubMed]

25. Jain, M.; Rivera, S.; Monclus, E.A.; Synenki, L.; Zirk, A.; Eisenbart, J.; Feghali-Bostwick, C.; Mutlu, G.M.; Budinger, G.R.; Chandel, N.S. Mitochondrial reactive oxygen species regulate transforming growth factor-beta signaling. *J. Biol. Chem.* **2013**, *288*, 770–777. [CrossRef]

26. Liu, D.; Black, B.L.; Derynck, R. TGF-β inhibits muscle differentiation through functional repression of myogenic transcription factors by Smad3. *Genes Dev.* **2001**, *15*, 2950–2966. [CrossRef]

27. Negmadjanov, U.; Godic, Z.; Rizvi, F.; Emelyanova, L.; Ross, G.; Richards, J.; Holmuhamedov, E.; Jahangir, A. TGF-β1-Mediated Differentiation of Fibroblasts Is Associated with Increased Mitochondrial Content and Cellular Respiration. *PLoS ONE* **2015**, *10*, e0123046. [CrossRef]
28. Widmann, P.; Reverter, A.; Weikard, R.; Suhre, K.; Hammon, H.M.; Albrecht, E.; Kühn, C. Systems Biology Analysis Merging Phenotype, Metabolomic and Genomic Data Identifies Non-SMC Condensin I Complex, Subunit G (NCAPG) and Cellular Maintenance Processes as Major Contributors to Genetic Variability in Bovine Feed Efficiency. *PLoS ONE* **2015**, *10*, e0124574.
[CrossRef] [PubMed]

29. Alexandre, P.A.; Naval-Sanchez, M.; Porto-Neto, L.R.; Ferraz, J.B.S.; Reverter, A.; Fukumasu, H. Systems Biology Reveals NR2F6 and TGFB1 as Key Regulators of Feed Efficiency in Beef Cattle. *Front. Genet.* **2019**, *10*, 230. [CrossRef]

30. Shelton, L.M.; Park, B.K.; Copple, I.M. Role of Nrf2 in protection against acute kidney injury. *Kidney Int.* **2013**, *84*, 1090–1095.
[CrossRef]

31. Chen, Q.M.; Maltagliati, A.J. Nrf2 at the heart of oxidative stress and cardiac protection. *Physiol. Genom.* **2018**, *50*, 77–97.
[CrossRef] [PubMed]

32. Kobayashi, A.; Kang, M.-I.; Watai, Y.; Tong, K.I.; Shibata, T.; Uchida, K.; Yamamoto, M. Oxidative and Electrophilic Stresses Activate Nrf2 through Inhibition of Ubiquitination Activity of Keap1. *Mol. Cell. Biol.* **2006**, *26*, 221–229. [CrossRef] [PubMed]

33. Zhou, N.; Lee, W.R.; Abasht, B. Messenger RNA sequencing and pathway analysis provide novel insights into the biological basis of chickens’ feed efficiency. *BMC Genom.* **2015**, *16*, 195. [CrossRef] [PubMed]

34. Furukawa, M.; Xiong, Y. BTB Protein Keap1 Targets Antioxidant Transcription Factor Nrf2 for Ubiquitination by the Cullin 3-Roc1 Ligase. *Mol. Cell. Biol.* **2005**, *25*, 162–171. [CrossRef]

35. Maher, J.; Yamamoto, M. The rise of antioxidant signaling—The evolution and hormetic actions of Nrf2. *Toxicol. Appl. Pharmacol.* **2010**, *244*, 4–15. [CrossRef]

36. Dinkova-Kostova, A.T.; Abramov, A.Y. The emerging role of Nrf2 in mitochondrial function. *Free. Radic. Biol. Med.* **2015**, *88*, 179–188. [CrossRef]

37. Acosta-Alvear, D.; Zhou, Y.; Blais, A.; Tsikitis, M.; Lents, N.H.; Arias, C.; Lennon, C.J.; Kluger, Y.; Dynlacht, B.D. XBP1 Controls Diverse Cell Type- and Condition-Specific Transcriptional Regulatory Networks. *Mol. Cell* **2007**, *27*, 53–66. [CrossRef]

38. Greene, E.; Flees, J.; Dadgar, S.; Mallman, B.; Orlowski, S.; Dhamad, A.; Rochell, S.; Kidd, M.; Laurendon, C.; Whitfield, H.; et al. Quantum blue reduces the severity of muscle myopathy via modulation of oxygen homeostasis-related genes in broiler chickens. *Front. Physiol.* **2019**, *10*, 1251. [CrossRef]

39. Jacinto, E.; Loewith, R.; Schmidt, A.; Lin, S.; Rüegg, M.A.; Hall, A.; Hall, M.N. Mammalian TOR complex 2 controls the actin cytoskeleton and is rapamycin insensitive. *Nat. Cell Biol.* **2004**, *6*, 1122–1128. [CrossRef]

40. Kim, D.-H.; Sarbassov, D.D.; Ali, S.M.; King, J.E.; Latek, R.R.; Erdjument-Bromage, H.; Tempst, P.; Sabatini, D.M. mTOR Interacts with Raptor to Form a Nutrient-Sensitive Complex that Signals to the Cell Growth Machinery. *Cell* **2002**, *110*, 163–175. [CrossRef]

41. Sarbassov, D.D.; Ali, S.M.; Kim, D.-H.; Guertin, D.A.; Latek, R.R.; Erdjument-Bromage, H.; Tempst, P.; Sabatini, D.M. Rictor, a Novel Binding Partner of mTOR, Defines a Rapamycin-Insensitive and Raptor-Independent Pathway that Regulates the Cytoskeleton. *Curr. Biol.* **2004**, *14*, 1296–1302. [CrossRef]

42. Bentzinger, C.F.; Romanino, K.; Cloëtta, D.; Lin, S.; Mascarenhas, J.B.; Oliveri, F.; Xia, J.; Casanova, E.; Costa, C.F.; Brink, M.; et al. Skeletal muscle-specific ablation of raptor, but not of rictor, causes metabolic changes and results in muscle dystrophy. *Cell Metab.* **2008**, *8*, 411–424. [CrossRef]

43. Nisoli, E.; Clementi, E.; Moncada, S.; Carruba, M.O. Mitochondrial biogenesis as a cellular signaling framework. *Biochem. Pharmacol.* **2004**, *67*, 1–15. [CrossRef]

44. Nisoli, E.; Clementi, E.; Paolucci, C.; Cozzi, V.; Tonello, C.; Sciorati, C.; Bracale, R.; Valerio, A.; Francolini, M.; Moncada, S.; et al. Mitochondrial Biogenesis in Mammals: The Role of Endogenous Nitric Oxide. *Science* **2003**, *299*, 896–899. [CrossRef]

45. Hosotani, M.; Kawasaki, T.; Hasegawa, Y.; Wakasa, Y.; Hoshino, M.; Takahashi, N.; Ueda, H.; Takaya, T.; Iwasaki, T.; Watanabe, T. Physiological and Pathological Mitochondrial Clearance Is Related to Pectoralis Major Muscle Pathogenesis in Broilers with Wooden Breast Syndrome. *Front. Physiol.* **2020**, *11*. [CrossRef]

46. Bottje, W.; Carstens, G.E. Association of mitochondrial function and feed efficiency in poultry and livestock species1. *J. Anim. Sci.* **2009**, *87*, E4–E63. [CrossRef]

47. Papah, M.B.; Brannick, E.M.; Schmidt, C.J.; Abasht, B. Gene expression profiling of the early pathogenesis of wooden breast disease in commercial broiler chickens using RNA-sequencing. *PLoS ONE* **2018**, *13*, e0207346. [CrossRef]

48. Lake, J.A.; Papah, M.B.; Abasht, B. Increased Expression of Lipid Metabolism Genes in Early Stages of Wooden Breast Links Myopathy of Broilers to Metabolic Syndrome in Humans. *Genes* **2019**, *10*, 746. [CrossRef]

49. Guo, J.; Sartor, M.; Karyala, S.; Medvedovic, M.; Kann, S.; Puga, A.; Ryan, P.; Tomlinson, C.R. Expression of genes in the TGF-B signaling pathway is significantly deregulated in smooth muscle cells from aorta of aryl hydrocarbon receptor knockout mice. *Toxicol. Appl. Pharmacol.* **2004**, *194*, 76–89. [CrossRef]

50. Larigot, L.; Juricek, L.; Dairou, J.; Coumoul, X. AhR signaling pathways and regulatory functions. *Biochim. Open* **2018**, *7*, 1–9.
[CrossRef]

51. Gutierrez-Vasquez, C.; Quintana, F. Regulation of the immune response by the aryl hydrocarbon receptor. *Immun. Cell Biol.* **2018**, *48*, 19–29. [CrossRef] [PubMed]

52. Fan, W.; He, N.; Lin, C.S.; Wei, Z.; Hah, N.; Waizenegger, W.; Ming-Xiao, J.; Liddle, C.; Yu, R.T.; Atkins, A.R.; et al. Errg promotes angiogenesis, mitochondrial biogenesis, and oxidative remodeling in PGC1a/b-deficient mice. *Cell Rep.* **2018**, *22*, 2521–2529.
[CrossRef] [PubMed]
<figure><img src="image_98.png" /></figure>
MDPI

St. Alban-Anlage 66

4052 Basel

Switzerland

Tel. +41 61 683 77 34

Fax +41 61 302 89 18

www.mdpi.com

Foods Editorial Office

E-mail: foods@mdpi.com

www.mdpi.com/journal/foods

<figure><img src="image_99.png" /></figure>
<figure><img src="image_100.png" /></figure>
MDPI

St. Alban-Anlage 66

4052 Basel

Switzerland

Tel: +41 61 683 77 34

www.mdpi.com

<figure><img src="image_101.png" /></figure>

ISBN 978-3-0365-5666-6