We perform a global sensitivity analysis with two coupled codes that model the interaction between the plasma boundary and the material surface in fusion energy applications. The coupled codes depend on 19 parameters, including uncertain operating conditions, uncertain material properties and uncertain numerical parameters. We construct a set of generalized polynomial chaos surrogate models that allow us to extract the desired Sobol' sensitivity indices as a post-processing step. We illustrate how the predicted sensitivity indices change when considering both codes in isolation, or when considering the coupled setting. Finally, we investigate if multifidelity polynomial chaos methods can be used to alleviate the computational burden in the coupled setting.