In this paper we study Hardy spaces H p,q (R d ), 0 < p, q < ∞, modeled over amalgam spaces (L p , q )(R d ). We characterize H p,q (R d ) by using first order classical Riesz transforms and compositions of first order Riesz transforms depending on the values of the exponents p and q. Also, we describe the distributions in H p,q (R d ) as the boundary values of solutions of harmonic and caloric Cauchy-Riemann systems. We remark that caloric Cauchy-Riemann systems involve fractional derivative in the time variable. Finally we characterize the functions in2010 Mathematics Subject Classification. 42B30, 42B35.according to [24, p. 89-90]is harmonic and f ∈ H p,q (R d ) if, and only if, the maximal function Atomic characterizations of the distributions inH p,q (R d ) were established in [1, Theorems 4.3, 4.4 and 4.6]. Dual spaces of H p,q (R d ) and the boundedness of certain pseudodifferential operators, Calderón-Zygmund operators and Riesz potentials in these Hardy spaces were studied in [2] and [3]. Recently, Sawano, Ho, D. Yang and S. Yang ([23]) have developed a real variable theory of Hardy spaces modeled over a general class of function spaces called ball quasi-Banach functions spaces. A quasi-Banach space X of measurable functions on R d is a ball quasi-Banach function space on R d when the following properties hold:i) f X = 0 implies that f = 0; ii) If |g| ≤ |f | and f, g ∈ X, then g X ≤ f X ; iii) If {f m } ∞ m=1 is an increasing sequence in X, f ∈ X and lim m→∞ f m (x) = f (x), a.e.x ∈ R d , then f m X −→ f X , as m → ∞;iv) For every x ∈ R d and r > 0, 1 B(x,r) ∈ X.