The famous results of Komlós, Major and Tusnády (see [15] and [17]) state that it is possible to approximate almost surely the partial sums of size n of i.i.d. centered random variables in L p (p > 2) by a Wiener process with an error term of order o(n 1/p ). Very recently, Berkes, Liu and Wu [3] extended this famous result to partial sums associated with functions of an i.i.d. sequence, provided a condition on a functional dependence measure in L p is satisfied. In this paper, we adapt the method of Berkes, Liu and Wu to partial sums of functions of random iterates. Taking advantage of the Markovian setting, we shall give new dependent conditions, expressed in terms of a natural coupling (in L ∞ or in L 1 ), under which the strong approximation result holds with rate o(n 1/p ). As we shall see our conditions are well adapted to a large variety of models, including left random walks on GL d (R), contracting iterated random functions, autoregressive Lipschitz processes, and some ergodic Markov chains. We also provide some examples showing that our L 1 -coupling condition is in some sense optimal.