You haven't said whether $K$, $\Sigma$, or $L_{1}$ are sparse or otherwise specially structured (e.g. Toeplitz.) This could be hugely important in finding the most efficient way to do this.
You also haven't given any indication of the size of the matrices involved. Furthermore, you haven't said anything about the larger algorithm that this computation is part of- it may have much more time consuming steps than this computation. This computation may not be a significant part of the overall computation that you're doing.
For this answer, I'll assume that the matrices are dense and have no other easily exploited special structure, and that the matrices are large enough that asymptotic worst case analysis of running time is appropriate.
One simple option is to compute
$M=L_{1}^{-1}$
in $O(n^{2})$ time, and then compute
$K^{-1}=M^{T}M$
in $O(n^{3})$ time.
You can then compute the trace of the product by multiplying each element of $K^{-1}$ times the corresponding element in $\Sigma$ and taking the sum. That takes $O(n^2)$ time.
Although this procedure takes $O(n^{3})$ time, then only order $n^3$ operation is a matrix-matrix multiplication which is typically very fast in practice compared with other operations such as matrix factorizations.
You could compute columns of $K^{-1}\Sigma$ in $O(n^{2})$ time (by forward and backward solves with the Cholesky factor $L_{1}$.) However, you'd need to perform $2n$ forward or backward solves, so the whole process would still take $O(n^3)$ time, and the constant hidden in the $O(n^3)$ would likely be larger in practice.