If very massive stars (M ≳ 100 M⊙) can form and avoid too strong mass-loss during their evolution, they are predicted to explode as pair-instability supernovae (PISNe). One critical test for candidate events is whether their nucleosynthesis yields and internal ejecta structure, being revealed through nebular-phase spectra at t ≳ 1 yr, match those of model predictions. Here, we compute theoretical spectra based on model PISN ejecta at 1–3 yr post-explosion to allow quantitative comparison with observations. The high column densities of PISNe lead to complete gamma-ray trapping for t ≳ 2 yr which, combined with fulfilled conditions of steady state, leads to bolometric supernova luminosities matching the 56Co decay. Most of the gamma-rays are absorbed by the deep-lying iron and silicon/sulphur layers. The ionization balance shows a predominantly neutral gas state, which leads to emission lines of Fe i, Si i, and S i. For low-mass PISNe, the metal core expands slowly enough to produce a forest of distinct lines, whereas high-mass PISNe expand faster and produce more featureless spectra. Line blocking is complete below ∼5000 Å for several years, and the model spectra are red. The strongest line is typically Ca ii λλ7291, 7323, one of few lines from ionized species. We compare our models with proposed PISN candidates SN 2007bi and PTF12dam, finding discrepancies for several key observables and thus no support for a PISN interpretation. We discuss distinct spectral features predicted by the models, and the possibility of detecting pair-instability explosions among non-superluminous supernovae.