The stem rust effector protein AvrSr50 escapes Sr50 recognition by a substitution in a single surface‐exposed residue

Summary Pathogen effectors are crucial players during plant colonisation and infection. Plant resistance mostly relies on effector recognition to activate defence responses. Understanding how effector proteins escape from plant surveillance is important for plant breeding and resistance deployment. Here we examined the role of genetic diversity of the stem rust (Puccinia graminis f. sp. tritici (Pgt)) AvrSr50 gene in determining recognition by the corresponding wheat Sr50 resistance gene. We solved the crystal structure of a natural variant of AvrSr50 and used site‐directed mutagenesis and transient expression assays to dissect the molecular mechanisms explaining gain of virulence. We report that AvrSr50 can escape recognition by Sr50 through different mechanisms including DNA insertion, stop codon loss or by amino‐acid variation involving a single substitution of the AvrSr50 surface‐exposed residue Q121. We also report structural homology of AvrSr50 to cupin superfamily members and carbohydrate‐binding modules indicating a potential role in binding sugar moieties. This study identifies key polymorphic sites present in AvrSr50 alleles from natural stem rust populations that play important roles to escape from Sr50 recognition. This constitutes an important step to better understand Pgt effector evolution and to monitor AvrSr50 variants in natural rust populations.

Introduction
During a lifetime, plants may encounter many diverse microbes that can affect their growth and development in beneficial or detrimental ways. To maintain their health, plants have evolved an immune system that recognises pathogenic microbes to then trigger an arsenal of defence responses (Jones & Dangl, 2006). Plant immunity relies on extracellular and intracellular receptor proteins that work together to sense either conserved or specific virulence factors in the pathogen and to efficiently activate disease-resistance responses Tang et al., 2017;Ngou et al., 2021). Nucleotide-binding leucine-rich repeat receptors (NLRs) are largely responsible for mediating the intracellular surveillance and constitute the most abundant family of diseaseresistance proteins identified so far (Jacob et al., 2013).
Nucleotide-binding leucine-rich repeat receptors are observed throughout aquatic and land plants and appear to have originated and rapidly diverged before the split of green algae (Gao et al., 2018;Ortiz & Dodds, 2018;Shao et al., 2019).
Nucleotide-binding leucine-rich repeat receptors share a common modular structure with a central nucleotide-binding site (NBS) and C-terminal leucine-rich repeat (LRR) domains, but are separated into three major groups depending on the presence of a Toll-interleukin 1 receptor (TIR), a coiled-coil (CC) or a RPW8/HeLo-like domain at the N-terminus of the proteins (Jubic et al., 2019). Nucleotide-binding leucine-rich repeat receptors recognise specific pathogen secreted proteins, called effectors, through two main strategies; direct binding to the effector  or indirectly by detecting either a modification of the effector host target protein, 'guardee', or modifications of a host protein that mimics the effector target, 'decoy' (van der Hoorn & Kamoun, 2008).
Stem rust disease, caused by the fungus Puccinia graminis f. sp. tritici (Pgt), is considered one of the main threats for wheat production and food security around the world. During the last 2 decades, the emergence of new Pgt virulent races, especially Ug99 detected in Uganda (1999) and virulent on most of the wheat cultivars grown in the field, has prompted the study and identification of effective stem rust resistance genes (Singh et al., 2011). Two main classes of genes, race-specific (NLRs; functional from seedling to adult growth stages) or adult plant resistance (APRs; functioning mainly in the adult stage), are typically used by wheat breeders to prevent plant disease (Ellis et al., 2014). With the exception of Sr60, which encodes a protein with two putative kinase domains (Chen et al., 2020), so far, all cloned race-specific stem rust resistance genes encode NLR proteins with a Nterminal CC domain (Periyannan et al., 2013;Saintenac et al., 2013;Mago et al., 2015;Steuernagel et al., 2016;Zhang et al., 2017Zhang et al., , 2021Chen et al., 2018;Arora et al., 2019;Upadhyaya et al., 2021). One approach to better understand how these immune receptors work and how they should be deployed in the field, is the identification of the matching effectors and their mode of recognition.
To date three Pgt avirulence effector proteins, AvrSr50, AvrSr35 and AvrSr27, have been cloned and recognition by the corresponding NLR proteins Sr50, Sr35 and Sr27 has been demonstrated Salcedo et al., 2017;Upadhyaya et al., 2021). The Sr50 resistance gene is derived from rye and is present in wheat as a translocated segment of rye chromosome 1 but has not been widely deployed in cultivated varieties (Mago et al., 2015). This gene is orthologous to Mla (mildew locus a) in barley, which confers resistance to Blumeria graminis f. sp. hordei (Bgh). Allelic MLA proteins also bind directly to sequenceunrelated effectors from the barley powdery mildew pathogen, B. graminis f.sp. hordei (Bgh), indicating that direct effector recognition seems to be common in the MLA-like family of receptors (Saur et al., 2019). Similarly, NLR proteins from flax also recognise flax rust avirulence proteins through direct interaction (Dodds et al., 2006;Catanzariti et al., 2010), suggesting that direct recognition is common in resistance to biotrophic fungi (Saur & H€ uckelhoven, 2021).
Structural analyses of fungal effectors can provide molecular insight into their function and interaction with host proteins, and how their sequence diversification can lead to escape from recognition. To date, the three-dimensional structures of only a few effectors from rust fungi have been determined: AvrL567 (Wang et al., 2007), AvrM (Ve et al., 2013) and AvrP  from the flax rust fungus Melampsora lini, and MLP124266 and MLP124017 from the poplar rust fungus Melampsora larici-populina (de Guillen et al., 2019). For AvrL567 and AvrP, the structures were important for guiding further mutagenesis studies that identified multiple surfaceexposed polymorphic residues involved in the recognition by their corresponding NLRs (Wang et al., 2007;Zhang et al., 2018). Currently no structures have been determined of effectors from the economically important stem rust pathogen Pgt.
Functional and structural characterisation of rust effector proteins is important to uncover the molecular mechanisms underlying effector gain of virulence and detection by NLR proteins. In this study, we determined the crystal structure of a virulent avrSr50 allele. Functional analysis identified a single polymorphism in a surfaced exposed amino-acid residue that enables the virulent avrSr50 variant to escape from Sr50 recognition. We further identified substantial variation in AvrSr50 alleles in a diverse collection of global Pgt isolates, with four virulence alleles arising from different mechanisms including insertion, amino-acid variation and introduction of a destabilising Cterminal extension fusion through stop codon loss. Sr50 confers resistance to wide array of stem rust isolates, including the Ug99 race group, and is an important source of resistance to introduce in commercial wheat crops. Collectively, our results provide valuable information about AvrSr50 diversity and loss of recognition by Sr50, which is essential for preventing the breakdown of resistance to the devastating stem rust fungus in the field.

Vector construction
For protein expression and purification, the cDNA of three AvrSr50 variants, AvrSr50 WT , AvrSr50 RKQQC and avrSr50 QCMJC  from this point forwards referred to as AvrSr50-A1, AvrSr50-C and avrSr50-B6, respectively, were PCR-amplified, excluding their signal peptide as identified by SignalP-5.0 (Almagro Armenteros et al., 2019), and cloned into a modified, Golden Gate-compatible, pOPIN expression vector using Golden Gate cloning (Supporting Information Table S1, Bentham et al., 2021). Primers consisted of a gene-specific sequence flanked by overhangs compatible for Golden Gate cloning (Table S2) and Golden Gate digestion/ligation reactions and cycling were carried out as described previously (Iverson et al., 2016). The resulting constructs contained an N-terminal His 6 -tag followed by a 3C protease cleavage site to facilitate protein purification by affinity chromatography. Plasmid integrity was confirmed using Sanger sequencing.

Heterologous expression of three AvrSr50 variants in Escherichia coli and protein purification
AvrSr50 proteins were expressed in T7 Express lysY competent E. coli (NEB,C3010I). Bacterial cultures were grown in Terrific Broth (TB) medium at 37°C with shaking at 220 RPM until OD 600 was 0.6-0.8, and then induced with 500 µM isopropyl-1thio-b-D-galactopyranoside (IPTG) and incubated for 16 h at 20°C. Cells were harvested by centrifugation and resuspended in a lysis buffer containing 50 mM HEPES pH 8.0, 300 mM NaCl, 1 mM dithiothreitol (DTT), 1 mM phenylmethanesulfonyl fluoride (PMSF) with 1 µg g À1 of cell pellet DNase I (Sigma). Cells were lysed using an Avestin Emulsifex C5 between 5000-10000 psi, followed by clarification of the lysate by centrifugation at 20 000 g for 35 min at 4°C. The clarified lysate was applied to a 5 ml HisTrap FF crude column (GE Healthcare, Chicago, IL, USA) and the column washed with 20 column volumes (c/v) of buffer consisting of 50 mM HEPES pH 8.0, 300 mM NaCl, 30 mM imidazole, and 1 mM DTT. The proteins were eluted using a linear gradient from 30 mM to 250 mM imidazole over 10 c/v. The fractions were visualised using Coomassie-stained SDS-PAGE, and the fractions containing the protein of interest were pooled. The N-terminal His 6 -tag was removed by overnight treatment with His 6 -tagged Human Rhinovirus 3C protease at 4°C in 50 mM HEPES pH 8.0, 300 mM NaCl, and 1 mM DTT. The 3C-cleaved proteins of interest were purified from the 3C cleavage mixture by passing the sample over a 5-ml HisTrap FF crude column pre-equilibrated with 50 mM HEPES pH 8.0, 300 mM NaCl and 30 mM imidazole, with 1 mM DTT. The column flowthrough was concentrated to 1 ml and purified further by sizeexclusion chromatography on a Superdex 75 16/600 column (GE Healthcare) pre-equilibrated with 10 mM HEPES pH 8.0, 150 mM NaCl and 1 mM DTT. The peak fractions were analysed by Coomassie-stained SDS-PAGE and fractions containing the protein of interest were concentrated to c. 30 mg ml À1 in a 10-kDa MWCO Amicon centrifugal concentrator (Merck, Darmstadt, Germany), before flash freezing in liquid nitrogen and storage at À80°C. Protein concentration was determined from absorbance at 280 nm using an extinction coefficient calculated by the Expasy ProtParam tool; https://web.expasy.org/protparam/.

Crystallisation, data collection and structure determination
Initial screening of crystallisation conditions for AvrSr50 variants was performed in MRC2 96-well plates at 20°C using the sitting-drop vapour-diffusion method and commercially available sparse matrix screens. For screening, 400 nl drops, which consisted of 200 nl protein solution equilibrated against 200 nl of reservoir solution, were prepared using a NT8 ® -Drop Setter robot (Formulatrix, Bedford, MA, USA). The drops were monitored and imaged using the Rock Imager system (Formulatrix). No crystals were observed for either AvrSr50-A1 and AvrSr50-C, however, crystals of avrSr50-B6 (33 mg ml À1 ) were observed in several different conditions. Crystals with the best morphology were observed in 3.5 M sodium formate pH 7.5 (Index C1) and 60% Tacsimate pH 7.0 (Index C5). Crystal optimisation was carried out in 24-well hanging-drop vapour-diffusion plate format, with 1 µl protein solution equilibrated against 1 µl of reservoir solution.
Before X-ray data collection, crystals were transferred into a cryoprotectant solution containing the reservoir solution with 20% glycerol. Some crystals were soaked for c. 60 s in reservoir solution supplemented with 1 M sodium bromide. Singlewavelength anomalous diffraction (SAD) and native datasets were collected on the MX1 (Cowieson et al., 2015) and MX2 (Aragão et al., 2018) beamlines at the Australian Synchrotron, respectively (Table S3). The datasets were processed in XDS (Kabsch, 2010) and scaled using AIMLESS (Evans & Murshudov, 2013) in the CCP4 suite (Winn et al., 2011). For SAD phasing, the CRANK2 pipeline was used (Skub ak & Pannu, 2013) in the CCP4 suite. The model was then refined using phenix.refine in the PHENIX package (Afonine et al., 2012), and iterative model building and refinement was carried out in Coot (Emsley et al., 2010). This model was used as a template for MR (McCoy et al., 2007). Automatic model building was carried out with AUTOBUILD (Terwilliger et al., 2007) and completed using phenix.refine (Afonine et al., 2012) and COOT (Davis et al., 2004). Structure validation was carried out using the MolProbity online server (Davis et al., 2004). The coordinates and structure factors were deposited in the Protein Data Bank (Berman et al., 2002) with the accession no. 7MQQ. Structures were visualised using PYMOL, and structural comparisons were carried out using the PDBeFold (Krissinel & Henrick, 2004) or the DALI server (Holm, 2020).

Transient expression and cell death induction assays in Nicotiana benthamiana and wheat protoplasts
Nicotiana benthamiana plants were grown in a growth chamber at 23°C with 16 h : 8 h, light : dark photoperiod. For agroinfiltration in N. benthamiana, pBIN19, pAM-PAT and pAHC17-pUbi vectors containing either AvrSr50 variants, Sr50 or CC-Sr50 were transformed into Agrobacterium strain GV3101_pMP90 (pBIN19 vectors) and GV3103 (pAM-PAT vectors) by electroporation. Single colonies were selected and grown in Luria-Bertani liquid medium containing 25 lg ml À1 of rifampicin, 15 lg ml À1 gentamicin and either 50 lg ml À1 of kanamycin (pBIN19 vectors) or 25 lg ml À1 of carbenicillin (pAM-PAT vectors) at 28°C for 24 h before infiltration. Bacteria were harvested by centrifugation and resuspended in infiltration medium (10 mM MES pH 5.6, 10 mM MgCl 2 and 150 lM acetosyringone). Bacteria culture densities were adjusted to different
Wheat (Triticum aestivum L.) seedlings of Fielder variety were grown in a growth chamber under shade at 25°C with a photoperiod of 16 h : 8 h, light : dark for 7-9 d. Isolation of wheat protoplasts was performed as described in Arndell et al. (2019). Briefly, the epidermis of the wheat leaves was peeled off and the leaf immersed in 0.6 M mannitol for 15 min to induce plasmolysis. Leaves were then transferred to enzyme solution containing (20 mM MES-KOH (pH 5.7), 0.8 M mannitol, 10 mM KCl, 1.5% (w/v) cellulase R10, 0.75% (w/v) macerozyme R10, 10 mM CaCl 2 , 10 mM KCl and 0.1% (w/v) BSA) and incubated in the dark for 3 h with gentle rotation at c. 40 rpm. Then the solution was filtered through a 100-lM cell strainer. Protoplasts were collected by centrifugation at 80 g for 3 min and washed twice with washing solution W5 (2 mM MES-KOH (pH 5.7), 5 mM KCl, 125 mM CaCl 2 , 154 mM NaCl), followed by gravity collection over 40 min on ice and resuspended in MMG solution (4 mM MES-KOH (pH 5.7), 0.4 M mannitol, 15 mM MgCl 2 ). Protoplasts were counted on a haemocytometer and adjusted to a concentration of c. 3.0 9 10 5 cells ml À1 using MMG solution. The luciferase reporter plasmid (12 lg) and the respective effector and resistance gene constructs (15 lg) were co-transfected into 200 ll of wheat protoplasts using PEG transfection solution (40% (w/v) PEG 4000, 0.2 M mannitol, 100 mM CaCl 2 ). After 16 h protoplasts were collected by centrifugation and dissolved in Luciferase Cell Culture Lysis buffer (Promega, catalogue no. E1531). Luciferase activity was measured using the Luciferase Assay System (Promega, catalogue no. E1501).
Interaction assays in yeast-two-hybrid system AD-AvrSr50 variants were co-transformed with the full-length BD-Sr50 vector in the Saccharomyces cerevisiae HF7c strain and grown on synthetic dropout (SD) medium without tryptophan (W) and leucine (L). For interaction assays, three independent colonies for each combination were grown in 5 ml of liquid SD medium (ÀWL) for 16 h at 30°C. Yeast cultures were adjusted to OD 600 = 1.0 and then diluted 1/10, 1/100 and 1/1000 times in sterile conditions. Aliquots (5 ll) of each sample were inoculated onto selective medium lacking tryptophan, leucine and histidine (ÀWLH) for interaction analysis and incubated at 30°C. Yeast growth was evaluated from 72 to 96 h after incubation.

Wheat transformation and stem rust resistance analyses
Stable transgenic wheat plants were generated by the transformation of pVecNeo binary vectors (Table S1) into the stem rust susceptible cultivar Fielder using Agrobacterium tumefaciensmediated transformation system as previously described (Richardson et al., 2014). For stem rust resistance analyses, 12-20 independent lines were generated for each construct and five T0 lines (Table S4) were infected with Pgt isolate 98-1,2,3,5,6. Transgene presence was confirmed by PCR using Sr50-specific primers (Sr50-F3 and Sr50-R3; Table S2).

Confocal microscopy
Wheat epidermal cells were observed on 10-14-d-old T1 plants and samples were prepared in perfluorodecalin. Wheat cell images were acquired using the Online Fingerprinting mode of ZEN 2012 software on a Zeiss 780 confocal microscope using a LD 940 1.1W lens. Reference spectra for Online Fingerprinting were derived from 14-channel (517-641 nm) Lambda mode data. The YFP spectra derived from samples were compared with manufacturer's data to confirm their identity. The validity of these spectra was checked periodically on experimental samples by Spectral Unmixing of Lambda mode data, thereby confirming that the residuals channel, representing fluorescence emission from uncharacterised sources, was consistently of low intensity.

Allele mining and phylogenetic analyses
AvrSr50 allelic sequences from six rust isolates (Pgt21-0, Ug99, QCMJC, RKQQC, 4a, MCCFC, SCCL) were reported previously . Allelic sequences from a further 40 global isolates were assembled from DNA or RNA sequence data (Illumina paired-end reads, European Nucleotide Archive, PRJEB22223). Reads from each isolate (Table S5) were initially mapped at low stringency (similarity fraction 0.88, length fraction 0.7) to the reference AvrSr50 allele from Pgt21-0 using the 'Map Reads to Reference' tool in CLC Genomics Workbench 10 or later (CLCGWB; Qiagen Aarhus) and the mapped reads were extracted. This approach captured reads representing both haplotypes, which were then de novo assembled (CLCGWB tool) at high stringency (similarity fraction 0.98, length fraction 0.95, Global alignment) to build the two allele models for each isolate. In some cases in which the second allele did not assemble fully, unmapped reads from the high stringency assembly were reassembled with lower stringency. Accuracy of haplotype sequences were cross checked by visual inspection of read mapping tracks in the CLCGWB browser. Translated amino-acid sequences were aligned using CLUSTAL and a maximum likelihood tree constructed using MEGA v.X (Kumar et al., 2018).

Sr50 recognises AvrSr50 and induces stem rust resistance in the cytosol of wheat cells
Nucleotide-binding leucine-rich repeat receptors have been shown to confer resistance and recognise effector proteins in different subcellular compartments. To investigate Sr50 distribution in wheat cells, we generated stable wheat transgenic lines expressing Sr50 fused to Yellow Fluorescence Protein-venus (YFPv) under the control of its native promoter (pSr50) (Fig. 1a). T0 transgenic lines showed a typical resistance response to a P.  Table S5). By contrast, the nontransformed Fielder control was highly susceptible, with extensive large pustule development and no chlorosis/necrosis. Confocal microscopy showed that Sr50:YFPv was mainly present in the cytosol, as previously observed for Sr33 (Cesari et al., 2016) (Fig. 1b). To determine whether Sr50 function is confined to the cytosol we generated additional transgenic lines expressing Sr50:YFPv fused to either a Nuclear Localisation Signal (NLS), a Nuclear Export signal (NES) or nonfunctional variant of these signal motifs (nls and nes respectively). T0 transgenic plants expressing Sr50:YFPv:NES, Sr50:YFPv:nes and Sr50:YFPv:nls showed resistance to Pgt similar to those expressing wild-type Sr50:YFPv ( Fig. 1a; Table S5). However, lines expressing the nuclear-targeted Sr50:YFPv:NLS version showed a similar highly susceptible phenotype to the nontransformed Fielder control (Fig. 1a). All lines were fully susceptible to an isolate of wheat leaf rust (Puccinia triticina), indicating that there was no nonspecific resistance (Fig. S1a). Confocal microscopy confirmed that Sr50:YFPv:NES, Sr50:YFPv:nls and Sr50:YFPv:nes proteins were present in the cytosol but not the nucleus and were indistinguishable from each other and the untagged Sr50: YFPv, while Sr50:YFPv:NLS exclusively accumulated in the nuclei (Fig. 1b). All the fusion proteins were detected by immunoblotting of proteins extracted from the transgenic lines (Fig. S1b). AvrSr50 was previously reported to localise in both the cytosol and nucleus of N. benthamiana cells after transient expression . To further investigate whether AvrSr50 recognition by Sr50 is indeed restricted to the cytosol in wheat, we isolated protoplasts from different Sr50-tagged transgenic lines and co-transformed them with constructs expressing AvrSr50 and the luciferase protein. We then evaluated cell death induction by quantifying the reduction of luciferase activity as an indicator of cell death and therefore AvrSr50 recognition. Substantially reduced luciferase accumulation was observed in the protoplasts from transgenic lines expressing Sr50 in the cytosol compared with protoplasts expressing Sr50 in the nucleus, which were similar to the nontransgenic control line (Fig. 1c). Similar levels of luciferase activity were detected between transgenic lines when co-expressed with GFP as a control (Fig. S1c). Together, these results confirm that a cytosolic Sr50 pool is required to (a)  The crystal structure of AvrSr50 reveals a single surfaceexposed residue that controls recognition by Sr50 We previously identified two divergent alleles of AvrSr50 from isolates RKQQC and QCMJC ) that encode proteins with either 9 or 12 amino-acid differences from the originally identified avirulence allele (Fig. 2a). We designated these alleles as avrSr50-B6 (from QCMJC), AvrSr50-C (from RKQCC) and AvrSr50-A1 (wild-type form), based on a phylogenetic analysis of AvrSr50 variants performed in this study (please refer to the following sections). We previously found that AvrSr50-A1 and AvrSr50-C induced cell death when co-transformed with Sr50 in N. benthamiana plants, while avrSr50-B6 did not, suggesting that the unique polymorphic sites in avrSr50-B6 led to a loss of recognition .
To advance our understanding of AvrSr50 recognition by Sr50, we expressed the three allelic variants in E. coli (residues 23-132) and purified the proteins to homogeneity using nickel affinity chromatography and size-exclusion chromatography. All three proteins were used for crystallisation trials, but we were only able to obtain crystals for avrSr50-B6. The structure of avrSr50-B6 was determined using a bromide ion-based SAD approach (Table S3), to a final resolution of 1.15 A. The final model contains a single molecule in the asymmetric unit consisting of residues 23-132 of avrSr50-B6 with three nonnative residues (Gly-Pro-Met) at the N-terminus remaining following cleavage of the His 6 -tag. The overall structure of avrSr50-B6 has seven b-strands arranged in a b-sandwich fold composed of two antiparallel b-sheets, with a C-terminal a-helix comprised of residues 121-132. The first b-sheet was comprised of strands b1, b2, b4, b5 and b7, and the second b-sheet consisted of b3 and b6 (Fig. 2b). The loop region between b2-b3 (residues Thr58 to Arg65) appeared to exist in multiple conformations, and while this region was included in the final model, we have reduced confidence concerning the absolute position of these residues compared with the rest of the structure (further details in Fig. S2). The five polymorphic residues unique to the nonrecognised avrSr50-B6 variant all map to the surface of the protein (Fig. 2c) and could therefore contribute to recognition specificity by affecting the interaction with the Sr50 receptor.
To evaluate the role of these residues in recognition we introduced the corresponding mutations in each of these five positions  Fig. 2 The crystal structure of avrSr50-B6. (a) Multiple sequence alignment of the three AvrSr50 allelic variants, B6, A1, and C generated using ESPript (Robert & Gouet, 2014). The secondary structure elements of avrSr50-B6 are displayed on the alignment. 'T' corresponds to a b-turn. (b) Cartoon representation showing the overall structure of avrSr50-B6, which consists of seven b-strands arranged in two antiparallel b-sheets, in which the first bsheet comprises strands b1, b2, b4, b5 and b7, and the second b-sheet consists of b3 and b6, and a C-terminal a-helix comprised of residues 121-132 (a1). Topology map of avrSr50-B6 in which the arrows represent b-strands and cylinder represents the a-helix. (c) Surface representation of avrSr50-B6 highlighting polymorphic residues (in cyan) unique to avrSr50-B6 when compared with AvrSr50-A1 and AvrSr50-C. Residue numbers are as indicated. individually to AvrSr50-A1 (V27I, I29T, V61E, V93I, Q121K), including one double mutation (V27I/I29T). These mutants were tested in transient expression assays in N. benthamiana leaves and wheat protoplasts. In wheat protoplasts, constructs expressing AvrSr50 variants were co-transformed with Sr50 and luciferase (all driven by the maize pUbq promoter), and cell death was quantified by relative luciferase expression. Coexpression of Sr50 with the wild-type AvrSr50-A1 version showed a strong reduction in luciferase activity compared with co-expression with the virulence allele avrSr50-B6, or expression of Sr50 or AvrSr50 variants alone (Fig. S3a). All the AvrSr50-A1 mutants induced a strong cell death response similar to wild-type AvrSr50-A1, except for the Q121K substitution, which showed elevated luciferase expression, although slightly less than the native B6 variant (Fig. 3a). In N. benthamiana, AvrSr50 variants were transiently co-expressed with Sr50 driven by either the Fig. 3 Recognition of AvrSr50 mutants by Sr50. (a) Cell death induction was estimated in wheat protoplasts by measuring relative luciferase activity when luciferase is co-transformed with pUbq-Sr50 and either pUbq-AvrSr50-A1 or pUbq-avrSr50-B6 variants carrying the indicated amino-acid substitutions.
Values are relative to samples co-transformed with avrSr50-B6 and luciferase constructs, which was set to 1. The plots represent the results for three biological replicates with error bars indicating SE. Differences between samples were assessed by analysis of variance (ANOVA). Samples marked by identical letters in the plots do not differ significantly (P < 0.05) in the Tukey test. (b, c) Cummings estimation plots showing cell death phenotypes observed upon co-expression of the AvrSr50-A1 (b) or avrSr50-B6 (c) variants with Sr50 (pUbq promoter) in N. benthamiana leaves scored on a scale of 0 (no symptoms) to 4 (strong cell death) (Supporting Information raw data Figs S7-S10). The upper panels show the distribution of the observed scores for each variant in which each dot represents the cell death score of an independent infiltration assay and the number of replicates per sample are indicated in the x-axis. The lower panels represent the mean differences for comparisons of variants against a shared control (A1 or B6). Each mean difference is depicted as a dot with 95% confidence intervals indicated by the ends of the vertical error bars and bootstrap sampling distributions indicated by the horizontal graph. (d) Yeast-two-hybrid protein interaction assay in yeast colonies co-transformed with AvrSr50-A1 and avrSr50-B6 variants fused to the GAL4 activation domain (AD) and Sr50 fused to GAL4 binding domain (BD). Yeast colonies with a OD = 600 and two subsequent dilutions of 1/10 and 1/ 100 were spotted on synthetic double dropout (DDO) medium (ÀWL, growth control) and on synthetic triple dropout TDO (ÀWLH, interaction selection) selective medium.

Research
New Phytologist pUbq or p35S promoters. The cell death response induced by each combination was scored on a 0-4 scale and analysed by estimation graphs (Ho et al., 2019). Similar to the protoplast assays, co-expression of the AvrSr50-A1 mutants with Sr50 consistently induced cell death responses similar to the wild-type protein, with the exception of the Q121K variant that showed no recognition (Figs 3b, S4a, S5-S10). As observed previously , expression of Sr50 from the strong p35S promoter resulted in some background cell death due to Sr50 autoactivity, but nevertheless allowed clear detection of the AvrSr50-induced response (Figs S4a, S5-S7). Conversely, the weaker activity of the pUbq promoter in N. benthamiana largely eliminated the background autoactivity of Sr50, but still allowed detection of the AvrSr50-induced responses (Figs 3b, S8-S10). Both constructs gave qualitatively similar results. Altogether, these results indicated that the Q121K substitution was sufficient to prevent recognition of AvrSr50 by the Sr50 receptor.
We also generated the reciprocal substitution (K121Q) in the virulence allele variant avrSr50-B6 to test if this was sufficient to restore recognition. Indeed, this single amino-acid change restored recognition of avrSr50-B6 in both N. benthamiana and wheat protoplast co-expression assays (Figs 3a,c, S4b, S11, S12). By contrast, the reciprocal I27V and E61V mutations in avrSr50-B6 did not restore recognition. Immunoblotting showed that all of the AvrSr50 protein variants were expressed in N. benthamiana (Fig. S11b). Previously, we have shown that AvrSr50-A1 binds to Sr50 in yeast, while the avrSr50-B6 allelic variant failed to interact . Here, we tested all the AvrSr50 mutants that we generated for binding to Sr50 in yeast-two-hybrid assays and observed a complete correlation between interaction in yeast and recognition in planta (Fig. 4b). That is, only the Q121K mutation disrupted this interaction for the wild-type (A1) variant, while the reciprocal K121Q mutation restored the interaction for the B6 variant. All proteins were well expressed in yeast (Fig. S3b). Taken together, these results supported the idea that Sr50 recognition of AvrSr50 occurs through a direct interaction and that the Q121K polymorphism in avrSr50-B6 is the only polymorphic residue that impacts on AvrSr50 recognition.

AvrSr50 shares structural homology to members of the cupin superfamily and carbohydrate-binding module moiety
The molecular function of AvrSr50 is unknown. Utilising the crystal structure of avrSr50-B6 we carried out a comparative structural analysis using the PDBeFold server (Krissinel & Henrick, 2004) and DALI server (Holm, 2020). The top hits from both servers show structural similarity to cupin superfamily members and carbohydrate hydrolases, in particular glucosidases and xylosidases (Tables S6, S7), despite sharing limited sequence identity (< 10%). In both cases, the structural similarity is limited to the b-sandwich fold and does not include the C-terminal helix. The cupin superfamily is one of the most functionally diverse families known, and includes catalytically inactive seed storage proteins, sugar-binding metal-independent epimerases and metal-dependent enzymes. Most of the known enzymatic cupin proteins, including the hit from PDBeFold (the cupin domain from Bordetella pertussis Tohama I (Fig. 4a; Table S6), contain a transition state metal, coordinated by three histidine residues and an acidic amino acid. However, structural and sequence investigation shows that these coordinating residues are not present in AvrSr50 (Fig. 4b), suggesting it is not likely to be enzymatic. The structure of avrSr50-B6 also shares similarity to HP0902 from Helicobacter pylori (the top hit from the DALI server; Table S7), a nonenzymatic sugar-binding cupin, which shares structural homology to Rml-C-like cupins (Sim et al., 2016) (Fig. 4a). While the substrate-binding location of HP0902 and Rml-C is conserved the residues that co-ordinate sugar binding are not, making it difficult to predict if and which sugars AvrSr50 may bind. Structural similarity between avrSr50-B6 and the carbohydrate hydrolases (Tables S6, S7) is limited to the CBM moiety, which are discrete modules in these larger enzymes that typically aid in maintaining proximity of the carbohydrate to the catalytic domain but have no catalytic activity of their own. Further investigations into these hits do not show conservation in the known residues involved in sugar co-ordination in AvrSr50, and further experimentation is required to determine if AvrSr50 is a sugarbinding protein.

Natural AvrSr50 alleles are highly polymorphic
We had previously characterised the genotype of seven Pgt isolates from Australia and the USA as well as Ug99 from Africa at the AvrSr50 locus, and identified six allelic variants . Here, we analysed publicly available genomic DNA sequence data generated from a larger set of 43 stem rust isolates from around the world (Lewis et al., 2018). After mining the sequence data, we were able to reconstruct the genotypes of 38 isolates, while the other five samples did not contain sufficient sequence depth to reassemble this genomic region in full (Table S4). Most isolates (31) displayed heterozygosity for two different sequence variants, while a single allele was identified in seven isolates consistent with either homozygosity, the presence of a deletion in one haplotype, or insufficient sequence data to resolve the second allele. In total, 14 different AvrSr50 alleles were identified among all the Pgt isolates, with a further two alleles identified in two isolates of P. graminis f. sp. avenae (Pga), which infects oat (Figs 5a, S13, S14). A phylogenetic analysis of these variants revealed two main clades in Pgt, designated A and B, with each containing several closely related alleles (Fig. 5a). The A group includes the originally identified AvrSr50 avirulence allele from Pgt21-0, designated A1, as well as four additional alleles. One of these, A2, has two forms that are identical in sequence except that one version contains a 26-kbp insertion (A2 ins ) within the gene and corresponds to the virulence allele detected in Pgt21-0 and Ug99 isolates . The A5 variant is identical in sequence to A2, except that it contains a single nucleotide change in the stop codon that results in translation of a larger protein with a C-terminal extension of 149 aa. This variant occurs in the UK-01 isolate, which is a representative of the Digalu race group (pathotype TKTTF, Clade IV-B, Lewis et al., 2018). The B group includes the virulence allele avrSr50-B6 from isolate QCMJC, as well as five other closely related variants. Two additional divergent variants outside these clades, C (from RKQQC) and D (from Iranian isolate IR-01) were also observed.
To determine the virulence phenotypes of these additional

New Phytologist
AvrSr50 variants, we generated expression constructs and tested these for recognition by Sr50 in N. benthamiana plants. We observed that the AvrSr50-B1, -B2, -B3, -B4 and -B5 variants all induced cell death in N. benthamiana when co-expressed with pUbq-Sr50 (Figs 5b, S15, S16), and therefore are likely to represent avirulence alleles. None of these variants carried the Q121K substitution, which prevents recognition by the B6 virulence allele and appears to be a unique and recently acquired virulence mutation. The V93I and Q121K substitutions are the only polymorphisms unique to the nonrecognised B6 variant, consistent with our finding above that only the position 121 polymorphism affected recognition. Finally, we observed a lack of recognition and cell death induction of the alleles AvrSr50-A5 and AvrSr50-D. For AvrSr50-A5, this correlated with a lack of expression of the fusion protein (Fig. S16b). This suggests that the C-terminal extension due to the stop codon mutation in this allele results in destabilisation of the protein, therefore preventing recognition. The AvrSr50-D variant contains a Q121L substitution, which may explain its lack of recognition, although this variant also includes numerous other unique variants that could also affect recognition. Because these latter two variants occur in rust strains that also contain a recognised allele of AvrSr50, we cannot directly correlate the in planta recognition with the virulence phenotype they confer in the rust fungus.

Discussion
Functional molecular studies of proteins involved in plantpathogen interactions are important for a better understanding of plant disease emergencies. In many cases, these studies have shown that the breakdown of plant resistance is explained by the evolution of pathogen avirulence effector proteins that escape recognition from major resistance genes. As pathogens exhibit shorter generational cycles compared with those of plants, evolution of effector proteins often occurs at higher rates than those of resistance genes (Karasov et al., 2020). A deeper understanding of how pathogens evolve to gain virulence and the genetic diversity of virulence loci is important for surveying the emergence of potential virulent isolates on natural pathogen populations. To address this, we studied the genetic variability of avirulent/virulent AvrSr50 alleles and identified unique polymorphic sites present in the virulent effector protein variant avrSr50-B6. Using structure-guided site-directed mutagenesis in combination with binding and transient expression assays, we then investigated the role of individual amino-acid substitutions found in avrSr50-B6 to escape Sr50 recognition and gain virulence. Here, we observed that one residue substitution (Q121K) in avrSr50-B6 is sufficient to explain the lack of binding and recognition of this resistanceescaping variant. Substitution of this residue in avrSr50-B6 was sufficient to fully restore recognition as detected in wheat protoplasts, N. benthamiana and yeast-two-hybrid assays. The reciprocal substitution in AvrSr50-A1 abolished recognition in these assays, although not completely in wheat protoplasts, which are likely to reflect slight differences in the quantitative outputs of these assays, or possibly a weak recognition of this variant in the Fielder wheat background. Conversely, we observed that the cumulative polymorphic sites present in the avirulent AvrSr50-C and other natural variants such as AvrSr50-B1, -B2, -B3 and -B4 are not sufficient to abolish recognition. Nevertheless, it is likely that other residues in the AvrSr50 protein contributed to the binding surface interacting with Sr50. This picture would be consistent with other studies in which multiple amino acids across the effector protein surface are involved in binding and recognition by the respective resistance protein (Rehmany, 2005;Wang et al., 2007;Catanzariti et al., 2010;Krasileva et al., 2010). An interesting example is AVR-Pik, an effector from the blast fungus Magnaporthe oryzae, which is recognised by the rice NLR pair, Pik1 and Pik2. Pik1 contains an integrated heavy metal associated (HMA) domain that directly interacts with AVR-Pik. Structural and biochemical studies of the effector variants and Pik-HMA domains have shown that polymorphisms in effector proteins and the HMA domain binding interface underpins binding and ultimately recognition specificity. Nevertheless, they also showed that a single polymorphism in AVR-PikC (Ala67Asp) severely reduces Pik-HMA binding (Maqbool et al., 2015) which nicely illustrates both the robustness of effector recognition and the crucial role of key residues in effector binding (Kanzaki et al., 2012;Maqbool et al., 2015;De la Concepcion et al., 2018Maidment et al., 2021). Similar observations have been reported for the interaction of the RGA5-HMA integrated domain and the M. oryzae effector proteins AVR-Pia and AVR1-CO39 Guo et al., 2018).
To date, at least three independent examples of matching R-Avr proteins in model fungal pathosystems; AvrL567/M. lini L, AVR-Pik-Pik and MLAa/AVR A from/flax, M. oryzae/rice and Bgh/barley (Dodds et al., 2006;Maqbool et al., 2015;Saur et al., 2019), respectively, have been characterised by a high genetic diversity of the interacting proteins. In these examples, gain of virulence mainly occurred through polymorphic changes in key surface-exposed amino acids, which led to escape from recognition. Similarly, here we observed that some virulent variants of the AvrSr50 effector protein have also arisen due to changes in surface-exposed residues involved in the interaction with the resistant protein Sr50. Nevertheless, the insertion of a large DNA fragment (26 kbp) resulting in gene knockout is at the origin of the most common virulence allele of AvrSr50 . Gain of virulence through insertion and deletions has also recently been observed in other Pgt effector proteins. For example, there are two main classes of AvrSr35 virulence alleles, one containing a 400-bp MITE insertion (Salcedo et al., 2017) and the other with a 57-kbp insertion detected in Ug99 (Li et al., 2019). AvrSr35 variants containing a few amino-acid variations from the reference avirulence allele have also been identified, but the effect of these substitutions on recognition specificity has not yet been determined. Similarly, several Pgt isolates in the race 21 clonal lineage have evolved virulence for Sr27 in the field through the deletion of the AvrSr27 locus, which contains two adjacent copies of this Avr gene . The apparent prevalence of Pgt gain of virulence through deletions and insertions over the accumulation of amino-acid sequence changes may be explained by the fact that Pgt populations are mainly clonal due to its heteroecious lifestyle. Therefore, virulent isolates are often the result of spontaneous mutations that occur within asexual lineages, rather than the segregation of pre-existing allelic diversity (Figueroa et al., 2020). Sr50 naturally occurs in rye (Mago et al., 2015) and is orthologous to the Mla resistance gene in barley, which is characterised by extensive functional diversification and encodes multiple allelic resistance proteins with different specificities of recognition (Saur et al., 2019). To date, Sr50 wheat introgression lines have not been widely deployed in commercial crops, implying that selection for the AvrSr50 virulence allele in the agricultural context has been limited. Nevertheless, we have observed four different alleles of AvrSr50 that can allow escape from recognition. One of these results from an insertion that disrupts the gene (A2 ins ), one from the loss of the stop codon (A5) and two from amino-acid polymorphisms (B6 and D). Two of these variants (A5 and D) occur in rust isolates (UK01 and IR01 respectively) that are heterozygous and also contain a recognised allele of AvrSr50, so we cannot directly correlate their lack of recognition in planta with a virulence phenotype in the rust fungus. Evolution of the AvrSr50 variants that we report in this study could have arisen under selection on wild relatives carrying resistance genes related to Sr50 or through occasional infection of rye, which is also a host for Pgt. In addition, it cannot be excluded that AvrSr50 alleles have evolved through the accumulation of neutral variation during long-term clonal reproduction of Pgt lineages in the absence of selection for Sr50 virulence. Apart from a single isolate (QCMJC) that was originally isolated from the sexual host barberry , all of the Pgt isolates analysed here contain at least one presumed avirulence allele for AvrSr50, consistent with the lack of reported virulence for Sr50 in Pgt to date. A single race, TTRTF, detected in Georgia in 2014 showed partial virulence on Sr50, but the genotype is not known (Olivera et al., 2019). However, several isolates in this group, including members of the Ug99, Pgt21 and 'Digalu' race groups (UK-01), are heterozygous and contain a nonrecognised variant at this locus suggesting that virulence could arise in these lineages under field selection by a mutation of a single haplotype.
Identifying and understanding the mode of function of fungal effector proteins is largely hindered by their high sequence diversity, the absence of sequence similarity to proteins of known function, and the lack of conserved domains or functional motifs (Lo Presti et al., 2015). The determination of the threedimensional structure of effector proteins has been useful by providing molecular insights into biochemical functions, as well as guiding mutagenesis studies to understand how they are recognised by the plant. Here, we determined the first crystal structure of an effector protein from Pgt. The structure adopts a bsandwich fold, with a C-terminal a-helix in which the residue responsible for escaping recognition from Sr50 localises. Out of interest, and given the recent developments in ab initio protein modelling, we predicted the structure of avrSr50-B6 using ALPHAFOLD2 (Jumper et al., 2021;Mirdita et al., 2021) with multiple sequence alignment (MSA) generated using MMSEQS2 (Steinegger & S€ oding, 2017). The top model of avrSr50-B6 from ALPHAFOLD (Fig. S17a) has a low confidence score (pLDDT of 40) and comparison with the experimentally determined crystal structure showed a root-mean-squared deviation (RMSD) of 2.63 A, across a small region of the protein (53 ca atoms (Fig. S17c). Visual inspection of the superimposed structures suggests that the overall structural elements are present but the topology of the predicted model is largely incorrect (Fig. S17a,b). The accuracy of ALPHAFOLD predictions is known to be dependent on the number of sequences in the MSA (Jumper et al., 2021), and on inspection the MSA used by ALPHAFOLD for avrSr50-B6 prediction was confined to four sequences. These data highlight the potential limitations of ALPHAFOLD for fungal effector prediction when low sequence representation exists for an effector of interest.
An emerging theme amongst fungal effector proteins is the conservation of structural folds among some effector families without detectable sequence similarity and sometimes originating from different taxa (for recent reviews please refer to Franceschetti et al., 2017;Mukhi et al., 2020). Based on structural similarity searches using DALI and PDBeFold servers, we did not identify any fungal effectors with high structural similarity to AvrSr50. Despite this, visual inspection of the structures of AvrSr50 and the RNase-like effector BEC1054 from Blumeria graminis (PDB: 6FMB) (Pennington et al., 2019) do appear similar. We inspected this further and show via pairwise structural alignment (PDBeFold) both the overall structure and topology are distinct (Fig. S18). These analyses demonstrate that the fold of AvrSr50 is novel amongst fungal effectors, perhaps suggestive of a new effector structural family.
Our structural homology searches revealed similarity of avrSr50-B6 to cupin superfamily members and the CBM of carbohydrate-active enzymes, such as glucosidases and xylosidases. To date, carbohydrate-binding activity in several effector proteins have been shown experimentally. For example, Avr4 from Cladosporium fulvum (Hurlburt et al., 2018) and SnTox1 from Parastagonospora nodorum (Liu et al., 2016), are able to bind chitin and protect the invading fungi from the hydrolytic activity of plant chitinases. Other effectors are known to modulate chitin-induced immunity by sequestering chitin to dampen the host defence response, such as the lysin motif (LysM) effectors including Ecp6 (S anchez- Vallet et al., 2013), MpChi from Moniliophthora perniciosa (Fiorin et al., 2018) and Tal6 from Trichoderma atroviride (Romero-Contreras et al., 2019). It will be interesting to investigate further this structural identity between AvrSr50 and carbohydrate-binding modules by determining the ability of AvrSr50 to bind to different classes of carbohydrates, and the biological relevance of this.

Conclusions
The elucidation of the molecular mechanisms underlying the coevolution between effector and resistance proteins is important to inform both surveillance of natural pathogen populations and deployment of resistance genes in the field. A complementary knowledge of genetic diversity of interacting proteins coupled to functional analyses is a steppingstone towards the development of durable resistance in crops. In this study, we used the natural

Research
New Phytologist genetic variability and structural studies of AvrSr50 to investigate how virulence variants of this protein arise that escape from Sr50 recognition. Functional analyses showed that one residue substitution, Q121K, in the coding sequence was sufficient for AvrSr50 to escape recognition and based on the structure this residue was shown to be surface exposed. These results suggest that single changes in avirulent effector proteins can result on the overcoming of resistance based on major resistance genes. Although Sr50 provides resistance to almost all known lineages of Pgt, several important lineages are heterozygous for avirulence at this locus and may be predisposed to mutate to overcome resistance in Sr50 wheat lines. We also determined the structure of AvrSr50, which is the first structure of a stem rust effector, and show that it has a novel fold among fungal effectors but shares structural homology to cupin superfamily members and carbohydrate-binding modules.

Supporting Information
Additional Supporting Information may be found online in the Supporting Information section at the end of the article.                Table S1 List of plasmids used in this study.

Table S2
List of primers used in this study.
Table S3 X-ray data collection, structure solution and refinement statistics for AvrSr50-B6.


Identification and validation of ferroptosis key genes in bone mesenchymal stromal cells of primary osteoporosis based on bioinformatics analysis

Primary osteoporosis has long been underdiagnosed and undertreated. Currently, ferroptosis may be a promising research direction in the prevention and treatment of primary osteoporosis. However, the specific mechanism of ferroptosis in primary osteoporosis remains a mystery. Differentially expressed genes (DEGs) were identified in bone mesenchymal stromal cells (BMSCs) of primary osteoporosis and heathy patients from the GEO databases with the help of bioinformatics analysis. Then, we intersected these DEGs with the ferroptosis dataset and obtained 80 Ferr-DEGs. Several bioinformatics algorithms (PCA, RLE, Limma, BC, MCC, etc.) were adopted to integrate the results. Additionally, we explored the potential functional roles of the Ferr-DEGs via GO and KEGG. Protein–protein interactions (PPI) were used to predict potential interactive networks. Finally, 80 Ferr-DEGs and 5 key Ferr-DEGs were calculated. The 5 key Ferr-DEGs were further verified in the OVX mouse model. In conclusion, through a variety of bioinformatics methods, our research successfully identified 5 key Ferr-DEGs associated with primary osteoporosis and ferroptosis, namely, sirtuin 1(SIRT1), heat shock protein family A (Hsp70) member 5 (HSPA5), mechanistic target of rapamycin kinase (MTOR), hypoxia inducible factor 1 subunit alpha (HIF1A) and beclin 1 (BECN1), which were verified in an animal model.

Introduction
As a complex bone metabolism disorder mainly characterized by bone loss and bone microstructure changes, primary osteoporosis has long been underdiagnosed and undertreated (1)(2)(3)(4). Approximately 40% of postmenopausal white women are affected by osteoporosis. The aging population will exacerbate the medical and socioeconomic effects of osteoporosis. A patient with osteoporosis has a 40% lifetime fracture risk, and most of them occur in the spine, hip or wrist (5). Osteoporotic fractures affect 8.9 million people around the world per year. It is expected that osteoporotic hip fractures among elderly men will increase by 310% and 240% in women by 2050 (6,7). Osteoporosis and its complications have become a worldwide health crisis.
As a recently identified regulated cell death (RCD), ferroptosis differs from apoptosis (8), necrosis (9), autophagy (10), pyroptosis (11), etc. The main characteristics of ferroptosis are iron-dependent lipid peroxidation and reactive oxygen species (ROS) accumulation (12,13). It involves multiple signaling pathways and their regulators (14) and has been extensively studied in the treatment of cancers, such as renal carcinoma and leukemia (15,16) and hepatocellular carcinoma (17). Ferroptosis has been observed to be associated with various noncancer diseases, such as neurological diseases (18, 19), viral infection, ischemia and reperfusion injury, and atherosclerosis (20). Interestingly, a number of studies have suggested that ferroptosis may be a promising research direction in the prevention and treatment of osteoporosis (21)(22)(23).
Bone mesenchymal stromal cells (BMSCs) are pluripotent mesenchymal cells that can differentiate into different lineages, such as osteoblasts and adipocytes (24). BMSCs play an important role in maintaining bone homeostasis. As a result, osteoporosis is always accompanied by the low osteogenic potential of circulating BMSCs (25), the two form a vicious cycle, exacerbating the progression of osteoporosis. The delicate balance between osteogenic and lipogenic differentiation of BMSCs will be disrupted during the onset of osteoporosis (26). Under certain conditions, ferroptosis can also occur during the process of directional differentiation of BMSCs (27). At present, some researchers have also studied the mechanisms related to ferroptosis and several signaling pathways had been confirmed to play a role in the process of ferroptosis related osteoporosis (28). It is the existing basic research that provides inspiration and guidance for this study. Despite this, the specific mechanisms and related signaling pathways of ferroptosis remain mysterious in osteoporosis.
Bioinformatics analysis of ferroptosis genes in primary osteoporosis has not been reported to date. Thus, we identified DEGs in BMSCs of primary osteoporosis and healthy patients from the GEO databases with the help of bioinformatics analysis. By intersecting these DEGs with the ferroptosis dataset, ferroptosis DEGs (Ferr-DEGs) were obtained. After that, protein-protein interactions (PPI) were used to predict potential interactive networks. Then, the key Ferr-DEGs we obtained were further verified by an experimental OVX model. We have discovered the potential functional role of ferroptosis in BMSCs for the first time and have laid a possible foundation for understanding the pathological mechanism of primary osteoporosis.

Microarray data collection
Information on primary osteoporosis patients was downloaded from the Gene Expression Omnibus (GEO; http:// www.ncbi.nlm.nih.gov/geo), a public database that collects gene sequencing results (29). We downloaded the raw data of GSE35958 uploaded by Benisch from GEO (30). Five samples of primary osteoporosis (GSM878104, GSM878105, GSM878106, GSM878107, GSM878108; op group) and four controls (GSM878100, GSM878101, GSM878102, GSM878103; ctrl group) in GSE35958 were collected. A detailed table of patients' parameters can be seen in Supplementary Table S1.

Datasets analysis
Bioconductor was used to download the annotated R 4.0.3 package, and R was used to convert the microarray probes to symbols. When multiple probes corresponded to one gene ID, the average value was taken for analysis. Principal component analysis (PCA) was used to determine the significant difference dimensions with a P value < 0.05 (31). Relative log expression (RLE) was used for quality control (32). We used Limma to identify the DEGs. P-Value <0.05 and |log2-fold change (FC)| >1 were considered statistically significant (33). In addition, the Ferroptosis Database (FerrDb; http://www.zhounan.org/ferrdb) provided us with a dataset containing 259 genes (34). Then, we intersected it with the DEGs of GSE35958 to obtain ferroptosis DEGs (Ferr-DEGs).

Functional enrichment analysis and visualization
Metascape (http://metascape.org) was used to determine gene functions (35). Biological process (BP), cellular component (CC) and molecular function (MF) are the three major components of gene ontology (GO), which are used to describe gene functions and interactions (36).. The Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway Database is an extensive database for mapping pathway annotation results (37,38). GO and KEGG analyses were performed to explore the functional roles of the Ferr-DEGs via Metascape, and Cytoscape (V3.7.2) was used for visualization (39,40). We then chose molecular complex detection (MCODE) to identify the interconnected central genes. Different colors are used to represent each MCODE network and the close interactions between molecules. During the interaction enrichment analysis, the minimum and maximum sizes of the selection network were 3 and 500, respectively.

Protein-protein interaction network analysis
Then, Ferr-DEGs were uploaded to STRING (https://stringdb.org/) to predict and construct PPI networks with a confidence level of >0.4 (41, 42). The node degree of proteins analyzed by STRING was used for ranking. Cytoscape was used to visualize the network. Nodes represent genes, and edges represent connections between genes. Genes were ranked by size according to the betweenness centrality (BC) score analyzed by the CytoNCA plugin. The cytoHubba plugin was used to rank genes by depth of the color corresponding to the weighted score based upon maximal clique centrality (MCC) algorithms (43,44).

Identification of key Ferr-DEGs
Protein and gene expression profiles at the tissue level (high, medium, low, NA) were obtained from the Human Protein Atlas database (HPA, https://www.proteinatlas.org). After calculation, we regarded the top ranked genes of each algorithm as hub genes and then took the trend of tissue expression and intersected most of those algorithms as the 5 most likely key Ferr-DEGs.

OVX model construction and cell culture
We selected twelve eight-week-old female C57BL/6 mice (provided by the Experimental Animal Center of Soochow University) for osteoporosis induction experiments (45). Mice were randomly divided into two groups: the ovariectomy (OVX) and sham groups. Eight weeks after surgery, all the specimens of the left femurs were removed for microCT and histological experiments. The right femurs were removed for BMSCs collection. Blood was collected by eyeballs removing, and the serum was separated by centrifugation at 2,000 r/m for 20 min and frozen at -80°C for future analysis. All procedures and experiments were approved by the Animal Ethics Committee of the First Affiliated Hospital of Soochow University.
BMSCs were obtained from the right femurs according to the protocol (46). Cells were seeded in a 10-cm 2 dish and incubated in a 37°C incubator with 5% CO 2 . BMSCs after three passages were used for subsequent experiments.

Micro-CT analysis
Left femur samples were scanned with a SkyScan 1176 micro-CT (Aartselaar, Belgium) (n = 6/group) after fixation in 10% buffered formalin for 48 h. Three-dimensional (3D) histomorphometric images were constructed by using NRecon software (SkyScan micro-CT, Aartselaar, Belgium). Bone mineral density (BMD), bone volume (BV), trabecular number (Tb. N), bone volume per tissue volume (BV/TV) and trabecular separation (Tb.sp) were used to evaluate bone mass. The region of interest (ROI) began with 100 pieces below the femur growth plate, and 200 pieces (6 µm each) were read per sample.

Lipid peroxidation assay
A Lipid Peroxidation Malondialdehyde (MDA) Assay Kit (Beyotime, China) was used to measure lipid peroxidation levels. Serum can be used directly for MDA detection. Cells were homogenized on ice in cell lysis buffer (Beyotime, China). After homogenization followed by centrifugation at 12,000×g for 10 minutes, the supernatant was collected for MDA assays. A 100 ml sample was added to each vial containing 200 ml of MDA detection working solution. The samples were mixed and incubated at 100°C for 15 mins. The samples were cooled to room temperature in a water bath and centrifuged at 1000×g for 10 minutes at room temperature. Then, 200 ml supernatant from each sample was pipetted into a 96-well plate, and the absorbance was measured at 532 nm with a microplate reader. MDA content was determined according to the standard curve.

Iron assay
Serum and cellular iron levels were measured by an iron Content Assay Kit (Solarbio, China). Serum can be used directly for iron assay. 100ml serum was mixed with 250 ml of working solution and incubated at 100°C for 5 mins. The samples were cooled to room temperature in a water bath, mixed with 62 ml chloroform, and centrifuged at 10000×g for 10 minutes at room temperature. Then, 210 ml supernatant from each sample was pipetted into a 96-well plate, and the absorbance was measured at 520 nm. BMSCs were homogenized on ice in cell lysis buffer (Beyotime, China). After homogenization followed by centrifugation at 12,000×g for 10 minutes, 20 ml supernatant was collected and mixed with 180 ml working buffer at 25°C for 10 min, and then the absorbance was measured at 510 nm. A standard curve was simultaneously generated according to the manufacturer's instructions.

Immunofluorescence staining
We used immunofluorescence (IF) staining to detect the expression of Ferr-DEGs between the two groups. Antigen repair of the sections was performed, and the sections were blocked with horse serum. Finally, the primary antibody was incubated at 4°C overnight, and fluorescent-labeled secondary antibodies (ab150079, abcam, UK) were incubated at room temperature for 1 hour. The localization and protein expression level of Becn1 (A7353, ABclonal, China) were observed by IF staining under an AxioCam HRC microscope (Carl Zeiss, Germany).

Statistical analysis
The quantitative data are presented as the mean ± SD. Student's t test was used to compare the differences between two groups. Each assay condition was performed in triplicate for all quantitative assays. A P value <0.05 was defined as statistically significant.

Data quality control and identification of Ferr-DEGs
A schematic diagram of this study and the main findings is shown in Figure 1. Samples from GEO were divided into two groups [primary osteoporosis (OP) group versus control (Ctrl) group]. PCA and RLE analysis were performed to determine available dimensions and screen correlated samples according to quality control standards and the normalization of raw data. PCA demonstrated significant differences between the op group and the ctrl group (Supplementary Figure S1A). The RLE plots showed that the normalization was acceptable (Supplementary Figure S1B).
We first analyzed the DEGs of GSE35958 by using the "Limma" package and presented it in a volcano plot  Figure S2). Then, we obtained a gene set including 259 genes from FerrDb and intersected them with the DEGs of GSE35958 to identify Ferr-DEGs. Eighty Ferr-DEGs were found and are shown in a heatmap and Venn diagram (Figures 2A, C). The Ferr-DEGs were further divided into driver, suppressor and marker according to FerrDb (Table 1).

Enrichment analysis of Ferr-DEGs
Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway Database were used to analyze the biological classification of Ferr-DEGs by using Meyascape. The respective results of BP, CC, MF and KEGG are shown in detail (Supplementary Figure S3). The results showed that the significantly enriched genes were involved in the cellular response to stress, response to starvation, ferroptosis, autophagy, response to iron ion, etc. To further understand the relationships between different terms, we visualized it as a network plot (Figures 3A, B). Details of the top 20 terms are listed in the table ( Figure 3C).
The MCODE algorithm was applied to explore those highly connected modules, and three of them were found ( Figures 4A-D), The genes in Module 1 were associated with nongenomic actions of 1,25 dihydroxyvitamin D3, transcriptional activation by NFE2 like bZIP transcription factor 2(NFE2L2) in response to phytochemicals and regulation of protein transport. The genes in Module 2 were associated with autophagy and the response to starvation. The genes in Module 3 were related to protein processing and cellular responses to stress ( Figure 4E).

PPI network analysis and construction of Ferr-DEGs
Because one of the 80 genes was not related to other molecules and did not form a molecular network, we obtained a PPI network containing 79 nodes and 348 edges. Node degree

Selection and analysis of key Ferr-DEGs
We took the top gene in each algorithm as hub genes of different algorithms. After calculation, 10 Ferr-DEGs were obtained for further research. According to the literature, BMSCs tend to differentiate into adipogenic cells rather than osteogenic cells in osteoporosis patients. Therefore, protein and gene expression profiles at the tissue level (high, medium, low, NA) were obtained from HPA. We used the differences in their expression levels in adipose and bone marrow tissues as a preliminary reference for key Ferr-DEGs (Supplementary Figure S4). At the protein level, the expression of SIRT1, MAPK8(Mitogen-activated protein kinase 8), NFE2L2, HIF1A and BECN1 in bone marrow was higher than that in adipose tissues. In contrast, the expression of ATG5 was higher in adipose tissues, and the expression level of HSPA5 was the same, while XBP1(X-box binding protein 1) and JUN were not detected in either tissue. At the gene level, the expression of SIRT1, MAPK8, XBP1, HSPA5, HIF1A and BECN1 in bone marrow was higher than that in adipose tissues, and in contrast, the expression of ATG5, NFE2L2, MTOR and JUN was higher in adipose tissues. Finally, we considered the trends of protein level and gene expression in tissues and intersected most of those algorithms as the 5 most likely key Ferr-DEGs, including SIRT1, HSPA5, MTOR, HIF1A and BECN1 (Figures 6A, B).

Validation of the key Ferr-DEGs in vivo
To verify the key Ferr-DEGs in vivo, an OVX mouse model was used for validation. The OVX model was successfully established in C57BL/6 mice 8 weeks after surgery. The microCT results showed that the BMD, BV, BV/TV and Tb. N of the OVX group were decreased, while Tb.sp was increased (Figure 7). Meanwhile, the mRNAs of BMSCs isolated from the femurs were used to detect the expression of those 5 key Ferr-DEGs. The results showed that the relative mRNA expression levels of Sirt1, Hspa5, Mtor, Hif1a and Becn1 in the OVX group were significantly lower than those in the sham group (Figures 8A-E). To further validate the expression of these five differentially expressed genes in vivo, we chose Becn1 as the representative target, and IF staining was used for evaluation. The results showed that BECN1 was significantly more highly expressed in the sham group than in the OVX group (Figures 8F, G).
To further demonstrate the involvement of ferroptosis in the development of osteoporosis, lipid peroxidation levels and iron concentrations in serum and BMSCs from OVX mice were measured. The results showed that MDA levels and iron concentrations were significantly increased both in serum and BMSCs from OVX group (Figures 9A-D). Lipid peroxidation and iron overload are the characteristic manifestations of ferroptosis.
In summary, these in vivo results are consistent with the bioinformatics analysis, demonstrating that there is a high correlation between ferroptosis genes and the progression of osteoporosis in a mouse model. However, the exact mechanism needs to be further elucidated.

Discussion
Different types of osteoporosis, such as diabetic osteoporosis, glucocorticoid-induced osteoporosis, and postmenopausal osteoporosis, were confirmed to be associated with ferroptosis (22, 47-50). However, the specific mechanisms and signaling pathways of ferroptosis involved in osteoporosis remain a mystery. We first intersected microarrays on primary osteoporosis patients in GEO with FerrDb, especially the detection and analysis of BMSCs, which has not been reported in previous studies. We aimed to identify the key genes of ferroptosis and explore its mechanisms in primary osteoporosis. Therefore, we compared the gene expression of primary osteoporosis patients with controls. To be more statistically persuasive, several bioinformatics algorithms (PCA, RLE, Limma, BC, MCC, etc.) were adopted to integrate the results. Additionally, PPI network construction, GO and KEGG were performed to explore the potential functional roles of the Ferr-DEGs. Finally, 80 Ferr-DEGs and 5 key Ferr-DEGs were calculated for further study.
As a newly discovered form of RCD, ferroptosis is accompanied by iron overload, lipid peroxidation and oxidative stress and is regulated by manifold genes. There appears to be a synergistic or antagonistic effect between ferroptosis and other RCDs because they have cross regulators in some pathways (51-53). Interestingly, according to our findings in this research, autophagy was regarded as the most crucial biological process after ferroptosis. Autophagy is an intracellular degradation system that maintains the stability of the intracellular environment (54, 55). Including three subtypes: microautophagy, macroautophagy and chaperonemediated autophagy (CMA) (56)(57)(58). Several stress-related proteins were reported to be important regulators of autophagydependent ferroptosis in cancer therapy (59). New evidence suggests that abnormal levels of autophagy and mitochondrial autophagy disrupt the balance of bone metabolism. Osteoporotic BMSCs exhibit impaired osteogenic and increased adipogenic differentiation due to decreased autophagy (60, 61). Our results indicate that ferroptosis is involved in the process of osteoporosis and that autophagy is also involved. This conclusion was further supported by the MCODE results, which suggested that the Ferr-DEGs in different modules were connected with nongenomic actions of 1,25-dihydroxyvitamin D3, autophagy, and response to starvation and stress. We verified the potential connection between ferroptosis and autophagy in osteoporosis via various bioinformatics methods. The interaction between ferroptosis and autophagy deserves further study.
As a group of cells with the potential for multidirectional differentiation, BMSCs can differentiate into osteoblasts, chondrocytes and adipocytes. There is an opposite relationship between the differentiation of BMSCs into osteoblasts and adipocytes (62, 63). When osteoporosis occurs, BMSCs tend to differentiate into adipocytes rather than osteoblasts (64, 65). To test and verify our hypothesis, based on the HPA database, we compared the protein and gene expression levels of 10 hub genes between human bone marrow and adipose tissue. We considered the trends of protein and gene expression in tissues and intersected most of those algorithms as the 5 most likely key Ferr-DEGs, including SIRT1, HSPA5, MTOR, HIF1A and BECN1. Ultimately, these 5 key Ferr-DEGs were considered to be a preliminary reference for determining whether they are differentially expressed when osteoporosis occurs.
SIRT1 is an important regulator of bone homeostasis. The results of HPA showed that the expression level of SIRT1 in adipose tissue is NA, which is not entirely consistent with what is known thus far as SIRT1 is a potential regulator in adipose tissue inflammation and metabolism. The conflicting results may be due A B FIGURE 6 The expression analysis of hub genes in human tissues according to Human Protein Atlas data base. to delays in updating the site. However, the relative expression levels of SIRT1 in adipose tissue and bone marrow tissue were consistent, which did not affect the conclusion of this study. Sirt1 KO mice presented with osteoporosis characterized by decreased osteogenesis and increased adipogenesis in BMSCs (66). In contrast, BMSCs overexpressing Sirt1 repressed the increased expression of superoxide dismutase 2 (Sod2) and forkhead box O3A (Foxo3a), which promoted the differentiation of BMSCs into osteoblasts and delayed senescence (67). A study in head and neck cancer cells showed that silencing SIRT1 expression can inhibit epithelial-mesenchymal transition and decrease ferroptosis, while SIRT1 agonists can promote ferroptosis (68). SIRT1 was also reported to be associated with mitophagy (69). HSPA5 is a chaperone protein mainly expressed in the Endoplasmic reticulum (ER) that maintains the stability of glutathione peroxidase 4 (GPX4) by forming the HSPA5-GPX4 complex, thereby causing resistance to ferroptosis (70, 71). There seems to be no current research on the potential role of HSPA5 in osteoporosis, which would be an interesting direction. HIF1A degradation, circadian rhythm and lipid peroxidation play a role in the regulation of ferroptosis, and HIF1A is a key factor (72).   cell line was significantly reduced after knockdown of autophagy related 7 (Atg7) and Becn1 (77). Our previous study also found that tet methylcytosine dioxygenase 2(Tet2) promoted bone loss in OVX mice. It positively regulates BECN1-dependent autophagy by inhibiting BCL2 expression and promoting osteoclast differentiation (78). However, whether BECN1 is involved in ferroptosis during BMSCs differentiation is still unclear and will be further researched. Acts as a central regulator of the cellular response to growth stimulation, MTOR can inhibit autophagydependent ferroptosis. In cancer cells, interactions between MTOR and GPX4 can regulate autophagy-dependent ferroptosis by inhibiting CMA (79), and MTOR inhibitors may promote GPX4 degradation by activating the CMA pathway. MTOR also has a protective effect on ferroptosis in BMSCs, and autophagy mediated by the MTOR pathway can regulate the regeneration function of BMSCs, thereby influencing the occurrence and development of osteoporosis in postmenopausal women (61). Collectively, there are few studies about the specific functions and related mechanisms of these 5 key Ferr-DEGs in ferroptosis during osteoporosis progression. Our results provide a very meaningful direction for future research.
To further verify our conclusions, an osteoporotic model was successfully established in OVX mice. We found that the expression of Sirt1, Hspa5, Mtor, Hif1a and Becn1 in BMSCs from the OVX group was lower than that in BMSCs from the control group, which is consistent with our bioinformatic analysis results. The results of immunofluorescence staining were also consistent with our expectations. Detection of lipid peroxidation and iron levels at serum and cellular levels also confirmed that ferroptosis may be involved in the development of primary osteoporosis. These outcomes suggest that ferroptosis of BMSCs may be involved in the pathological process of primary osteoporosis. Finally, we believe that our results will offer new insight into the role of ferroptosis in primary osteoporosis and identify the 5 key Ferr-DEGs as potential biomarkers for primary osteoporosis diagnosis and treatment.

Conclusion
In summary, through a variety of bioinformatics methods, our research successfully identified 5 key Ferr-DEGs associated with primary osteoporosis and ferroptosis, namely, SIRT1, HSPA5, MTOR, HIF1A and BECN1. In addition, autophagy may also be involved in ferroptosis-related primary osteoporosis. In the future, we intend to collect enough bone tissues from patients with primary osteoporosis for further research. We aimed to clarify the regulatory roles of these genes to explore the differentiation mechanisms of BMSCs.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics statement
The animal study was reviewed and approved by Animal Ethics Committee of the First Affiliated Hospital of Soochow University.


s2fieldsofstudy: Biology Medicine

text:
Comprehensive Evaluation of Anti-PD-1, Anti-PD-L1, Anti-CTLA-4 and Their Combined Immunotherapy in Clinical Trials: A Systematic Review and Meta-analysis

Immunotherapy with immune checkpoint inhibitor (ICI) drugs is gradually becoming a hot topic in cancer treatment. To comprehensively evaluate the safety and efficacy of ICI drugs, we employed the Bayesian model and conducted a network meta-analysis in terms of progression-free survival (PFS), overall survival (OS) and severe adverse events (AEs). Our study found that treatment with ipilimumab was significantly worse than standard therapies in terms of PFS, whereas treatment with cemiplimab significantly improved PFS. The results also indicated that cemiplimab was the best choice for PFS. Treatment with nivolumab, pembrolizumab and nivolumab plus ipilimumab significantly improved OS compared to standard therapies. In terms of OS, cemiplimab was found to be the best choice, whereas avelumab was the worst. In terms of severe AEs, atezolizumab, avelumab, durvalumab, nivolumab, and pembrolizumab all significantly reduced the risk of grade 3 or higher AEs compared to standard therapy. The least likely to be associated with severe AEs were as follows: cemiplimab, avelumab, nivolumab, atezolizumab, and camrelizumab, with nivolumab plus ipilimumab to be the worst. Therefore, different ICI drug therapies may pose different risks in terms of PFS, OS and severe AEs. Our study may provide new insights and strategies for the clinical practice of ICI drugs.

INTRODUCTION
Immunotherapy has become one of the most important breakthroughs in the treatment of cancer in recent years, and its development has promoted changes in many cancer treatment methods. As a series of co-inhibitory and co-stimulatory receptors and ligands, immune checkpoint inhibitors (ICI) drugs can block negative regulatory factors expressed by immune or tumor cells to enhance their immune function against cancer cells, mainly programmed death-1 (PD-1), programmed death-ligand-1 (PD-L1) and cytotoxic T lymphocyte antigen-4 (CTLA-4) (Rosenberg et al., 2004). In 2011, the CTLA-4 inhibitor ipilimumab was approved by the US Food and Drug Administration for the treatment of advanced melanoma (Hodi et al., 2010). Subsequently, several ICI drugs were also approved for the treatment of cancer (Topalian et al., 2012;Gong et al., 2018). Since then, interest for immunotherapy with ICI drugs has been increasing. Many studies focused on the prognosis and treatment for different cancers (Wu et al., 2015).
Chemotherapy is the first-line treatment for advanced cancer, and patients undergoing chemotherapy often experience severe adverse events (AEs). Although ICI drugs have achieved good anticancer effects in the treatment of many solid tumors, they may still cause severe treatment-related or drug-related AEs. Progression-free survival (PFS) and overall survival (OS) are usually efficacy end-points. In terms of PFS and OS, the therapeutic effects of ICI drugs remain unclear compared with standard therapies. Due to the limitations of randomized clinical trials, the overall safety evaluation of different ICI drugs for cancer treatment is not comprehensive, especially in terms of PFS and OS.
We conducted a systematic review and network meta-analysis of the therapeutic effects of ICI drugs targeting PD-1, PD-L1, and CTLA-4, focusing on PFS, OS and treatment-related severe AEs in patients receiving ICI drug monotherapy, combination therapies and standard therapies (chemotherapy, targeted therapies and their combination therapies included). This study comprehensively evaluated the safety and efficacy of different ICI drugs and their combination therapies, aiming to provide better guidance for the clinical application of various ICI drugs.

Search Methods and Study Selection
We searched PubMed, Embase, and Cochrane Library for English-language studies between January 2000 and September 2021, using keywords such as ipilimumab, tremelimumab, pembrolizumab, nivolumab, cemiplimab, camrelizumab, toripalimab, tislelizumab, spartalizumab, atezolizumab, avelumab, durvalumab, PD-1, PD-L1, and CTLA-4. The search strategy was described in Supplementary Table S1. The study search, selection and data extraction were independently conducted by two reviewers (ZX and ZZ), and discrepancies were evaluated by an independent reviewer (JL). The three authors (ZX, JL and ZZ) reviewed and discussed the full text of studies that may be eligible, and differences of opinions were resolved by consensus.
Only high-quality head-to-head phase 2 and 3 randomized controlled trials (RCTs) comparing two or more treatments including ICI drug monotherapy, ICI drug combination therapies and standard therapies were included. Some RCTs only presented interim results, as insufficient information may affect the final analysis, we selected the most recent results as much as possible. Data provided include at least one of the following: hazard ratios (HRs) of PFS, OS and treatment-related severe AEs. We excluded reviews, conference abstracts and posters. The study was performed based on the Preferred Reporting Items for Systematic Reviews and Meta-analyses (PRISMA) guideline (Hutton et al., 2015;Wang et al., 2021). This study was approved by International prospective register of systematic reviews (PROSPERO) (registered ID: CRD42021278158).

Data Extraction
The authors (ZX and ZZ) independently extracted data according to the PRISMA guidelines. The first author, year of publication, national clinical trial identification number, trial name, phase, number of patients, type of cancer, drug used, follow-up time, number of severe AEs, HRs, and confidence interval (CI) of PFS and OS were summarized in standardized Tables 1-3. PFS is considered to be the primary endpoint of randomized clinical trials evaluating patients with solid tumors (Korn and Crowley, 2013). OS is defined as the time from the start of treatment to death or the last follow-up. The HRs of PFS and OS represent HRs between treatment 1 and 2. In assessing AEs, we chose treatment-related or drug-related AEs as the main results. If there were no treatment or drug-related AEs in studies, we included all AEs. The classification of AEs is often used to evaluate the type and severity of AEs in clinical trials. According to AE classification, grade 3 or higher AEs are considered as severe AEs. The risk of severe AEs is the focus of the evaluation of therapeutic effectiveness, so the number of AEs surveyed and severe AEs were both extracted (Xu et al., 2018).

Adverse Events Analysis
We used gemtc and pcnetmeta packages in R v4.0.3 and called JAGS v4.3.0 to perform statistical analysis in a Bayesian framework based on Markov Chain Monte Carlo (MCMC) methods, and generated the graph depicting the network geometry .
Firstly, we made a rough comparison between the fit of the consistent model with the inconsistent model. Secondly, the inconsistency on the specific comparison was tested by node splitting analysis. p < 0.05 was considered as indicating a significant inconsistency. Outstanding consistency is the key to robust results, as evidenced by the consistency between direct and indirect results. We compared the results of network metaanalysis (indirect results) with those of pairwise analysis (direct results) to explore the sources of inconsistency. Additionally, if there existed significant heterogeneity, we used the random-effect model. Otherwise, we used the fixed-effect model (Dias et al., 2011). We used non-information prior distributions and overdispersed initial values (scaling 2.5) in 3 chains to fit the model. 56 independent randomized controlled         experiments yielded 100,000 iterations (including 20,000 optimization iterations) with 10 refinement intervals for each chain. This method was used to generate a posterior distribution of model parameters. The convergence of iterations was evaluated by using the Gelman-Rubin-Brooks statistics, all of which converge near 1. Based on the odds ratio (OR) advantage ratio and posterior probability, we ranked probabilities of each treatment as the safest, followed by the second, third and so on.

Progression-free Survival and Overall Survival Analysis
For the consistency and heterogeneity analysis of PFS and OS, we chose to use R's netmeta package in the Frequentist framework to make a preliminary judgment using the traditional frequency method, avoiding the artificial bias caused by complex prior settings, settings of dummy variables and variance-covariance matrices of regression models in Bayesian statistics, which would simplify the operator's parameter setting. The I 2 test was used to evaluate the heterogeneity between studies, with the significance level set as p < 0.05. I 2 greater than 25, 50 or 75% indicated low, medium and high heterogeneity respectively. If significant heterogeneity exists, the random-effect model was used. Otherwise, we employed the fixed-effect model (Higgins and Thompson, 2002).
Since Bayesian statistics are more accurate and the results are highly consistent with those in the frequency model, we subsequently chose the Bayesian framework by using the MCMC method in WinBUGS v1.4.3 for network metaanalysis. We used the consistency model (due to I 2 < 25) to calculate HRs and 95% CIs. We simulated 3 different chains, each with 45,000 built-in samples, resulting in 15 iterations with a refinement rate of 15 (3 different chains with 15,000 iterations and 45,000 burn-in samples and 50 thinning rates). The model fitting was further determined according to the deviation information criterion. The output was a posterior distribution of relative effect size, and we got the estimated average of HR and 95% CI (95% CI as the 2.5th and 97.5th percentiles) (van Valkenhoef et al., 2012). The ranking probability distribution was calculated, ranking the probabilities of each treatment as the safest, followed by the second, third and so on.

Progression-free Survival
In analyzing PFS, no significant heterogeneity (I 2 = 19%) or inconsistency was observed (p = 0.97) ( Supplementary Table  S2). Therefore, the Bayesian fixed-effect model was used. HRs and 95% CI from the network meta-analysis are shown in Figure 3A. Treatment with ipilimumab was significantly worse in terms of PFS than standard therapies, whereas treatment with cemiplimab significantly improved PFS. According to the probability ranking diagram, the results showed that cemiplimab was the best choice in terms of PFS, camrelizumab ranked the second safest and nivolumab plus ipilimumab ranked the third safest, whereas ipilimumab was the worst ( Figure 3B). Additionally, treatment with ipilimumab was significantly worse than most other treatments in terms of PFS. Interestingly, nivolumab plus ipilimumab significantly improved PFS compared to ipilimumab, which suggested that treatment with combinations of ICI drugs may benefit PFS compared to monotherapy. The results calculated according to the frequency model were highly consistent with the results of the Bayesian fixed-effect model. Additionally, we performed subgroup analyses based on treatment of different cancer types, particularly lung cancer and melanoma. The safety profile and probability ranking diagram for lung cancer and melanoma are shown in Supplementary Figures S1, S4 in the Supplement respectively. Cemiplimab was also the best choice in terms of PFS in treating lung cancer, and nivolumab plus ipilimumab ranked the second safest. Compared with standard therapies, HR (95% CI) for cemiplimab was 0.77 (0.61-0.96). Treatment with cemiplimab also significantly improved PFS compared to nivolumab.
Tremelimumab was considered the worst choice in terms of PFS in treating lung cancer. In terms of melanoma, our results showed that nivolumab plus ipilimumab was the best choice for PFS. HR (95% CI) for nivolumab plus ipilimumab was 0.69 (0.52-0.92) compared with standard therapies. In addition, HRs (95% CI) for nivolumab and pembrolizumab were 0.78 (0.65-0.94) and 0.80 (0.64-0.99) respectively. The probability ranking diagram of melanoma indicated that ipilimumab was the worst choice for PFS.

Overall Survival
In analyzing OS, no consistency (I 2 = 0%) or inconsistency (p = 0.60) ( Supplementary Table S2) was observed, and so the Bayesian fixed-effects model was applied. HRs and 95% CI are  shown in Figure 4A. Treatment with nivolumab, pembrolizumab and nivolumab plus ipilimumab significantly improved OS compared to standard therapies. According to the probability ranking diagram, cemiplimab was the best choice in terms of OS, and durvalumab ranked the second safest, whereas avelumab was the worst ( Figure 4B). Of note, nivolumab plus ipilimumab may improve OS compared with nivolumab and ipilimumab monotherapy, which was similar to durvalumab plus tremelimumab compared with durvalumab and tremelimumab monotherapy. The results calculated based on the frequency model were also highly similar to the results of the Bayesian fixed-effect model.
We also conducted subgroup analyses for lung cancer and melanoma. The safety profiles and probability ranking diagrams of lung cancer and melanoma are shown in Supplementary  Figures S5, S8 in the Supplement. For lung cancer treatment, cemiplimab was the best option for OS and durvalumab ranked the second safest. Safety profile of lung cancer suggested that compared with standard therapies, HRs (95% CI) for nivolumab, nivolumab plus ipilimumab and pembrolizumab were 0.91 (0.82-0.99), 0.87 (0.76-0.99), and 0.89 (0.80-0.99) respectively. Standard therapies were considered to be the worst option for lung cancer in terms of OS whereas. For melanoma treatment, nivolumab plus ipilimumab was the best option. Safety profiles showed that HR (95% CI) for nivolumab was 0.84 (0.71-0.99) compared with standard therapies. Our results also indicated that standard therapies were the worst choice for melanoma in terms of OS.

Severe Adverse Events
In the network meta-analyses of severe AEs, high heterogeneity was found (Supplememntary Table S3), and the random-effect model was employed. Safety profile in the consistency model is shown in Figure 5A. Atezolizumab, avelumab, durvalumab, nivolumab, and pembrolizumab all significantly reduced the risk of grade 3 or higher AEs compared to standard therapies.
Through node splitting analysis, significant inconsistency could not be detected for most comparisons (Supplementary Table S4). There was significant inconsistency between ipilimumab and nivolumab plus ipilimumab, and ipilimumab and standard therapies (p < 0.05). The comparison between nivolumab plus ipilimumab and standard therapies also showed a degree of inconsistency (p = 0.08). In the direct comparisons, patients receiving nivolumab plus ipilimumab were more likely to have severe AEs than those receiving ipilimumab, and patients receiving ipilimumab were more likely to have severe AEs than those receiving standard therapies. However, patients receiving standard therapies were more likely to have severe AEs than those receiving nivolumab plus ipilimumab. The comparison between the above three groups may be the main reason for the inconsistency. In the safety profile, efficacy of treatment for grade 3-5 adverse events is represented as odd ratios (ORs) with 95% confidence intervals. All comparisons are made as column versus row. Statistically significant results are in bold. Probability ranking diagram shows the probability of the safety of different therapies ranking the first to the last for severe adverse events.
Frontiers in Pharmacology | www.frontiersin.org We performed subgroup analyses of lung cancer, melanoma, and squamous cell carcinoma treatment. The respective safety profiles and probability ranking diagram are shown in Supplementary Figures S9, S14 in the Supplement. Interestingly, the results suggested that nivolumab was the joint best choice for lung cancer, melanoma and squamous cell carcinoma. Standard therapies, based on the probability ranking diagram, were considered to be the worst for lung cancer and squamous cell carcinoma treatment, and nivolumab plus ipilimumab was the worst for melanoma treatment.

DISCUSSION
In order to comprehensively evaluate the safety and efficacy of ICI drug monotherapy and combination therapies, we conducted a network meta-analysis combining HRs of PFS and OS, and the risk of severe AEs, and performed subgroup analyses particularly for lung cancer and melanoma. The application of bioinformatics is often used to analyze published data (Wu et al., 2017). To our knowledge, this study is the first comprehensive report comparing PFS HRs, OS HRs, and corresponding treatmentrelated severe AEs among ICI drug monotherapy, combination therapies and standard therapies.
In terms of PFS and OS, we first tested the heterogeneity and consistency of network meta-analysis based on the frequency method. No significant heterogeneity and consistency were found, indicating that this network meta-analysis was consistent in PFS and OS. We used the frequency model and the Bayesian model separately. The results of the frequency model and the Bayesian model agree well. In view of the greater accuracy of the Bayesian model, our final results were presented by the Bayesian model.
In terms of PFS, treatment with ipilimumab was significantly worse than standard therapies, whereas treatment with cemiplimab significantly improved PFS. The results also indicated that cemiplimab was the best choice for PFS. Treatment with nivolumab, pembrolizumab and nivolumab plus ipilimumab significantly improved OS compared to standard therapies. For OS, cemiplimab was considered to be the best choice, whereas avelumab was the worst. Since few studies compared cemiplimab and camrelizumab with other therapies, nivolumab plus ipilimumab ranked the third safest in PFS and durvalumab ranked the second safest in terms of OS. In comparing ICI drug combination therapies with monotherapy, we found that nivolumab plus ipilimumab significantly improved PFS compared to ipilimumab. Additionally, nivolumab plus ipilimumab may improve OS compared with nivolumab and ipilimumab monotherapy, similar to durvalumab plus tremelimumab compared with durvalumab, and tremelimumab monotherapy. However, further studies are needed to compare the safety and efficacy of ICI drug combination treatment with monotherapy. Although abundant included studies may lead to low significance of the overall results in terms of PFS and OS, we thought that the results would be more reliable and further performed subgroup analysis.
Severe AEs were examined as the measure of the toxicity of different ICI drug therapies and standard therapies. In terms of severe AEs, there was a large inconsistency in the comparison between the above three groups, which is similar to inconsistency reported by Chen et al (Xu et al., 2018), but the degree of inconsistency was more obvious in the comparison between ipilimumab and standard therapies in this study. We considered the main reasons as follows: In spite that the inclusion and exclusion criteria in our study are similar to those of Chen et al, eligible studies with inconsistent results were relatively abundant, which may lead to higher inconsistency. Of note, Chen et al combined durvalumab plus tremelimumab and nivolumab plus ipilimumab as the two ICI drug group, however, we grouped them in our study to assess the safety and efficacy more precisely.

Strengths and Limitations
The main strengths of our studies are as follows: we used the Bayesian model to conduct network meta-analysis and then employed the frequency model for inconsistency test and result verification in terms of PFS and OS. We found that the results of the frequency model were highly consistent with those of the Bayesian model, and we represented our final results from the Bayesian model, which greatly enhanced the reliability of conclusions. To comprehensively investigate the safety and efficacy of various ICI drugs, we analyzed three different indicators PFS, OS and severe AEs. Of note, we included enough studies to ensure the accuracy of the results. Despite that some studies represent the data of the same RCTs at different times, we chose the most recent results as much as possible.
This study also has several limitations. Firstly, enrolled studies showed high heterogeneity. In order to avoid publication bias, we tested the heterogeneity and used different models accordingly. Secondly, the number of RCTs that meet the requirements for inclusion is different among ICI drugs at present, and there is obvious inconsistency in severe AEs, which require more studies for higher-level verification. Thirdly, despite randomization of the eligible studies, there are still characteristic imbalances between the groups in trials.

CONCLUSION
In the present study, the Bayesian model was used to comprehensively assess survival data and the risk of severe AEs for ICI drugs, which showed that different ICI drug therapies may pose different risks in terms of PFS, OS and severe AEs. Our study may provide new insights and strategies for the clinical practice of ICI drugs.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.