The authors study the function \[ \zeta(k_1,\dots, k_{n-1};s)= \sum_{0< m_1< m_2<\dots< m_n} \frac{1} {m_1^{k_1}\cdots m_{n-1}^{k_{n-1}} m_n^s} \] and show that the poly-Bernoulli numbers introduced in a previous paper [M. Kaneko, J. Théor. Nombres Bordx. 9, 221–228 (1997; Zbl 0887.11011)] are expressed as special values at negative arguments of certain combinations of these functions. As a consequence of their study, the authors obtain a series of relations among multiple zeta values.