Output of discrete-time LTI system guaranteed to be same form as input?

I know that in the continuous-time context, if I supply a complex exponential input to a Linear Time Invariant system, the output will be of the same form as the input - for example, if the input is $x(t)=Ae^{j \omega_0 t+\phi}$, the output will be of the form $y(t)=|H(j\omega_o)|\centerdot A \centerdot e^{j(H(s)\omega_0 t+\phi +\theta_H)}$ where $H(s)|_{s=j\omega_0} = |H(j\omega_0)|e^{j\theta_H}$ is the transfer function of the system evaluated at the frequency $j\omega_0$. My point is that the system does not modify the input frequency.

Does the same hold true in the discrete-time context? If I supply an input $x[n]=a^n$ to an LTI system, will the output also contain $a^n$? Or is it possible to supply $a^n$ and get the output $b^n +c^n$ where $a\ne b\ne c$ (for example, $x[n]=(1/6)^n$, $y[n]=B(1/2)^n + C(1/3)^n$)?

Edit: I've posted a related question that contains the issue that actually prompted this question...

Answer

Despite the two other answers, I'll add one more, because I have the feeling that your final paragraph has not yet been answered satisfactorily.

A discrete-time LTI system with impulse response $h[n]$ has a transfer function

$$H(z)=\sum_{n=-\infty}^{\infty}h[n]z^{-n}\tag{1}$$

Its output $y[n]$ for a given input sequence $x[n]$ is given by the convolution sum

$$y[n]=\sum_{k=-\infty}^{\infty}h[k]x[n-k]\tag{2}$$

For $x[n]=a^n$ we get from $(2)$

$$y[n]=\sum_{k=-\infty}^{\infty}h[k]a^na^{-k}=a^n\sum_{k=-\infty}^{\infty}h[k]a^{-k}=a^nH(a)\tag{3}$$

where Eq. $(1)$ was used in the last equality. Eq. $(3)$ shows that for $x[n]=a^n$, the output of a discrete-time LTI system is just a weighted version of the input signal. Hence, the sequence $a^n$ is an eigenfunction of a discrete-time LTI system.