Defect physics of Eu-doped GaN is investigated using first-principles hybrid density-functional defect calculations. This includes the interaction between the rare-earth dopant and native defects (Ga and N vacancies) and other impurities (O, Si, C, H, and Mg) unintentionally present or intentionally incorporated into the host material. While the trivalent Eu 3+ ion is often found to be predominant when Eu is incorporated at the Ga site in wurtzite GaN, the divalent Eu 2+ is also stable and found to be predominant in a small range of Fermi-level values in the band-gap region. The Eu 2+ /Eu 3+ ratio can be tuned by tuning the position of Fermi level and through defect association. We find co-doping with oxygen can facilitate the incorporation of Eu into the lattice. The unassociated EuGa is an electrically and optically active defect center and its behavior is profoundly impacted by local defect-defect interaction. Among the Eu-related defects, we identify complexes such as EuGa-ON, EuGa-SiGa, EuGa-Hi, EuGa-MgGa, and EuGa-ON-MgGa as efficient defect-related Eu 3+ centers for non-resonant excitation. This work calls for a re-assessment of certain assumptions regarding specific defect configurations previously made for Eu-doped GaN and further investigation into the origin of the photoluminescence hysteresis observed in (Eu,Mg)-doped samples.