We study the problem of testing conditional independence for discrete distributions. Specifically, given samples from a discrete random variable (X, Y, Z) on domain [l(1)] x [l(2)] x [n], we want to distinguish, with probability at least 2/3, between the case that X and Y are conditionally independent given Z from the case that (X, Y, Z) is epsilon-far, in l(1)-distance, from every distribution that has this property. Conditional independence is a concept of central importance in probability and statistics with a range of applications in various scientific domains. As such, the statistical task of testing conditional independence has been extensively studied in various forms within the statistics and econometrics communities for nearly a century. Perhaps surprisingly, this problem has not been previously considered in the framework of distribution property testing and in particular no tester with sublinear sample complexity is known, even for the important special case that the domains of X and Y are binary. The main algorithmic result of this work is the first conditional independence tester with sublinear sample complexity for discrete distributions over [l(1)] x [l(2)] x [n]. To complement our upper bounds, we prove information-theoretic lower bounds establishing that the sample complexity of our algorithm is optimal, up to constant factors, for a number of settings. Specifically, for the prototypical setting when l(1), l(2) = O(1), we show that the sample complexity of testing conditional independence (upper bound and matching lower bound) is Theta(max(n(1/2)/epsilon(2), min(n(7/8)/epsilon, n(6/7)/epsilon(8/7)))). To obtain our tester, we employ a variety of tools, including (1) a suitable weighted adaptation of the flattening technique [DK16], and (2) the design and analysis of an optimal (unbiased) estimator for the following statistical problem of independent interest: Given a degree-d polynomial Q: R-n -> R and sample access to a distribution p over [n], estimate Q(p(1), ... , p(n)) up to small additive error. Obtaining tight variance analyses for specific estimators of this form has been a major technical hurdle in distribution testing (see, e.g., [CDVV14]). As an important contribution of this work, we develop a general theory providing tight variance bounds for all such estimators. Our lower bounds, established using the mutual information method, rely on novel constructions of hard instances that may be useful in other settings.