The central topic of this work are masses and mixing parameters of the η-η' system, which are investigated within the framework of Wilson twisted mass lattice QCD, using gauge configurations provided by the European Twisted Mass Collaboration. We present the first calculation with Nf=2+1+1 dynamical quark flavors performed at three different values of the lattice spacing and multiple values of the light quark mass, corresponding to charged pion masses ranging from ∼ 230 MeV to ∼ 500 MeV. Moreover, we use selected ensembles which differ only by the value of the strange quark mass while all other parameters are kept fixed in order to obtain information on the strange quark mass dependence of our observables. This allows us to carry out chiral and continuum extrapolations with well-controlled systematics for the mass of the η meson. Using the standard method, the statistical error for the η' turns out significantly larger due to the large contributions of quark disconnected diagrams and autocorrelation effects. However, employing an improved analysis method based on an excited state subtraction in the connected pieces of the correlation function matrix, it becomes feasible to obtain a result for the η' mass with controlled systematics as well. The values for both masses Mη=551(8)stat(6)sys MeV and Mη' = 1006(54)stat(38)sys(+64)ex MeV turn out to be in excellent agreement with experiment. Considering matrix elements in the quark-flavor basis, one expects the mixing in the η-η' system to be described reasonably well by a single mixing angle φ and two decay constants fl, fs. The required accuracy of the matrix elements is again guaranteed by the aforementioned, improved analysis method, yielding a value of φ=46.0(0.9)stat(2.7)sys° for the mixing angle extrapolated to the physical point. In addition we obtain results for the ratios flfPS = 0.859(07)stat(64)sys and fs/fK = 1.166(11)stat(31)sys. We find that our data is indeed described well by a single mixing angle, indicating that the η' is mostly a flavor singlet state. Moreover, our results confirm that the charm quark does not contribute to any of the two states within errors. Apart from the flavor singlet sector, we also perform calculations of masses for the remaining light pseudoscalar octet mesons. Matching these masses to two-flavor Wilson chiral perturbation theory allows for a determination of the low energy constants W'6, W'8 and their linear combination c2 which controls the O(a2) mass splitting between charged and neutral pion. We study the dependence of these low energy constants on the number of dynamical quark flavors and for different choices of the lattice action