تحلیل میزان ریسک معرفی ارقام با پتانسیل عملکرد بالا و پایدار در ذرت

. مقدمه

در به­نژادی گیاهان زراعی، ارزیابی آنها در آزمایش‌های چندمحیطی یک گام مهم برای درک پایداری و عملکرد آنها در تصمیم‌گیری‌های مبتنی بر داده است. پاسخ‌ متفاوت ژنوتیپ‌ها در محیط‌های متمایز ناشی از برهمکنش ژنوتیپ در محیط یک چالش شناخته­شده در برنامه­های به­نژادی گیاهان زراعی است. آزمایش‌های چندمحیطی برای درک اثرات برهمکنش ژنوتیپ×محیط بر عملکرد محصول و برای معرفی‌ رقم به یک منطقه هدف معین استفاده می‌شود (Crossa, 1990; Falconer & Mackay, 1996).

معمولاً تنها میانگین عملکرد نیست که شایستگی نسبی ارقام مختلف را تعیین می‌کند، بلکه پایداری عملکرد نیز بسیار مهم می‌باشد. معرفی ارقام در طی دوره‌های مختلف صرفاً بر اساس متوسط عملکرد تا تلفیق روش‌های تجزیه پایداری برای توصیه بهتر رقم صورت ‌گرفته است (Reckling et al., 2021).  

روش‌های تجزیه پایداری مبتنی بر عملکرد مختلف از جمله رویکردهای تک‌متغیره وچندمتغیره برای در شناسایی ارقام برتر در صورت وجود برهمکنش ژنتیکی در دسترس هستند. روش‌های تک‌متغیره مانند مدل رگرسیون فینلی- ویلکینسون، مدل رگرسیون ابرهارت-راسل و واریانس پایداری شوکلا (Shukla, 1972) اطلاعاتی را ارائه می‌دهند که به به­نژادگران اجازه می‌دهد ارقام را بر اساس سطوح سازگاری رتبه‌بندی یا طبقه‌بندی کنند. با این حال، این روش‌ها به‌صراحت نشان نمی‌دهند که چگونه یک شاخص قابل ‌استفاده بر اساس میانگین عملکرد و پایداری ایجاد شود. بعداً مدل‌های چندمتغیره مانند تجزیه به مؤلفه‌های اصلی (PCA)، اثر باقیمانده‌ها در تجزیه واریانس، تجزیه به مؤلفه‌های اصلی، اثر ژنوتیپ به­علاوه اثر ژنوتیپ در محیط به اصطلاح
 GGE Biplot(Yan et al., 2000) ، مدل رگرسیون سایت (Crossa & Cornelius, 1997) و مدل تلفیقی ضربی تحلیلی عاملی (Piepho, 1996) به­طور گسترده مورد استفاده قرار گرفته­اند. تجزیه و تحلیل پایداری عملکرد در سال‌های اخیر اهمیت بیشتری پیدا کرده است، زیرا افزایش تغییرات اقلیمی با کاهش پایداری عملکرد محصول همراه می­باشد
(Muller et al., 2018; Najafi et al., 2018; Ray et al., 2015, Tigchelaar et al., 2018) . برای کشاورزان، پایداری عملکرد که قابل پیش­بینی­کردن درآمد اقتصادی را ممکن ساخته و ریسک را کاهش می­دهد (Reckling et al., 2021) بسیار مهم می­باشد. به­طور معمول، به‌نژادگران گیاهی در انتخاب و توصیه ژنوتیپ‌های آزمایشی به یک محیط هدف همیشه با یک ریسک نهفته‌ای مواجه هستند. اگرچه روش‌های پایداری چندمتغیره، انتخاب هم­زمان پایداری و متوسط عملکرد را در اولویت قرار می­دهند، اما محدودیت‌هایی دارند که باید در نظر گرفته شوند. این روش‌ها اندازه‌گیری عدم قطعیت در پاسخ فنوتیپی به تغییرات محیطی را در نظر نمی‌گیرند که عامل مهمی برای بهبود سازگاری محصول و ارائه توصیه‌ دقیق رقم است.

با توجه به برهمکنش ژنوتیپ در محیط، هر پیش‌بینی عملکرد از تجزیه و تحلیل آماری آزمایش‌های چندمحیطی در معرض یک عدم قطعیت قرار دارد. به­طور خاص، اندازه­گیری یا کمی­کردن میزان ریسک مرتبط با عملکرد در برنامه­های به­نژادی از اهمیت خاصی برخوردار است. ریسک عملکرد به عدم قطعیت مرتبط با عملکرد یک رقم معین در یک محیط اشاره دارد. این عدم قطعیت ناشی از برهمکنش ژنوتیپ×محیط است. برای مدیریت ریسک عملکرد، کمی­کردن احتمال اینکه یک رقم معین از سایر ارقام در مجموعه‌ای از محیط‌های معین بهتر عمل کند، مفید است. این را می­توان با محاسبه احتمال عملکرد بهتر یک رقم از سایر ارقام در یک مجموعه معین از ارقام (Dias et al., 2022) و یا با محاسبه احتمالی (ریسک) که عملکرد یک رقم از یک حد معین کمتر باشد (Eskridge, 1990; Mead et al., 1986; Piepho, 1996) انجام داد. استفاده از روش‌های مبتنی بر احتمال می‌تواند گام مهمی برای افزایش توانایی به‌نژادگران در شناسایی ارقام پایدار برای طیف گسترده‌ای از محیط‌های متمایز یا مناطق تولید باشد. چندین مطالعه در چارچوب متداول (Frequentist framework) برای محاسبه معیاری از ریسک برای رتبه­بندی ژنوتیپ‌ها در آزمایش‌های چندمحیطی استفاده شده است (Annicchiarico, 1992; Eskridge et al., 1991; Mead et al., 1986;
Barah et al., 1981). در این میان، چارچوب بیزی به­ویژه به دلیل توانایی آن در ادغام اطلاعات از آزمایش‌های قبلی، روشی سودمند و جالب (Crossa et al., 2011; Gelman et al., 2013) بوده و با بررسی برهمکنش ژنوتیپ × محیط
 (Crossa et al., 2011) و واریانس‌های پایداری (Shukla, 1972)، درک بهتری از مکانیسم‌های زیربنایی برهمکنش ژنوتیپ × محیط و پیش‌بینی عملکرد ژنوتیپ‌های جدید در محیط‌هایی که هنوز کشت نشده اند، را ممکن می‌سازد.

اخیراًDias et al.  (2022) یک روش نوین بیزی را پیشنهاد کردند که امکان در نظر گرفتن این ریسک‌ها را هنگام انجام انتخاب فراهم می‌کند و از توزیع پسین برای به‌دست­آوردن برآوردهای هامیلتونی مونت­کارلو از احتمال‌های عملکرد و پایداری استفاده می‌کند. ایده‌ اصلی آن‌ها ارزیابی قابلیت پیش‌بینی عملکرد یک ژنوتیپ مورد بررسی از طریق احتمال قرار گرفتن در بین ژنوتیپ‌های منتخب در سازگاری عمومی‌ (حاشیه‌ای، بین محیط‌ها) یا خصوصی (شرطی، درون محیط‌ها یا محیط‌های بزرگ) و احتمال اینکه یک کاندید انتخابی عملکردی ثابت در بین محیط‌ها داشته باشد، است. این روش همچنین احتمالات زوجی را فراهم می‌کند که برای مقایسه مستقیم ژنوتیپ‌های آزمایشی یا ژنوتیپ‌های آزمایشی با ارقام شاهد مفید است.

برای بکارگیری هر روش آماری کارآمد، یک نرم­افزارکامپیوتری مناسب که بتواند نتایج را به صورت بصری ارائه دهد بسیار مهم است. در این راستا، یک بسته جدید و کاربرپسند به نام بسته ProbBreed ارائه شد که ابزاری برای محاسبه ریسک معرفی رقم در آزمایش‌های چندمحیطی (MET) با استفاده از نظریه احتمال بر اساس روش(2022) Dias et al.  است. این بسته به به‌نژادگران اجازه می‌دهد تا احتمال عملکرد یک ژنوتیپ معین از سایر ارقام رقیب را تحت چارچوب بیزی محاسبه کنند. این مدل اصولاً وضعیتی را تقلید می‌کند که معمولاً تولیدکنندگان با آن مواجه هستند، یعنی اولا انتخاب رقم(هایی) که احتمالاً برای فصل زراعی بعدی عملکرد خوبی خواهند داشت و دوماً احتمالات (حاشیه­ای یا مشروط) با توجه به شدت انتخاب محاسبه می‌شود که این نیز بخشی از روال به‌نژادگران گیاهی است. علاوه بر این، در مقایسه با روش‌های مبتنی بر بای‌پلات (Yan et al., 2000)، این روش ساده‌تر است زیرا تنها معیار مورد استفاده برای انتخاب، احتمالات محاسبه شده است.

Malikouski et al. (2024) از مدل احتمالی بیزی برای انتخاب ژنوتیپ‌های برتر و پایدار لیمو ترش تاهیتی استفاده کردند. آن‌ها با استفاده از الگوریتم نمونه‌برداری همیلتونی مونت­کارلو، احتمال عملکرد برتر (ارزش ژنوتیپی برتر) و احتمال پایداری برتر برای هر ژنوتیپ برآورد و اعلام کردند که مدل احتمالی بیزی قابلیت کاربرد در انتخاب ژنوتیپ‌های برتر و پایدار داشته و اطلاعات بیشتری درباره احتمالات مرتبط با عملکرد و پایداری ژنوتیپ فراهم می‌کند. در بررسی آن‌ها، ژنوتیپ‌هایG15، G4، G18  و G11 به عنوان بهترین‌ها از نظرعملکرد شناخته شدند؛ در­حالی­که  G24، G7،G13  و G3 به عنوان پایدارترین‌ها شناسایی شدند. بر اساس نتایج حاصله، به استفاده از مدل‌های احتمالی بیزی در مطالعات آتی تأکید کردند.

این مطالعه، با هدف استفاده از مفاهیم مدل‌های بیزی و روش‌های احتمال تجزیه و تحلیل پایداری برای بررسی برهمکنش ژنوتیپ در محیط (GEI)، به منظور انتخاب هیبرید یا هیبریدهای برتر ذرت، هیبریدهای پایدار، انتخاب هم­زمان برای عملکرد و پایداری و برآورد میزان ریسک انتخاب ارقام انجام شد.

2. روش­شناسی پژوهش

1-2. مواد آزمایشی و طرح آزمایشی

در این بررسی تعداد 37 هیبرید ذرت امیدبخش به همراه یک شاهد تجاری (هیبرید شماره 38) ( جدول 1 ) در قالب طرح بلوک‌های تصادفی در پنج مکان مختلف (کرج (E1)، شیراز (E2)، کرمانشاه (E3)، دزفول (E4)، و گرگان (E5)) مورد ارزیابی قرار گرفتند (جدول 2). هر کرت آزمایشی شامل دو ردیف 5/6 متری با فاصله 75/0 متر از هم، 36/0 متر بین کپه‌ها با دو بوته در هر کپه بود. تراکم کاشت حدود 74500 بوته در هکتار بود. برای اطمینان از سبز شدن کامل بذر، سه بذر در هر کپه کاشته شد. پس از تنک در مرحله چهار تا پنج برگی (حدود 18 روز پس از کاشت)، در هر کپه کاشت تنها دو بوته نگهداری شد. تمام اقدامات توصیه­شده زراعی مانند آبیاری، کنترل علف‌های هرز و کوددهی برای رشد خوب گیاه ذرت در همه مکان‌ها رعایت شد.

مناطق مورد اجرا در ارتفاع بین 143 تا 1600 متر از سطح دریا و در عرض جغرافیایی بین 77/29 تا 50/36 قرار داشتند (جدول 2).

2-2. تحلیل آماری

احتمالات معیارهای عملکرد و پایداری مورد استفاده در این مطالعه توسطDias et al.  (2022) پیشنهاد شده و در بسته نرم­افزاری ProbBreed (Chaves et al., 2024) پیاده­سازی شده است.

جدول 1. اسامی هیبریدهای ذرت مورد مطالعه.

جدول 2. مخشصات جغرافیایی مکان‌های مورد بررسی.

1-2-2. احتمال عملکرد برتر

مجموعه داده­ای را در نظر بگیرید که در آن J ژنوتیپ (j = 1, 2, ..., J) در K محیط ((k= 1, 2, . . ., K با فنوتیپ مشاهده­شده y ارزیابی شده‌اند. فرض کنید Ω زیرمجموعه­ای از ژنوتیپ‌های انتخاب­شده با عملکرد بالا باشد که با توجه به شدت انتخاب تعیین شده باشند (مثلاً هشت ژنوتیپ برتر از مجموع 38 ژنوتیپ با شدت انتخاب 20 درصد). یک ژنوتیپ j در صورتی متعلق به Ω خواهد بود که مقدار حاشیه‌ای ژنوتیپی آن ( ) در مقایسه با همتایانش به اندازه کافی زیاد (یا کم) باشد. با استفاده از نمونه‌های گسسته مونت­کارلو از توزیع‌های خلفی مدل‌ بیزی، وقوع S آزمایش‌ (S = 1, 2,..., S) با مدل‌های بیزی شبیه‌سازی شد. سپس، احتمال تعلق  ژنوتیپ به Ω، نسبت رویدادهای موفقیت‌آمیز آن  به تعداد کل رویدادهای نمونه‌برداری­شده طبق رابطه 1 به­دست آمد:

رابطه (1)                                

که  یک متغیر شاخص است که می‌تواند دو مقدار را به خود بگیرد، اگر  در Sام نمونه باشد آنگاه مقدار یک و در غیر این صورت مقدار صفر به خود می‌گیرد. به­طور مشابه، احتمال شرطی عملکرد برتر را می‌توان به هریک از محیط‌ها تعمیم داد. اگر  زیرمجموعه­ای از ژنوتیپ‌های برتر در محیط kام باشد، به­طوری­که می‌توان احتمال  را به صورت رابطه 2 محاسبه کرد:

 رابطه (2)                     

که در آن  یک متغیر نشانگر است که موفقیت (1) را در صورتی که   در  وجود داشته باشد، و در غیر این صورت شکست (0)  را نشان می‌دهد. توجه داشته باشید که در محاسبه احتمالات شرطی (به عنوان مثال، مشروط به K ام محیط یا مگا محیط)، برهمکنش  ژنوتیپ با  محیط در نظر گرفته می‌شود.

احتمالات جفتی عملکرد برتر را نیز می‌توان در بین یا درون محیط‌ها محاسبه کرد. این آماره احتمال برتری ژنوتیپ J ام نسبت به ژنوتیپ‌های دیگر یا شاهد تجاری مقایسه می‌کند که براساس رابطه­های 3 و 4 محاسبه می‌شوند:

رابطه (3)                         

رابطه (4)              

توجه داشته باشید که معادلات 3 و 4 برای مواقعی که جهت انتخاب مثبت باشد (یعنی هدف افزایش مقدار صفت است) ارایه شده­اند. در صورتی که اگر انتخاب منفی باشد، به سادگی می‌توان جهت > را به < تغییر داد. معادله 3 احتمالات زوجی را بین محیط‌ها و معادله 4 درون محیط‌ها محاسبه می‌کند.

2-2-1-1. احتمال پایداری برتر

احتمال عملکرد برتر، ژنوتیپ‌های با عملکرد بالا را مشخص می‌کند. برای پایداری، احتمال پایداری برتر باید محاسبه گردد. این معیار را می‌توان مستقیماً با روش شوکلا(Shukla, 1972)  مقایسه کرد. یک ژنوتیپ پایدار ژنوتیپی است که واریانس کمی از اثرات برهمکنش ژنوتیپ × محیط [ ] را دارد. با استفاده از همان اصول احتمال که قبلاً توضیح داده شد، احتمال پایداری برتر به صورت رابطه 5 ارائه می‌شود:

رابطه (5)                     

که  نشان می‌دهد که آیا  در Ω وجود دارد (1) یا خیر (0). توجه داشته باشید که این احتمال تنها در میان محیط‌ها قابل محاسبه است زیرا به  بستگی دارد. احتمالات دوتایی پایداری برتر نیز برای پایداری از طریق رابطه 6 محاسبه می‌شود:

رابطه (6)                     

2-2-1-2. احتمال توام عملکرد برتر و پایداری

احتمال توام وقوع رویدادهای مستقل برابر با حاصل ضرب احتمالات فردی آنها است. اثر اصلی ژنوتیپی و واریانس اثر برهمکنش ژنوتیپ × محیط به دلیل طراحی مدل‌های خطی مستقل هستند. بنابراین، احتمال توام عملکرد برتر و پایداری به صورت رابطه 7 محاسبه می‌شود:

رابطه (7)                     

2-2-2. تشخیص همگرایی مدل (Convergence diagnostics)

ضریب کاهش مقیاس (Scale reduction factor ( )) برای ارزیابی سودمندی همگرایی زنجیره مارکوف مونت­کارلو(MCMC) استفاده شد. این آماره نشان می‌دهد که آیا زنجیره‌ها به اندازه کافی ادغام شده­اند و پارامترهای مدل تخمین­زده­شده، به توزیع پایدار رسیده یا نه. هر چه به عدد یک نزدیکتر باشد، کیفیت اختلاط و همگرایی بیشتر است (Gelman et al., 2013). مقادیر بیشتر نشان می‌دهد که به تکرارهای بیشتر نیاز است. همچنین تجزیه و تحلیل گرافیکی انجام شد تا به صورت بصری ارایه شود که چگونه داده‌های تولیدشده توسط مدل با فرآیند تولید واقعی داده‌های مشاهده­شده هماهنگ است. این فرایند با ایجاد نمونه‌های ( ) از مدل‌ برازش­شده با استفاده از نمونه‌گیری اجدادی از توزیع توام شرطی و سپس ترسیم این نمونه‌ها در برابر داده‌های مشاهده‌شده می‌باشد. علاوه بر این، از مقادیر p پیش‌بینی پسین برای سنجش میزان شباهت معیارهای آماری (حداکثر، حداقل، میانه، میانگین و انحراف استاندارد) داده‌های تولیدشده توسط مدل‌ برازش­شده با داده‌های مشاهده­شده استفاده شد
 (Gelman et al., 2013). به عنوان مثال، هنگام در نظر گرفتن حداکثر آمار، مقدار p بیزی را به صورت زیر تعریف شد:
 که در آن T آزمون آماری است. درجه شباهت بین آماره‌های به­دست­آمده از داده‌های تولید شده و آماره حاصل از داده‌های مشاهده­شده با نزدیک­شدن مقادیر p بیزی به 5/0 افزایش می‌یابد.

3. یافته­های پژوهش و بحث

3-1. ارزیابی مدل برای رویکرد بیزی: اثرات پسین و نکویی برازش (Posterior effects and goodness-of-fit diagnostics) در آزمایش‌های چندمحیطی

مدل‌های بیزی با استفاده از الگوریتم‌هایی مانند نمونه‌گیری بدون برگشت (No-U- trun Sampler) و مونت­کارلو هامیلتونی (Hamiltonian Monte Carlo )، کارایی محاسباتی را بهبود می‌بخشند و نیاز به تعیین تعداد به‌روزرسانی‌های جهشی را از بین می‌برند. این امر از رفتار گام‌های تصادفی (Random-walk behavior) جلوگیری کرده و در نتیجه کارایی محاسباتی را افزایش می‌دهد (Hoffman & Gelman, 2014, Nishio & Arakawa, 2019).

توزیع احتمال پسین پارامترهای مدل در این مطالعه با استفاده از زبان برنامه‌نویسی احتمالی استخراج شد. برای این منظور، از رابط (Hoffman & Gelman, 2014; Carpenter et al., 2017) RStan بهره گرفته شد. در این فرآیند، چهار زنجیره، هر کدام شامل ۲۰۰۰ نمونه، تنظیم شد. به منظور بهبود کیفیت نتایج، ۵۰ درصد از نمونه‌های اولیه به عنوان نمونه‌های سوخته (burn-in) حذف شدند که این امر منجر به مجموع ۴۰۰۰ نمونه برای تحلیل‌های بعدی شده و برای محاسبه احتمال‌های پسین پیشنهاد شد (جدول 3).

ضریب کاهش مقیاس بالقوه (R̂) برای ارزیابی کیفیت نمونه‌های زنجیره مارکوف مونت­کارلو (MCMC) استفاده شد. این آماره واریانس درون و بین زنجیره‌ها اندازه‌گیری می‌کند و نشان‌دهنده درجه تعادل بین زنجیره‌ها برای پارامترهای مدل می‌باشد. در واقع، این آماره، نسبت میانگین واریانس نمونه‌ها در زنجیره‌ها به واریانس بین زنجیره‌ها است. بنابراین مقادیر نزدیک‌تر به یک نشان‌دهنده مشابه بودن این واریانس‌ها است که نشان‌دهنده هم‌گرایی خوب بین نمونه‌های زنجیره‌های مختلف در فضای پارامتر است،که مطلوب است (Fabreti & Hohna, 2022). برعکس، مقدار R̂ بیشتر از یک نشان­دهنده اختلاط ضعیف زنجیره‌ها است
(Gelman et al., 2013). مدل ما دارای آماره R̂ برابر با یک بود که نشان­دهنده همگرایی موثر مدل‌های ارزیابی شده است (جدول 3).

برای ارزیابی میزان شباهت داده‌های تولید شده توسط مدل به فرآیند تولید واقعی داده‌های مشاهده‌شده، نمونه‌هایی از مدل برازش شده با استفاده از نمونه‌گیری اجدادی از توزیع توام مشروط تولید و نتایج آن‌ها در مقابل داده‌های مشاهده‌شده ترسیم گردید. آماره‌های مورد استفاده شامل حداکثر، حداقل، میانه، میانگین و انحراف معیار داده‌های تولید شده بود که در مقایسه با داده‌های مشاهده‌شده از طریق مقادیر p پیش‌بینی پسین ارزیابی شدند (Gelman et al., 2013). بیزین “p-values” به عنوان احتمال یک آماره آزمون (مثلاً میانگین) در داده‌های تولیدشده نسبت به داده‌های واقعی محاسبه می‌شود. اگر داده‌های تولیدشده شباهت زیادی به داده‌های مشاهده‌شده داشته باشند، انتظار می‌رود که مقادیر p بیزی از مقادیر حدی (99/0 یا 01/0) دور باشد
 (Gelman et al., 2013). بنابراین مقادیر p بیزی نزدیک به 5/0 مطلوب بوده و نشان دهنده شباهت بیشتر بین داده‌های تولیدشده و مشاهده­شده است (Dias et al., 2022). در این بررسی، آماره آزمون با مقایسه داده‌های مشاهده­شده و تولیدشده، مقادیر p بیزی را برای آماره‌های مورد نظر (حداکثر، حداقل، میانه، میانگین و انحراف معیار) در محدوده 2694/0 تا 7722/0 قرار داشت (جدول 3). این یافته‌ها نشان می‌دهد که داده‌های تولیدشده توسط مدل به­طور قابل توجهی با داده‌های مشاهده­شده مشابهت داشت.

اطلاعات حاصل از توزیع پسین، بیشینه‌های پسین و داده‌های تولیدشده توسط مدل برای ساخت نمودارهایی استفاده می‌شود که نمای کلی از نکویی برازش مدل را ارائه می‌دهند. بالاترین بازه چگالی پسین (HPDI) با استفاده از احتمالات محاسبه شد تا دامنه پارامتر مورد نظر که درصد معینی از احتمال پسین را در برمی‌گیرد، تعیین شود (شکل ۱).

جدول 3. پارامتر نکویی برازش " بیزین p-values " برای آماره‌های آزمون (حداکثر، حداقل، میانه، میانگین و انحراف استاندارد)، تعداد مؤثر پارامترها، ضریب کاهش مقیاس بالقوه (R̂) و اندازه نمونه مؤثر.

نمودار هیستوگرام (شکل ۱-الف) و چگالی (شکل ۱-ب) که بصری‌سازی اثر پسین توزیع را فراهم می‌کنند، ابزارهای مفیدی برای تشخیص مشکلات همگرایی هر زنجیره هستند. زمانی که توزیع پسین یک مؤلفه واریانس متقارن باشد، نتایج و استنباط‌های حاصل از آن به انتخاب معیار گرایش مرکزی (مانند میانگین یا میانه) حساس نخواهد بود (شکل 1-الف). شکل ۱-ب به ارزیابی این موضوع کمک می‌کند که آیا مدل قادر به تولید داده‌هایی با توزیع مشابه داده‌های واقعی بوده است یا خیر. نکته حائز اهمیت این است که چگالی داده‌های تولیدشده توسط مدل بیزی باید روند چگالی داده‌های واقعی را دنبال کند. این امر نشان‌دهنده اثربخشی مدل در بازتولید توزیع داده‌های مشاهده‌شده از طریق داده‌های تولیدشده است و قابلیت اعتماد آن در شناسایی الگوهای زیربنایی را برجسته می‌سازد (شکل ۱-ب).

شکل 1. هیستوگرام اثرات اصلی ژنوتیپی پسین (الف) و نمودار چگالی داده‌های تولیدشده در مقایسه با توزیع داده‌های واقعی (ب).

3-2. برآورد مؤلفه‌های واریانس بیزی برای داده‌های آزمایش‌های چندمحیطی

برآورد مؤلفه‌های واریانس برای اثرات ژنوتیپ (G)، مکان (L) و برهمکنش آن‌ها (G×L)، به همراه انحراف معیار (SD)، خطای استاندارد ساده (Naive SE) و بازه‌های چگالی پسین بالا (حد پایین HPD 05/0 و حد بالا HPD 95/0)، با استفاده از مدل بیزی انجام شد (جدول ۴). مؤلفه‌های واریانس برآوردشده برای اثرات مذکور به­ترتیب 55/0 برای ژنوتیپ، 528/71 برای مکان و 736/0 برای برهمکنش G×L بود (جدول ۴). این برآوردها اطلاعات خوبی در مورد تغییرات مربوط به هر اثر در مدل بیزی ارائه می‌دهند. در این میان "مکان" بالاترین میزان واریانس از واریانس کل را به خود اختصاص داد. این امر نشان می‌دهد که عامل "مکان" نقش بسیار زیادی در ایجاد تغییرات کل در مدل ایفا می‌کند.

3-3. احتمال عملکرد و پایداری برتر ژنوتیپی هیبریدهای ذرت

احتمال عملکرد برتر و پایداری نامزدهای انتخاب با شدت انتخاب ۲۰ درصد به منظور افزایش میانگین عملکرد دانه در پانل ژنوتیپ‌های انتخاب شده محاسبه گردید. شکل ۲-الف نمای کلی از عملکرد نامزدهای انتخاب را در محیط‌ها ارائه می‌دهد که در یک نمودار کاترپیلار (Caterpillar plot) شامل اثرات اصلی ژنوتیپی پسین و بازه‌های HPD مربوطه است. توزیع پسین ارزش‌های ژنتیکی هیبریدهای ذرت(gj)  الگوی همپوشانی متغیری را بین بازه‌های چگالی پسین بالا(HDI)  آن‌ها نشان می­دهد. هیبریدهای شماره 15 و 13 بالاترین ارزش‌های ژنتیکی را داشتند و هیبریدهای شماره 37 و 1 پایین‌ترین عملکرد ژنتیکی را داشتند (شکل ۲-الف). بیشینه مقادیر پسین معادل BLUPهای تجربی حاشیه‌ای مدل‌های خطی مختلط در روش متداول با فرض اثرات ژنوتیپی مستقل هستند. از­آنجایی­که احتمال عملکرد بر اساس ارزش‌های ژنتیکی برآوردشده محاسبه شد، انتظار می‌رفت که رتبه‌بندی هیبریدهای عملکرد برتر با نتایج حاصل از مقادیر بازه چگالی پسین بالا (HDI) همسو باشد (شکل 2-الف و ب). هیبریدهای شماره 15، 13، 36، 12 و 24 همچنین بیشترین احتمال عملکرد برتر را نشان دادند (شکل ۱-ب). هیبریدهای شماره 29، 4، 25، 1 و 37 در هیچ یک از نمونه‌ها انتخاب نشدند، در­حالی­که ۳۳ هیبرید باقی‌مانده در حداقل برخی نمونه‌ها انتخاب شدند (شکل ۱-ب). شایان ذکر است که احتمال داشتن عملکرد برتر برای هیبرید نامزد شماره 15 برابر با 5/95 درصد بود (شکل ۳-ب). این امر بیانگر ریسک 5/4 درصدی داشتن عملکرد ضعیف برای این هیبرید است، مشروط بر شدت انتخاب [ ]. تجزیه و تحلیل احتمال پایداری برتر نشان داد که در مکان‌ها، هیبریدهای شماره 11، 12 و 1 بالاترین احتمال پایداری برتر را داشتند، به این معنی که عملکرد آنها تغییرات کمتری در محیط‌های مختلف داشتند (شکل 2-پ).

جدول ۴. برآورد مقادیر اجزای واریانس داده‌های چندمحیطی، انحراف معیار (sd)، خطای استاندارد ساده (Naive SE) و بازه چگالی پسین بالا (HPD) (حد پایین HPD 0/05 و حد بالا HPD 0/95)

به­طور کلی، در این مطالعه احتمال شناسایی یک هیبرید ذرت با کمترین تغییرات در پایداری نسبتاً پایین بود. در میان هیبریدهای مورد مطالعه تنها هیبرید شماره 11 (پایدارترین هیبرید در بین هیبریدهای مورد مطالعه) با احتمال حدود ۴۵ درصد شانس وقوع پایداری را در محیط‌های مورد مطالعه داشت (شکل 2-پ).

این تجزیه و تحلیل تأثیر معیارهای احتمالی مختلف بر نتایج را برجسته می‌کند و بر ضرورت تعیین معیارهای هدف روشن قبل از انجام تجزیه و تحلیل‌ها تأکید می‌کند. اگر عملکرد مهم‌ترین اولویت باشد، شکل ۱-ب باید مورد استناد قرار گیرد و برعکس، اگر پایداری هدف نهایی باشد، باید اولویت به شکل 2-پ داده شود.

در این تحلیل، علاوه بر احتمال‌های فردی، احتمال‌های جفتی برای مقایسه بین ژنوتیپ‌ها محاسبه شده است (شکل 3-الف و ۳-ب). این تحلیل تفسیر احتمالی از اینکه یک هیبرید چقدر احتمال دارد که در بین هیبریدهای برتر قرار گیرد، بر اساس شدت انتخاب از پیش تعیین­شده، ارائه می‌دهد. احتمال‌های جفتی عملکرد برتر در یک نقشه حرارتی نمایش داده شده‌اند تا امکان مقایسه‌ای عملکرد ژنتیکی هیبریدها فراهم شود (شکل ۳-الف). مطابق با احتمال حاشیه‌ای عملکرد برتر، هیبرید شماره 15 نیز به عنوان بهترین رقم در مقابل سایرین معرفی شده است. هیبرید شماره 37 بدترین عملکرد مقایسه‌ای جفتی را داشت. فرض کنید که هیبرید شماره 15 یک هیبرید آزمایشی امیدبخش باشد و ما می‌خواهیم بررسی کنیم که آیا عملکرد بهتری نسبت به شاهد تجاری ( هیبرید شماره 38) دارد یا خیر. در تمام مکان‌ها، هیبرید شماره 15 در ۹۸ درصد مواقع بهتر از هیبرید شماره 38 عمل کرد (شکل 3-الف) و در ۶۰ درصد مواقع پایداری عملکرد کمتری نسبت به هیبرید شماره 38 داشت (شکل 3-ب). بنابراین، هیچ مدرکی برای فرض­کردن اینکه هیبرید شماره 15 در هر دو جنبه پایداری و عملکرد بهتر از شاهد تجاری ( هیبرید شماره 38) بود، وجود ندارد. در نهایت، اگر به‌نژادگر بخواهد ژنوتیپ‌هایی را شناسایی کند که به­طور همزمان عملکرد و پایداری بالایی داشته باشند، شکل 3-پ باید مورد تحلیل قرار گیرد، زیرا شامل احتمال توام عملکرد و پایداری برتر است. مهم است که توجه شود که به سختی می‌توان یک ژنوتیپ را یافت که در تمام معیارهای احتمالی بهترین ژنوتیپ باشد. بر اساس انتخاب همزمان برای عملکرد و پایداری برتر با استفاده از احتمال توام عملکرد و پایداری برتر، هیبریدهای شماره 12، 13 و 15 به عنوان بهترین هیبریدها شناسایی شدند.

شکل ۲. بالاترین بازه چگالی پسین (HPD) اثرات اصلی ژنوتیپی پسین (الف)، احتمال عملکرد برتر در محیط‌های مختلف (ب) و احتمال پایداری برتر در مکان‌های مختلف (پ). نقاط موجود در (الف) نمایانگر بیشینه‌های پسین هستند و خطوط ضخیم و نازک در (الف) به ترتیب بازه‌هایHPD، 95 درصد و 5/97 درصد را نمایش می‌دهند. محور x در شکل­های (ب) و (پ) بر اساس احتمال‌های محاسبه­شده به­ترتیب نزولی مرتب شده است.

در فرآیند به­نژادی محصولات زراعی، اطلاعات ضروری برای انتخاب، درک این نکته است که یک هیبرید بالقوه در کجا ممکن است برتری یا عدم توفیق در عملکرد داشته باشد. احتمال عملکرد برتر در محیط‌ها می‌تواند در این فرآیند تصمیم‌گیری به به‌نژادگر کمک کند (شکل ۴). هیبرید شماره 15، به عنوان بهترین هیبرید در محیط‌های مختلف معین شده، ممکن است در مکان E3  شکست بخورد، اما احتمال موفقیت آن در سایر مکان‌ها بیشتر است (شکل ۴). علاوه بر این، اگر مکان E3 برای تولید محصول مهم باشد، هیبریدهای شماره 13 و 36 نیز می‌توانند به این مکان توصیه شوند.

شکل ۳. (الف) احتمال‌های جفتی عملکرد برتر در مکان‌های مختلف، (ب) پایداری برتر در مکان‌های مختلف و (پ) احتمال توام عملکرد و پایداری برتر. نقشه‌های حرارتی در شکل­های (الف) و (ب) احتمال اینکه ژنوتیپ‌های موجود در محور x نسبت به ژنوتیپ‌های موجود در محور y برتر باشند را نشان می­دهد.

شکل 4. نقشه‌های حرارتی احتمال خصوصی عملکرد برتر در داخل مکان‌ها.

در برنامه‌های به‌نژادی، در نظر گرفتن اثر متقابل ژنوتیپ در محیط (GEI) در ارزیابی ژنتیکی محصولات زراعی بسیار حیاتی است؛ به این دلیل که بیان ژن به عوامل محیطی مختلف در محیط‌های گوناگون به­طور متفاوت واکنش نشان می‌دهد
 (Ferreira et al., 2020). رویکرد احتمالی بیزی روش‌های متداول مورد استفاده برای مطالعه برهمکنش ژنوتیپ × محیط را ساده‌سازی کرده و با استفاده از روش‌های احتمالاتی، دیدگاه‌های جدیدی را در مورد انعطاف­پذیری فنوتیپی ارائه می‌دهد. روش‌های قبلی توسط (1986) Mead et al. بعداً توسطPiepho  (1996) اصلاح شد. روش­های Eskridge (1990)،
Eskridge & Mumm (1992) و (2000) Piepho، رویکردهای مبتنی بر احتمال هستند که گزینه‌های جایگزین دیگری را برای کمک‌های لازم به به­نژادگران در انتخاب و توصیه ارقام ارائه می‌دهند. مدل‌های پیشنهادی به­طور کامل از توزیع پسین به­دست­آمده با الگوریتم نمونه­گیری بدون برگشت (Gelman et al., 2013) برای محاسبه احتمال‌های پسین استفاده می‌کنند. این ویژگی به­طور مستقیم عدم قطعیت هر پارامتر را لحاظ می‌کند. به طور نظری، به دلیل ویژگی‌های زنجیره مارکوف زنجیره مونت­کارلو (MCMC)، همانطور که اندازه نمونه زنجیره MCMC به طور نامحدود افزایش می‌یابد، احتمال‌ها با استفاده از توزیع احتمال تجمعی پسین ناشناخته از مدل‌های بیزی غیر همبسته به مقدار دقیق نزدیک می‌شوند (Goodfellow et al., 2016;
Gelman et al., 2013).(2021) Dias et al.  مفاهیم نظری و عملی استفاده از روش‌های احتمالاتی برای کاهش ریسک تصمیم‌گیری با استفاده از مدل‌های بیزی غیر همبسته برای تقسیم برهمکنش ژنوتیپ × محیط در یک برنامه به‌نژادی را بررسی کردند. آن‌ها همخوانی خوبی بین رویکرد احتمالاتی بیزی و مدل‌های کلاسیک برهمکنش ژنوتیپ × محیط در نتایج به­دست­آمده، یافتند. برای بهره‌گیری از مزایای ذکر شده، رویکرد احتمالاتی بیزی برای درک برهمکنش ژنوتیپ × محیط به منظور ارائه توصیه‌های ارقام در محیط‌های مختلف، با یک مجموعه داده فنوتیپی ذرت با ۳۸ هیبرید ذرت که در پنج محیط در ایران ارزیابی شده بودند، به کار گرفته شد. هیبریدهایی با بالاترین احتمال عملکرد برتر، به­ویژه هیبریدهای شماره 15 و 13، به عنوان نامزدهای قوی برای معرفی به کشاورزان در نظر گرفته شدند. این ژنوتیپ‌ها که از طریق روش‌های احتمالاتی شناسایی شده‌اند، دارای آللی هستند که سازگاری را به شرایط محیطی متغیر در محیط‌های مختلف القا کرده و به آن‌ها امکان حفظ عملکرد خوب و پایدار را می‌دهند. تحلیل احتمال‌های پایداری برتر هیبریدها نشان داد که هیبریدهای شماره 11، 12 و 1 بالاترین احتمال‌های پایداری برتر را در مکان‌های مختلف داشتند که نشان­دهنده عملکرد پایدارتر آنها در محیط‌های مختلف است.

یکی از گلوگاه‌های آزمایش‌های چندمحیطی (MET) تعیین معیار یا معیارهای قابل اعتماد برای مقایسه ریسک‌های مرتبط با انتخاب یا معرفی ارقام مختلف است. به­طور معمول، عملکرد یک رقم شاهد به عنوان مبنایی برای چنین مقایسه‌هایی در نظر گرفته می‌شود (Eskridge & Mumm, 1992). با این حال، عدم پایداری عملکرد رقم شاهد در محیط‌های مختلف می‌تواند به نتیجه‌گیری گمراه‌کننده درباره مناسب­بودن رقم و ارزیابی ریسک منجر شود (Lin & Binns, 1988). برای رفع این حساسیت،
Lin & Binns  (1988) شاخصی را بر اساس اختلاف با حداکثر عملکرد ژنوتیپ‌ها در هر محیط پیشنهاد کردند. این رویکرد روشی استاندارد‌تر برای ارزیابی عملکرد ارقام با تمرکز بر عملکرد نسبی به جای اندازه‌گیری‌های مطلق ارائه می‌دهد.

مقایسه‌های مدل‌ احتمالاتی توام در میان ژنوتیپ‌ها، با ارائه چارچوبی جامع برای ارزیابی سازگاری ارقام اعم از سازگاری عمومی و خصوصی، امکان درک دقیق‌تری از عملکرد ارقام با در نظر گرفتن پتانسیل عملکرد و پایداری در محیط‌ها فراهم می‌کنند. با استفاده از احتمال‌های توام، محققان می‌توانند بهتر تفسیر کنند که یک هیبرید چقدر احتمال دارد در بین گزینه‌های برتر بر اساس شدت انتخاب از پیش تعیین­شده قرار گیرد. شاخص ایمنی-اول بیزی (Eskridge, 1990) (Safety-first index) ارزیابی احتمالاتی از عملکرد ژنوتیپ را با در نظر گرفتن هم پتانسیل عملکرد و هم پایداری در محیط‌های مختلف ارائه می‌دهد (شکل ۳-الف). ژنوتیپ‌هایی با مقادیر بالا در این شاخص ترجیح داده می‌شوند؛ زیرا نشان‌دهنده ترکیبی از عملکرد بالا و عملکرد پایدار هستند که آن‌ها را برای سازگاری عمومی یا اختصاصی مناسب می‌کند. با انتخاب همزمان برای عملکرد و پایداری برتر بر اساس احتمال توام عملکرد و پایداری برتر، هیبریدهای شماره 12، 13 و 15به عنوان هیبریدهای برتر شناخته شدند. روش مطلوب قادر به انتخاب ژنوتیپ‌های پرمحصول و پایدار در محیط‌ها همواره یک چالش بوده است. در راستای این هدف، Kang (1988) روش مجموع رتبه‌ها را پیشنهاد کرد،Eskridge  (1990) قاعده ایمنی-اول را معرفی کرد وAnnicchiarico  (1992) شاخص ریسکی را ارائه داد که به­طور همزمان میانگین ارقام و پایداری را در نظر می‌گیرد.

انجام مقایسه‌های جفتی همه‌جانبه بین ژنوتیپ‌ها برای به­دست­آوردن درک جامعی از عملکرد کلی مقایسه‌ای آن‌ها، به­ویژه در آزمایش‌های پیشرفته‌تر قبل از معرفی رقم، ضروری است. این رویکرد به محققان امکان می‌دهد نقاط قوت و ضعف نسبی هر ژنوتیپ را ارزیابی کنند و فرآیند انتخاب و تصمیم‌گیری آگاهانه‌تری داشته باشند. احتمال‌هایی که با استفاده از نمونه‌های پسین محاسبه می‌شوند، عدم قطعیت پارامتر را در نظر می‌گیرند و ارزیابی قابل اعتمادتری از عملکرد ژنوتیپ ارائه می‌دهند. این چارچوب احتمالاتی فرآیند به‌نژادی را با تلفیق بینش‌های ژنتیکی با سازگاری محیطی بهبود می‌بخشد که در نهایت به توسعه ارقام محصول پایدار و پرمحصول منجر می‌شود. همانطور که قبلاً ذکر شد، بر اساس احتمال حاشیه‌ای عملکرد برتر، هیبرید شماره 15 به عنوان رقم برتر در بین هیبریدهای مورد مطالعه شناسایی شد. این هیبرید آزمایشی امیدبخش بوده و بررسی بیشتر برای تعیین عملکرد آن نسبت به شاهد تجاری ( هیبرید شماره 38) ضروری است. در بین مکان‌های مورد مطالعه، هیبرید شماره 15 در ۹۸ درصد مواقع بهتر از هیبرید شماره 38 عمل کرد (شکل ۳-الف). با این حال، مهم است که توجه داشته باشیم هیبرید شماره 15 در ۶۰ درصد موارد پایداری کمتری نسبت به هیبرید شماره 38 داشت (شکل ۳-ب). در نتیجه، شواهدی برای فرض این که هیبرید شماره 15 به­طور کامل از نظر هم پایداری و هم عملکرد بهتر از شاهد تجاری (هیبرید شماره 38) است، وجود ندارد. برای به‌نژادگرانی که به دنبال شناسایی ژنوتیپ‌هایی هستند که در هر دو جنبه عملکرد و پایداری برتر هستند، تحلیل شکل ۳-پ حیاتی است. این شکل احتمال توام عملکرد و پایداری برتر را نشان داده و نمای جامعی از عملکرد ژنوتیپ ارائه می‌کند. درک این نکته ضروری است که احتمال ندارد یک ژنوتیپ در تمام معیارهای احتمالی بهترین باشد که پیچیدگی ارزیابی ژنوتیپ در شرایط محیطی متغیر را برجسته می‌سازد. مدل بیزی اطلاعاتی در مورد عدم قطعیت یک ژنوتیپ متعلق به مجموعه ژنوتیپ‌های با سازگاری خصوصی یا عمومی ارائه می‌دهد. این اطلاعات در فرآیند انتخاب مهم است، زیرا به به‌نژادگران این امکان را می‌دهد که تصمیمات آگاهانه‌ای بر اساس احتمال قرار گرفتن یک ژنوتیپ در بین بهترین‌ها اتخاذ کنند.

سوال اصلی مورد بررسی در این تحقیق، کاربرد روش‌های احتمالاتی مبتنی بر مدل‌های بیزی برای بهبود دقت توصیه‌های رقم برای سازگاری‌های عمومی و خصوصی است. این جنبه بسیار مهمی در فرآیند تصمیم‌گیری هر برنامه به‌نژادی است. احتمال سازگاری خصوصی بینش‌های ارزشمندی در مورد روندهای سازگاری ژنوتیپ با شناسایی هیبریدهایی با انعطاف­پذیری مطلوب متناسب با مناطق به‌نژادی خاص ارائه می‌دهد. این رویکرد احتمالاتی به به‌نژادگران امکان می‌دهد تصمیمات آگاهانه‌تری درباره مناسب­بودن ارقام برای محیط‌های هدف اتخاذ کنند. علاوه بر این، تحلیل‌های سازگاری عمومی جفتی درک عملکرد مقایسه‌ای بین ژنوتیپ‌های هدف در محیط‌های مختلف را بهبود می‌بخشند. این تحلیل‌ها اغلب از طریق نقشه‌های حرارتی کاربرپسند ارائه می‌شوند که تفسیر و تصمیم‌گیری را تسهیل می‌کنند. در نهایت، این یافته‌ها درک روشن‌تر و جامع‌تری از برهمکنش ژنوتیپ درمحیط با شدت انتخاب تعریف­شده ارائه می‌دهند. با تلفیق این اطلاعات در فرآیند تصمیم‌گیری، به‌نژادگران می‌توانند توصیه‌های رقم آگاهانه‌تری را در آزمایش‌های چندمحیطی ارائه دهند.

4. نتیجه ­گیری

با استفاده از مدل‌های بیزی مبتنی بر احتمال می‌توان احتمال عملکرد برتر یک ژنوتیپ و احتمال‌های جفتی عملکرد برتر بین ژنوتیپ‌ها را در و درون محیط‌های مختلف تخمین زد. هیبریدهای 15، 13، 36، 12 و 24 به عنوان هیبریدهای برتر از نظر عملکرد و هیبریدهای شماره 11، 12 و 14به عنوان هیبریدهای برتر از نظر پایداری‌ شناخته شدند. بنابراین، می‌توان نتیجه­گیری کرد مدل‌های بیزی احتمالی می‌توانند به توصیه‌های دقیق‌تری در مورد ژنوتیپ‌های ذرت که در چند محیط ارزیابی شده‌اند، کمک کنند، زیرا این مدل‌ها امکان تفسیر مستقیم و دقیقی از عملکرد و پایداری ژنوتیپ‌های نامزد فراهم می‌آورند.

4. منابع

Annicchiarico, P. (1992). Cultivar adaptation and recommendation from alfalfa trials in Northern Italy. Journal of Genetics and Breeding (Italy), 46, 269-278.

Barah, B.C., Binswanger, H.P., Rana, B.S., & Rao, N.G.P. (1981). The use of risk aversion in plant breeding: Concept and application. Euphytica, 30, 451–458.

Carpenter, B., Gelman, A., Hofman, M.D., Lee, D., Goodrich, B., Betancourt, M., Brubaker, M., Guo, J., Li, P., & Riddell, A. (2017). Stan: Probabilistic programming language. Journal of Statistical Software, 76(1), 1–32.

Chaves, S.F.S., Evangelista, J.S.P.C., Trindade, R.S., Dias, L.A.S., Guimarães, P.E., Guimarães, L.J.M., Alves, R.S., Bhering, L.L., & Dias, K.O.G. (2023). Employing factor analytic tools for selecting high-performance and stable tropical maize hybrids. Crop Science, 63, 1114–1125.

Crossa, J. (1990). Statistical analyses of multilocation trials. Advances in Agronomy, 44, 55–85.

Crossa, J., & Cornelius, P.L. (1997). Sites regression and shifted multiplicative models clustering of cultivar trials sites under heterogeneity of error variance. Crop Science, 37, 406–415.

Crossa, J., Perez-Elizalde, S., Jarquin, D., Cotes, J.M., Viele, K., Liu, G., & Cornelius, P.L. (2011). Bayesian estimation of the additive main effects and multiplicative interaction model. Crop Science, 51, 1458–1469.

Dias, K.O.G., Santos, J.P.R., Krause, M.D., Piepho, H.P., Guimarães, L.J.M., Pastina, M.M., & Garcia, A.F. (2022). Leveraging probability concepts for cultivar recommendation in multi-environment trials. Theoretical and Applied Genetics, 135(6), 1385–1399.

Eskridge, K., Byrne, P., & Crossa, J. (1991). Selection of stable varieties by minimizing the probability of disaster. Field Crops Research, 27(2), 169–181.

Eskridge, K.M. (1990). Selection of stable cultivars using a safety-first rule. Crop Science, 30(2), 369–374.

Eskridge, K.M., & Mumm, R.F. (1992). Choosing plant cultivars based on the probability of outperforming a check. Theoretical and Applied Genetics, 84(4), 494–500.

Fabreti, L.G., & Höhna, S. (2022). Convergence assessment for Bayesian phylogenetic analysis using MCMC simulation. Methods in Ecology and Evolution, 13(1), 77–90.

Falconer, D.S., & Mackay, T.F.C. (1996). Introduction to quantitative genetics (4th ed.). Longmans Green.

Gelman, A., Carlin, J.B., Stern, H.S., Dunson, D.B., Vehtari, A., & Rubin, D.B. (2013). Bayesian Data Analysis (3rd ed.). Chapman and Hall/CRC.

Goodfellow, I., Bengio, Y., & Courville, A. (2016). Machine learning basics (Chapter 5) in Deep Learning (MIT Press).

Hoffman, M.D., & Gelman, A. (2014). The No-U-Turn sampler: Adaptively setting path lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research, 15(1), 1593–1623.

Kang, M.S. (1988). A rank-sum method for selecting high-yielding stable corn genotypes. Cereal Research Communications, 16(2), 113–115.

Lin, C.S., & Binns, M.R. (1988). A superiority measure of cultivar performance for cultivar × location data. Canadian Journal of Plant Science, 68(2), 193–198.

Malikouski, R.G., Ferreira, F.M., Chaves, S.F.D.S., Couto, E.G.D.O., Dias, K.O.D.G., & Bhering, L.L. (2024). Recommendation of Tahiti acid lime cultivars through Bayesian probability models. Plos one.

Mead, R., Riley, J., Dear, K., & Singh, S.P. (1986). Stability comparison of intercropping and monocropping systems. Biometrics, 42(2), 253–266.

Müller, C., Elliott, J., Pugh, T.A.M., Ruane, A.C., Ciais, P., Balkovic, J., Deryng, D., Folberth, C., Cesar Izaurralde, R., Jones, C.D., Khabarov, N., Lawrence, P., Liu, W., Reddy, A.D., Schmid, E., & Wang, X. (2018). Global patterns of crop yield stability under additional nutrient and water inputs. Plos One.

Najafi, E., Devineni, N., Khanbilvardi, R.M., & Kogan, F. (2018). Understanding the changes in global crop yields through changes in climate and technology. Earth’s Future, 6(3), 410–427.

Nishio, M., & Arakawa, A. (2019). Performance of Hamiltonian Monte Carlo and No-U-Turn sampler for estimating genetic parameters and breeding values. Genetics Selection Evolution, 51(1), 73.

Piepho, H.P. (1996). A simplified procedure for comparing the stability of cropping systems. Biometrics, 52(2), 315–320.

Ray, D.K., Gerber, J.S., MacDonald, G.K., & West, P.C. (2015). Climate variation explains a third of global crop yield variability. Nature Communications, 6, 5989.

Tigchelaar, M., Battisti, D.S., Naylor, R.L., & Ray, D.K. (2018). Future warming increases probability of globally synchronized maize production shocks. Proceedings of the National Academy of Sciences, 115(26), 6644-6649.

Reckling, M., Ahrends, H., Chen, T.W., Eugster, W., Hadasch, S., Knapp, S., & Döring, T.F. (2021). Methods of yield stability analysis in long-term field experiments: A review. Agronomy for Sustainable Development, 41(1), 1-28.

Shukla, G.K. (1972). Some statistical aspects of partitioning genotype environmental components of variability. Heredity, 29, 237-245.

Stan Development Team. (2020). RStan: The R interface to Stan. R package version 2.21-8.

Yan, W., Hunt, L., Sheng, Q., & Szlavnics, Z. (2000). Cultivar evaluation and mega-environment investigation based on the GGE biplot. Crop Science, 40, 597-605.