The aim of this paper is twofold. Firstly, we derive upper and lower non-Gaussian bounds for the densities of the marginal laws of the solutions to backward stochastic differential equations (BSDEs) driven by fractional Brownian motions. Our arguments consist of utilising a relationship between fractional BSDEs and quasilinear partial differential equations of mixed type, together with the Nourdin-Viens formula. In the linear case, upper and lower Gaussian bounds for the densities and the tail probabilities of solutions are obtained with simple arguments by their explicit expressions in terms of the quasi-conditional expectation. Secondly, we are concerned with Gaussian estimates for the densities of a BSDE driven by a Gaussian process in the manner that the solution can be established via an auxiliary BSDE driven by a Brownian motion. Using the transfer theorem we succeed in deriving Gaussian estimates for the solutions. In addition, a representation theorem for this setting is also obtained.On the other hand, in the seminal paper [34] Pardoux and Peng initiated the theory of nonlinear backward stochastic differential equations (BSDEs), which is of increasing importance in stochastic control and mathematical finance (see, e.g., [14] and most recently [37]). This class of equations is of the following formis a jointly measurable, and B = (B t ) t≥0 is a Brownian motion adapted to {F t } t∈[0,T ] or simply {F t } t∈[0,T ] is taken as the natural filtration of B. Recall that a solution to the BSDE (1.1) is a pair of predictable processes (y, z) with suitable integrability conditions such that (1.1) holds P-a.s.. To date, there is a wealth of existence and uniqueness results under various assumptions on the generators f including e.g. the cases of Lipschitz or (super-)quadratic growth [13,14,23,34,35]. When dealing with applications, one needs to investigate the existence and regularity of densities for the marginal laws of (y, z). As far as we know, there are comparably only a few works to study this problem. The first results have been derived by Antonelli and Kohatsu-Higa [3], in which they study the existence and the estimates of the density for y t at a fixed time t ∈ [0, T ] via the Bouleau-Hirsch criterion. Then, based on the Nourdin-Viens formula, Aboura and Bourguin [1] have proved the existence of the density for z t under the condition that the generator f is linear w.r.t. its z variable, and moreover obtained the estimates on the densities of the laws of y t and z t . Recently, Mastrolia, Possamaï and Réveillac in [26] have studied the existence of densities for marginal laws of the solution (y, z) to (1.1) with a quadratic growth generator, and derived the estimates on these densities. Afterwards, Mastrolia [25] has extended the results to the case of non-Markovian BSDEs.One of the main objective of the present paper concerns the problem of density estimates for the following BSDEwhere η t = η 0 + t 0 b s ds + t 0 σ s dB H s with η 0 , b s and σ s being respectively a constant and deterministic function...